diff --git a/core/geo_utils.py b/core/geo_utils.py index 7d69c81..ce0cbf2 100644 --- a/core/geo_utils.py +++ b/core/geo_utils.py @@ -4,6 +4,7 @@ from math import floor import geopandas from pyproj import Transformer +from shapely import prepare from shapely.geometry import Point, Polygon TRANSFORMER_OS_GRID_TO_WGS84 = Transformer.from_crs("EPSG:27700", "EPSG:4326") @@ -12,15 +13,28 @@ TRANSFORMER_CI_UTM_GRID_TO_WGS84 = Transformer.from_crs("+proj=utm +zone=30 +ell cq_zone_data = geopandas.GeoDataFrame.from_features(geopandas.read_file("datafiles/cqzones.geojson")) itu_zone_data = geopandas.GeoDataFrame.from_features(geopandas.read_file("datafiles/ituzones.geojson")) - +for idx in cq_zone_data.index: + prepare(cq_zone_data.at[idx, 'geometry']) +for idx in itu_zone_data.index: + prepare(itu_zone_data.at[idx, 'geometry']) # Finds out which CQ zone a lat/lon point is in. def lat_lon_to_cq_zone(lat, lon): + lon = ((lon + 180) % 360) - 180 for index, row in cq_zone_data.iterrows(): polygon = Polygon(row["geometry"]) test_point = Point(lon, lat) if polygon.contains(test_point): return int(row["name"]) + + # Might have problems around the antemeridian, so if we didn't find a match, try offsetting the point by + or - + # 360 degrees longitude to try the other side of the Earth + if lon < 0: + test_point = Point(lon + 360, lat) + else: + test_point = Point(lon - 360, lat) + if polygon.contains(test_point): + return int(row["name"]) return None