Skip to content

Commit

Permalink
add get_bbox
Browse files Browse the repository at this point in the history
  • Loading branch information
akorosov committed Jan 10, 2022
1 parent 2c53cde commit 06a8b88
Show file tree
Hide file tree
Showing 2 changed files with 36 additions and 6 deletions.
19 changes: 16 additions & 3 deletions geodataset/geodataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -359,21 +359,34 @@ def get_lonlat_arrays(self):
# if lon and lat are arrays
return [self.get_variable_array(name) for name in self.lonlat_names]

def get_area_euclidean(self, pyproj_map):
def get_area_euclidean(self, mapping):
"""
Calculates element area from netcdf file
Assumes regular grid in the given projection
Parameters:
-----------
pyproj_map : pyproj.Proj
mapping : pyproj.Proj
Returns:
--------
* area (float)
"""
lon, lat = self.get_lonlat_arrays()
x, y = pyproj_map(lon, lat)
x, y = mapping(lon, lat)
dy = np.max([np.abs(np.mean(y[:, 2]-y[:, 1])), np.abs(np.mean(y[1, :]-y[0, :]))])
dx = np.max([np.abs(np.mean(x[:, 2]-x[:, 1])), np.abs(np.mean(x[1, :]-x[0, :]))])
return np.abs(dx*dy)

def get_bbox(self, mapping):
""" Get bounding box (extent)
Parameters:
* mapping: pyproj mapping
Returns:
* bbox = [xmin, xmax, ymin, ymax], where x,y are coordinates specified by mapping
"""

lon, lat = self.get_lonlat_arrays()
x, y = mapping(lon, lat)
return [x.min(), x.max(), y.min(), y.max()]
23 changes: 20 additions & 3 deletions geodataset/tests/test_geodataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -334,17 +334,34 @@ def test_projection(self):
__exit__=MagicMock(return_value=None),
get_lonlat_arrays=DEFAULT,
)
def test_get_lonlat_names_raises(self, **kwargs):
def test_get_area_euclidean(self, **kwargs):
''' test f4 with _FillValue defined '''
p = pyproj.Proj(3411)

kwargs['get_lonlat_arrays'].return_value = (
GeoDatasetRead.get_lonlat_arrays.return_value = (
np.array([[1,2,3],[1,2,3],[1,2,3]]),
np.array([[1,1,1],[2,2,2],[3,3,3]]))

with GeoDatasetRead() as ds:
area = ds.get_area_euclidean(p)
self.assertAlmostEqual(area, 23354252971.32609)

@patch.multiple(GeoDatasetRead,
__init__=MagicMock(return_value=None),
__exit__=MagicMock(return_value=None),
get_lonlat_arrays=DEFAULT,
)
def test_get_bbox(self, **kwargs):
''' test f4 with _FillValue defined '''
p = pyproj.Proj(3411)
GeoDatasetRead.get_lonlat_arrays.return_value = (
np.array([[1,2,3],[1,2,3],[1,2,3]]),
np.array([[1,1,1],[2,2,2],[3,3,3]]))

with GeoDatasetRead() as ds:
bbox = ds.get_bbox(p)
self.assertAlmostEqual(bbox,
[8420199.606917838, 9005961.652806347,
-8418368.037664523, -7832478.150085783])

if __name__ == "__main__":
unittest.main()

0 comments on commit 06a8b88

Please sign in to comment.