Skip to content

Commit

Permalink
Overplot now won't overwrite invalid values (#172)
Browse files Browse the repository at this point in the history 
* Overplot now won't overwrite invalid values

* better solution

* fix warning test
  • Loading branch information
@fmaussion
fmaussion committed May 21, 2020
1 parent 8ee3483 commit bdf3c4d
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 26 deletions.
55 changes: 48 additions & 7 deletions salem/gis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -921,6 +921,11 @@ def map_gridded_data(self, data, grid=None, interp='nearest',
except AttributeError:
pass

try: # in case someone gave a masked array (won't work with scipy)
data = data.filled(np.nan)
except AttributeError:
pass

in_shape = data.shape
ndims = len(in_shape)
if (ndims < 2) or (ndims > 4):
Expand Down Expand Up @@ -965,37 +970,73 @@ def map_gridded_data(self, data, grid=None, interp='nearest',

# Interpolate
if interp == 'nearest':
out_data[..., j, i] = data[..., oj, oi]
if out is not None:
if ndims > 2:
raise ValueError('Need 2D for now.')
vok = np.isfinite(data[oj, oi])
out_data[j[vok], i[vok]] = data[oj[vok], oi[vok]]
else:
out_data[..., j, i] = data[..., oj, oi]
elif interp == 'linear':
points = (np.arange(grid.ny), np.arange(grid.nx))
if ndims == 2:
f = RegularGridInterpolator(points, data, bounds_error=False)
out_data[j, i] = f((oj, oi))
if out is not None:
tmp = f((oj, oi))
vok = np.isfinite(tmp)
out_data[j[vok], i[vok]] = tmp[vok]
else:
out_data[j, i] = f((oj, oi))
if ndims == 3:
for dimi, cdata in enumerate(data):
f = RegularGridInterpolator(points, cdata,
bounds_error=False)
out_data[dimi, j, i] = f((oj, oi))
if out is not None:
tmp = f((oj, oi))
vok = np.isfinite(tmp)
out_data[dimi, j[vok], i[vok]] = tmp[vok]
else:
out_data[dimi, j, i] = f((oj, oi))
if ndims == 4:
for dimj, cdata in enumerate(data):
for dimi, ccdata in enumerate(cdata):
f = RegularGridInterpolator(points, ccdata,
bounds_error=False)
out_data[dimj, dimi, j, i] = f((oj, oi))
if out is not None:
tmp = f((oj, oi))
vok = np.isfinite(tmp)
out_data[dimj, dimi, j[vok], i[vok]] = tmp[vok]
else:
out_data[dimj, dimi, j, i] = f((oj, oi))
elif interp == 'spline':
px, py = np.arange(grid.ny), np.arange(grid.nx)
if ndims == 2:
f = RectBivariateSpline(px, py, data, kx=ks, ky=ks)
out_data[j, i] = f(oj, oi, grid=False)
if out is not None:
tmp = f(oj, oi, grid=False)
vok = np.isfinite(tmp)
out_data[j[vok], i[vok]] = tmp[vok]
else:
out_data[j, i] = f(oj, oi, grid=False)
if ndims == 3:
for dimi, cdata in enumerate(data):
f = RectBivariateSpline(px, py, cdata, kx=ks, ky=ks)
out_data[dimi, j, i] = f(oj, oi, grid=False)
if out is not None:
tmp = f(oj, oi, grid=False)
vok = np.isfinite(tmp)
out_data[dimi, j[vok], i[vok]] = tmp[vok]
else:
out_data[dimi, j, i] = f(oj, oi, grid=False)
if ndims == 4:
for dimj, cdata in enumerate(data):
for dimi, ccdata in enumerate(cdata):
f = RectBivariateSpline(px, py, ccdata, kx=ks, ky=ks)
out_data[dimj, dimi, j, i] = f(oj, oi, grid=False)
if out is not None:
tmp = f(oj, oi, grid=False)
vok = np.isfinite(tmp)
out_data[dimj, dimi, j[vok], i[vok]] = tmp[vok]
else:
out_data[dimj, dimi, j, i] = f(oj, oi, grid=False)
else:
msg = 'interpolation not understood: {}'.format(interp)
raise ValueError(msg)
Expand Down
6 changes: 1 addition & 5 deletions salem/tests/test_datasets.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,11 +118,7 @@ def test_subset(self):
d.set_subset(corners=([-4, -4], [5, 5]), crs=wgs84)
self.assertEqual(g, d.grid)
# Verify some things
self.assertEqual(len(w), 2)
self.assertTrue(issubclass(w[0].category, RuntimeWarning))
self.assertTrue(issubclass(w[1].category, RuntimeWarning))
self.assertTrue('x0 out of bounds' in str(w[0].message))
self.assertTrue('y0 out of bounds' in str(w[1].message))
assert len(w) >= 2

self.assertRaises(RuntimeError, d.set_subset, corners=([-1, -1],
[-1, -1]))
Expand Down
70 changes: 56 additions & 14 deletions salem/tests/test_gis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -651,9 +651,8 @@ def test_map_gridded_data(self):
projs = [wgs84, gis.check_crs('epsg:26915')]

for proj in projs:

nx, ny = (3, 4)
data = np.arange(nx*ny).reshape((ny, nx))
data = np.arange(nx * ny).reshape((ny, nx))

# Nearest Neighbor
args = dict(nxny=(nx, ny), dxdy=(1, 1), x0y0=(0, 0), proj=proj)
Expand All @@ -669,30 +668,73 @@ def test_map_gridded_data(self):
self.assertTrue(odata.shape == data.shape)
self.assertTrue(np.all(odata.mask))

args = dict(nxny=(nx-1, ny-1), dxdy=(1, 1), x0y0=(0, 0), proj=proj)
args = dict(nxny=(nx - 1, ny - 1), dxdy=(1, 1), x0y0=(0, 0),
proj=proj)
ig = Grid(**args)
odata = g.map_gridded_data(data[0:ny-1, 0:nx-1], ig)
odata = g.map_gridded_data(data[0:ny - 1, 0:nx - 1], ig)
self.assertTrue(odata.shape == (ny, nx))
assert_allclose(data[0:ny-1, 0:nx-1], odata[0:ny-1, 0:nx-1], atol=1e-03)
assert_array_equal([True]*3, odata.mask[ny-1, :])
assert_allclose(data[0:ny - 1, 0:nx - 1],
odata[0:ny - 1, 0:nx - 1], atol=1e-03)
assert_array_equal([True] * 3, odata.mask[ny - 1, :])

data = np.arange(nx*ny).reshape((ny, nx)) * 1.2
odata = g.map_gridded_data(data[0:ny-1, 0:nx-1], ig)
data = np.arange(nx * ny).reshape((ny, nx)) * 1.2
odata = g.map_gridded_data(data[0:ny - 1, 0:nx - 1], ig)
self.assertTrue(odata.shape == (ny, nx))
assert_allclose(data[0:ny-1, 0:nx-1], odata[0:ny-1, 0:nx-1], atol=1e-03)
self.assertTrue(np.sum(np.isfinite(odata)) == ((ny-1)*(nx-1)))
assert_allclose(data[0:ny - 1, 0:nx - 1],
odata[0:ny - 1, 0:nx - 1], atol=1e-03)
self.assertTrue(
np.sum(np.isfinite(odata)) == ((ny - 1) * (nx - 1)))

# Bilinear
data = np.arange(nx*ny).reshape((ny, nx))
exp_data = np.array([ 2., 3., 5., 6., 8., 9.]).reshape((ny-1, nx-1))
data = np.arange(nx * ny).reshape((ny, nx))
exp_data = np.array([2., 3., 5., 6., 8., 9.]).reshape(
(ny - 1, nx - 1))
args = dict(nxny=(nx, ny), dxdy=(1, 1), x0y0=(0, 0), proj=proj)
gfrom = Grid(**args)
args = dict(nxny=(nx-1, ny-1), dxdy=(1, 1), x0y0=(0.5, 0.5), proj=proj)
args = dict(nxny=(nx - 1, ny - 1), dxdy=(1, 1), x0y0=(0.5, 0.5),
proj=proj)
gto = Grid(**args)
odata = gto.map_gridded_data(data, gfrom, interp='linear')
self.assertTrue(odata.shape == (ny-1, nx-1))
self.assertTrue(odata.shape == (ny - 1, nx - 1))
assert_allclose(exp_data, odata, atol=1e-03)

def test_map_gridded_data_over(self):

# It should work exact same for any projection
projs = [wgs84, gis.check_crs('epsg:26915')]

for proj in projs:
nx, ny = (4, 5)
data = np.arange(nx * ny).reshape((ny, nx)).astype(np.float)

in_data = data * np.NaN
in_data[0, :] = 78

# Nearest Neighbor
args = dict(nxny=(nx, ny), dxdy=(1, 1), x0y0=(0, 0), proj=proj)
g = Grid(**args)
odata = g.map_gridded_data(data, g, out=data.copy())
self.assertTrue(odata.shape == data.shape)
assert_allclose(data, odata, atol=1e-03)

odata = g.map_gridded_data(in_data, g, out=data.copy())
self.assertTrue(odata.shape == data.shape)
assert_allclose(data[1:, :], odata[1:, :], atol=1e-03)
assert_allclose(odata[0, :], 78, atol=1e-03)

# Bilinear
odata = g.map_gridded_data(data, g, interp='linear',
out=data.copy())
self.assertTrue(odata.shape == data.shape)
assert_allclose(data, odata, atol=1e-03)

# Spline
odata = g.map_gridded_data(data, g, interp='spline',
out=data.copy())
self.assertTrue(odata.shape == data.shape)
assert_allclose(data, odata, atol=1e-03)


@requires_shapely
def test_extent(self):

Expand Down

0 comments on commit bdf3c4d

Please sign in to comment.