w1cdn/spothole
1
0
Fork 0
mirror of https://git.ianrenton.com/ian/spothole.git synced 2026-03-15 12:24:29 +00:00
Code Issues 13 Packages Projects Releases Wiki Activity

Fix CQ/ITU lookups for zones that cross the antemeridian

Browse Source
This commit is contained in:
Ian Renton
2026-02-14 07:46:16 +00:00
parent 76f289d66e
commit 4c22861666
1 changed files with 15 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
core/geo_utils.py
View File

@@ -4,6 +4,7 @@ from math import floor
import geopandas import geopandas
from pyproj import Transformer from pyproj import Transformer
from shapely import prepare
from shapely.geometry import Point, Polygon from shapely.geometry import Point, Polygon
TRANSFORMER_OS_GRID_TO_WGS84 = Transformer.from_crs("EPSG:27700", "EPSG:4326") 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")) 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")) 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. # Finds out which CQ zone a lat/lon point is in.
def lat_lon_to_cq_zone(lat, lon): def lat_lon_to_cq_zone(lat, lon):
lon = ((lon + 180) % 360) - 180
for index, row in cq_zone_data.iterrows(): for index, row in cq_zone_data.iterrows():
polygon = Polygon(row["geometry"]) polygon = Polygon(row["geometry"])
test_point = Point(lon, lat) test_point = Point(lon, lat)
if polygon.contains(test_point): if polygon.contains(test_point):
return int(row["name"]) 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 return None