[SciPy-User] calculating 3D surface elevation volume difference

Chimmula, Sumani chimmulas at usgs.gov
Wed Dec 6 14:30:50 EST 2017


I'm trying to calculate difference between two interpolated grids. The
grids are irregular and of different shapes, so cannot use np.subtract
to calculate the difference. Any help on how to calculate the
difference?

Another important information is these are spatial grids with
overlapping spatial extents. If we just broadcast/extend array shapes,
would that be enough? I do not want to subtract the grids at wrong x/y
locations. So basically like in GIS/Arcmap how we calculate raster
difference at each cell at a particular xy cell location.

#create two interpolated grids from two set of xyz points and
calculate difference
ixx = np.array(self.ix)
    iyy = np.array(self.iy)
    izz = np.array(self.iz)

    iminx = min(ixx.flat)
    imaxx = max(ixx.flat)
    iminy = min(iyy.flat)
    imaxy = max(iyy.flat)

    # post surface - get extents
    pxx = np.array(self.px)
    pyy = np.array(self.py)
    pzz = np.array(self.pz)

    pminx = min(pxx.flat)
    pmaxx = max(pxx.flat)
    pminy = min(pyy.flat)
    pmaxy = max(pyy.flat)

    ixi = np.arange(iminx, imaxx, self.cellsize)
    iyi = np.arange(iminy, imaxy, self.cellsize)
    iX,iY = np.meshgrid(ixi,iyi)
    iextent = (min(ixi), max(ixi), min(iyi), max(iyi))

    pxi = np.arange(pminx, pmaxx, self.cellsize)
    pyi = np.arange(pminy, pmaxy, self.cellsize)
    pX,pY = np.meshgrid(pxi,pyi)
    pextent = (min(pxi), max(pxi), min(pyi), max(pyi))

    igrid = griddata((ixx, iyy), izz, (iX, iY), method='linear')
    pgrid = griddata((pxx, pyy), pzz, (pX, pY), method='linear')


# would like to calcualte difference, but this wont work, due to
different shapes
changegrid = np.subtract(pgrid, igrid)

igrid shape is (181, 130), pgrid shape is (184, 126).

Thanks.


More information about the SciPy-User mailing list