Skip to content

Commit

Permalink
Bugfix/get geometries extent antimeridian (#524)
Browse files Browse the repository at this point in the history 
* Ensure get_country_geometries deals correctly with extent crossing the anti-meridian

* get_land_gemoetry to rely on get_country_geometries to avoid redundancies

* Ensure coord_on_land properly deals with extent of land_geom and lon

* Adding a simple unit test for coord_on_land with points beyond the anti-meridian

* Fix in case land_geom is empty (no land within extent)

* provide keyword arg in unit test test_get_land_geometry_country_pass

* u_coord: some refactoring and linting

* u_coord: fix tests for get_land_geometry

* coordinates: shapely.geometry.Polygon still referenced
in `points_to_raster`

* test_coordinates: complementary testcases for get_country_geometries

Co-authored-by: Thomas Vogt <[email protected]>
Co-authored-by: emanuel-schmid <[email protected]>
  • Loading branch information
@bguillod @emanuel-schmid
3 people authored Aug 5, 2022
1 parent 8a7e369 commit eed7058
Show file tree
Hide file tree
Showing 3 changed files with 114 additions and 51 deletions.
26 changes: 26 additions & 0 deletions climada/hazard/test/test_tc_tracks.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -897,6 +897,32 @@ def test_get_landfall_idx(self):
self.assertEqual(sea_land_idx.tolist(), [0,4])
self.assertEqual(land_sea_idx.tolist(), [2,6])

def test_track_land_params(self):
"""Test identification of points on land and distance since landfall"""
# 1 pt over the ocean and two points within Fiji, one on each side of the anti-meridian
lon_test = np.array([170, 179.18, 180.05])
lat_test = np.array([-60, -16.56, -16.85])
on_land = np.array([False, True, True])
lon_shift = np.array([-360, 0, 360])
# ensure both points are considered on land as is
np.testing.assert_array_equal(
u_coord.coord_on_land(lat = lat_test, lon = lon_test),
on_land
)
# independently on shifts by 360 degrees in longitude
np.testing.assert_array_equal(
u_coord.coord_on_land(lat = lat_test, lon = lon_test + lon_shift),
on_land
)
np.testing.assert_array_equal(
u_coord.coord_on_land(lat = lat_test, lon = lon_test - lon_shift),
on_land
)
# also when longitude is within correct range
np.testing.assert_array_equal(
u_coord.coord_on_land(lat = lat_test, lon = u_coord.lon_normalize(lon_test)),
on_land
)

# Execute Tests
if __name__ == "__main__":
Expand Down
105 changes: 65 additions & 40 deletions climada/util/coordinates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -647,33 +647,9 @@ def get_land_geometry(country_names=None, extent=None, resolution=10):
geom : shapely.geometry.multipolygon.MultiPolygon
Polygonal shape of union.
"""
resolution = nat_earth_resolution(resolution)
shp_file = shapereader.natural_earth(resolution=resolution,
category='cultural',
name='admin_0_countries')
reader = shapereader.Reader(shp_file)
if (country_names is None) and (extent is None):
LOGGER.info("Computing earth's land geometry ...")
geom = list(reader.geometries())
geom = shapely.ops.unary_union(geom)

elif country_names:
countries = list(reader.records())
geom = [country.geometry for country in countries
if (country.attributes['ISO_A3'] in country_names) or
(country.attributes['WB_A3'] in country_names) or
(country.attributes['ADM0_A3'] in country_names)]
geom = shapely.ops.unary_union(geom)

else:
extent_poly = Polygon([(extent[0], extent[2]), (extent[0], extent[3]),
(extent[1], extent[3]), (extent[1], extent[2])])
geom = []
for cntry_geom in reader.geometries():
inter_poly = cntry_geom.intersection(extent_poly)
if not inter_poly.is_empty:
geom.append(inter_poly)
geom = shapely.ops.unary_union(geom)
geom = get_country_geometries(country_names, extent, resolution)
# combine all into a single multipolygon
geom = geom.geometry.unary_union
if not isinstance(geom, MultiPolygon):
geom = MultiPolygon([geom])
return geom
Expand Down Expand Up @@ -701,14 +677,29 @@ def coord_on_land(lat, lon, land_geom=None):
if lat.size == 0:
return np.empty((0,), dtype=bool)
delta_deg = 1
lons = lon.copy()
if land_geom is None:
# ensure extent of longitude is consistent
bounds = lon_bounds(lons)
lon_mid = 0.5 * (bounds[0] + bounds[1])
# normalize lon
lon_normalize(lons, center=lon_mid)
# load land geometry with appropriate same extent
land_geom = get_land_geometry(
extent=(np.min(lon) - delta_deg,
np.max(lon) + delta_deg,
extent=(bounds[0] - delta_deg,
bounds[1] + delta_deg,
np.min(lat) - delta_deg,
np.max(lat) + delta_deg),
resolution=10)
return shapely.vectorized.contains(land_geom, lon, lat)
elif not land_geom.is_empty:
# ensure lon values are within extent of provided land_geom
land_bounds = land_geom.bounds
if lons.max() > land_bounds[2] or lons.min() < land_bounds[0]:
# normalize longitude to land_geom extent
lon_mid = 0.5 * (land_bounds[0] + land_bounds[2])
lon_normalize(lons, center=lon_mid)

return shapely.vectorized.contains(land_geom, lons, lat)

def nat_earth_resolution(resolution):
"""Check if resolution is available in Natural Earth. Build string.
Expand Down Expand Up @@ -741,11 +732,16 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
starts including the projection information. (They are saving a whopping 147 bytes by omitting
it.) Same goes for UTF.
If extent is provided, longitude values in 'geom' will all lie within 'extent' longitude
range. Therefore setting extent to e.g. [160, 200, -20, 20] will provide longitude values
between 160 and 200 degrees.
Parameters
----------
country_names : list, optional
list with ISO 3166 alpha-3 codes of countries, e.g ['ZWE', 'GBR', 'VNM', 'UZB']
extent : tuple (min_lon, max_lon, min_lat, max_lat), optional
extent : tuple, optional
(min_lon, max_lon, min_lat, max_lat)
Extent, assumed to be in the same CRS as the natural earth data.
resolution : float, optional
10, 50 or 110. Resolution in m. Default: 10m
Expand Down Expand Up @@ -777,18 +773,47 @@ def get_country_geometries(country_names=None, extent=None, resolution=10):
if country_names:
if isinstance(country_names, str):
country_names = [country_names]
out = out[out.ISO_A3.isin(country_names)]
country_mask = np.isin(
nat_earth[['ISO_A3', 'WB_A3', 'ADM0_A3']].values,
country_names,
).any(axis=1)
out = out[country_mask]

if extent:
bbox = Polygon([
(extent[0], extent[2]),
(extent[0], extent[3]),
(extent[1], extent[3]),
(extent[1], extent[2])
])
bbox = gpd.GeoSeries(bbox, crs=out.crs)
bbox = gpd.GeoDataFrame({'geometry': bbox}, crs=out.crs)
if extent[1] - extent[0] > 360:
raise ValueError(
f"longitude extent range is greater than 360: {extent[0]} to {extent[1]}"
)

if extent[1] < extent[0]:
raise ValueError(
f"longitude extent at the left ({extent[0]}) is larger "
f"than longitude extent at the right ({extent[1]})"
)

# rewrap longitudes unless longitude extent is already normalized (within [-180, +180])
lon_normalized = extent[0] >= -180 and extent[1] <= 180
if lon_normalized:
bbox = box(extent[0], extent[2], extent[1], extent[3])
else:
# split the extent box into two boxes both within [-180, +180] in longitude
lon_left, lon_right = lon_normalize(np.array(extent[:2]))
extent_left = (lon_left, 180, extent[2], extent[3])
extent_right = (-180, lon_right, extent[2], extent[3])
bbox = shapely.ops.unary_union(
[box(e[0], e[2], e[1], e[3]) for e in [extent_left, extent_right]]
)
bbox = gpd.GeoSeries(bbox, crs=DEF_CRS)
bbox = gpd.GeoDataFrame({'geometry': bbox}, crs=DEF_CRS)
out = gpd.overlay(out, bbox, how="intersection")
if ~lon_normalized:
lon_mid = 0.5 * (extent[0] + extent[1])
# reset the CRS attribute after rewrapping (we don't really change the CRS)
out = (
out
.to_crs({"proj": "longlat", "lon_wrap": lon_mid})
.set_crs(DEF_CRS, allow_override=True)
)

return out

Expand Down
34 changes: 23 additions & 11 deletions climada/util/test/test_coordinates.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -135,25 +135,25 @@ def test_dist_sqr_approx_pass(self):
self.assertAlmostEqual(
7709.827814738594,
np.sqrt(u_coord._dist_sqr_approx(lats1, lons1, cos_lats1, lats2, lons2)) * ONE_LAT_KM)

def test_geodesic_length_geog(self):
"""Test compute_geodesic_lengths for geographic input crs"""

LINE_PATH = DEMO_DIR.joinpath('nl_rails.gpkg')
gdf_rails = gpd.read_file(LINE_PATH).to_crs('epsg:4326')
lengths_geom = u_coord.compute_geodesic_lengths(gdf_rails)

self.assertEqual(len(lengths_geom), len(gdf_rails))
self.assertTrue(
np.all(
(abs(lengths_geom - gdf_rails['distance'])/lengths_geom < 0.1) |
(abs(lengths_geom - gdf_rails['distance'])/lengths_geom < 0.1) |
(lengths_geom - gdf_rails['distance'] < 10)
)
)

def test_geodesic_length_proj(self):
"""Test compute_geodesic_lengths for projected input crs"""

LINE_PATH = DEMO_DIR.joinpath('nl_rails.gpkg')
gdf_rails = gpd.read_file(LINE_PATH).to_crs('epsg:4326')
gdf_rails_proj = gpd.read_file(LINE_PATH).to_crs('epsg:4326').to_crs('EPSG:28992')
Expand All @@ -163,10 +163,10 @@ def test_geodesic_length_proj(self):

for len_proj, len_geom in zip(lengths_proj,lengths_geom):
self.assertAlmostEqual(len_proj, len_geom, 1)

self.assertTrue(
np.all(
(abs(lengths_proj - gdf_rails_proj['distance'])/lengths_proj < 0.1) |
(abs(lengths_proj - gdf_rails_proj['distance'])/lengths_proj < 0.1) |
(lengths_proj - gdf_rails_proj['distance'] < 10)
)
)
Expand Down Expand Up @@ -859,21 +859,21 @@ def test_get_coastlines_pass(self):
def test_get_land_geometry_country_pass(self):
"""get_land_geometry with selected countries."""
iso_countries = ['DEU', 'VNM']
res = u_coord.get_land_geometry(iso_countries, 110)
res = u_coord.get_land_geometry(country_names=iso_countries, resolution=10)
self.assertIsInstance(res, shapely.geometry.multipolygon.MultiPolygon)
for res, ref in zip(res.bounds, (5.85248986800, 8.56557851800,
109.47242272200, 55.065334377000)):
self.assertAlmostEqual(res, ref)

iso_countries = ['ESP']
res = u_coord.get_land_geometry(iso_countries, 110)
res = u_coord.get_land_geometry(country_names=iso_countries, resolution=10)
self.assertIsInstance(res, shapely.geometry.multipolygon.MultiPolygon)
for res, ref in zip(res.bounds, (-18.16722571499986, 27.642238674000,
4.337087436000, 43.793443101)):
self.assertAlmostEqual(res, ref)

iso_countries = ['FRA']
res = u_coord.get_land_geometry(iso_countries, 110)
res = u_coord.get_land_geometry(country_names=iso_countries, resolution=10)
self.assertIsInstance(res, shapely.geometry.multipolygon.MultiPolygon)
for res, ref in zip(res.bounds, (-61.79784094999991, -21.37078215899993,
55.854502800000034, 51.08754088371883)):
Expand Down Expand Up @@ -958,7 +958,7 @@ def test_get_country_geometries_country_pass(self):

def test_get_country_geometries_country_norway_pass(self):
"""test correct numeric ISO3 for country Norway"""
iso_countries = ['NOR']
iso_countries = 'NOR'
extent = [10, 11, 55, 60]
res1 = u_coord.get_country_geometries(iso_countries)
res2 = u_coord.get_country_geometries(extent=extent)
Expand Down Expand Up @@ -1003,6 +1003,18 @@ def test_get_country_geometries_all_pass(self):
self.assertIsInstance(res, gpd.geodataframe.GeoDataFrame)
self.assertAlmostEqual(res.area[0], 1.639510995900778)

def test_get_country_geometries_fail(self):
"""get_country_geometries with offensive parameters"""
with self.assertRaises(ValueError) as cm:
u_coord.get_country_geometries(extent=(-20,350,0,0))
self.assertIn("longitude extent range is greater than 360: -20 to 350",
str(cm.exception))
with self.assertRaises(ValueError) as cm:
u_coord.get_country_geometries(extent=(350,-20,0,0))
self.assertIn("longitude extent at the left (350) is larger "
"than longitude extent at the right (-20)",
str(cm.exception))

def test_country_code_pass(self):
"""Test set_region_id"""

Expand Down

0 comments on commit eed7058

Please sign in to comment.