mirror of
https://git.ianrenton.com/ian/spothole.git
synced 2026-04-30 10:45:57 +00:00
267 lines
11 KiB
Python
267 lines
11 KiB
Python
import json
|
|
import logging
|
|
import re
|
|
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")
|
|
TRANSFORMER_IRISH_GRID_TO_WGS84 = Transformer.from_crs("EPSG:29903", "EPSG:4326")
|
|
TRANSFORMER_CI_UTM_GRID_TO_WGS84 = Transformer.from_crs("+proj=utm +zone=30 +ellps=WGS84", "EPSG:4326")
|
|
|
|
with open("datafiles/cqzones.geojson") as f:
|
|
cq_zone_data = geopandas.GeoDataFrame.from_features(json.load(f)["features"])
|
|
with open("datafiles/ituzones.geojson") as f:
|
|
itu_zone_data = geopandas.GeoDataFrame.from_features(json.load(f)["features"])
|
|
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'])
|
|
|
|
|
|
def lat_lon_to_cq_zone(lat, lon):
|
|
"""Finds out which CQ zone a lat/lon point is in."""
|
|
|
|
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
|
|
|
|
|
|
def lat_lon_to_itu_zone(lat, lon):
|
|
"""Finds out which ITU zone a lat/lon point is in."""
|
|
|
|
lon = ((lon + 180) % 360) - 180
|
|
for index, row in itu_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
|
|
|
|
|
|
def lat_lon_for_grid_centre(grid):
|
|
"""Convert a Maidenhead grid reference of arbitrary precision to the lat/long of the centre point of the square.
|
|
Returns None if the grid format is invalid."""
|
|
|
|
lat, lon, lat_cell_size, lon_cell_size = lat_lon_for_grid_sw_corner_plus_size(grid)
|
|
if lat is not None and lon is not None and lat_cell_size is not None and lon_cell_size is not None:
|
|
return [lat + lat_cell_size / 2.0, lon + lon_cell_size / 2.0]
|
|
else:
|
|
return None
|
|
|
|
|
|
def lat_lon_for_grid_sw_corner(grid):
|
|
"""Convert a Maidenhead grid reference of arbitrary precision to the lat/long of the southwest corner of the square.
|
|
Returns None if the grid format is invalid."""
|
|
|
|
lat, lon, lat_cell_size, lon_cell_size = lat_lon_for_grid_sw_corner_plus_size(grid)
|
|
if lat is not None and lon is not None:
|
|
return [lat, lon]
|
|
else:
|
|
return None
|
|
|
|
|
|
def lat_lon_for_grid_ne_corner(grid):
|
|
"""Convert a Maidenhead grid reference of arbitrary precision to the lat/long of the northeast corner of the square.
|
|
Returns None if the grid format is invalid."""
|
|
|
|
lat, lon, lat_cell_size, lon_cell_size = lat_lon_for_grid_sw_corner_plus_size(grid)
|
|
if lat is not None and lon is not None and lat_cell_size is not None and lon_cell_size is not None:
|
|
return [lat + lat_cell_size, lon + lon_cell_size]
|
|
else:
|
|
return None
|
|
|
|
|
|
def lat_lon_for_grid_sw_corner_plus_size(grid):
|
|
"""Convert a Maidenhead grid reference of arbitrary precision to lat/long, including in the result the size of the
|
|
lowest grid square. This is a utility method used by the main methods that return the centre, southwest, and
|
|
northeast coordinates of a grid square.
|
|
The return type is always a tuple of size 4. The elements in it are None if the grid format is invalid."""
|
|
|
|
# Make sure we are in upper case so our maths works. Case is arbitrary for Maidenhead references
|
|
grid = grid.upper()
|
|
|
|
# Return None if our Maidenhead string is invalid or too short
|
|
length = len(grid)
|
|
if length <= 0 or (length % 2) != 0:
|
|
return None, None, None, None
|
|
|
|
lat = 0.0 # aggregated latitude
|
|
lon = 0.0 # aggregated longitude
|
|
lat_cell_size = 10.0 # Size in degrees latitude of the current cell. Starts at 10 and gets smaller as the calculation progresses
|
|
lon_cell_size = 20.0 # Size in degrees longitude of the current cell. Starts at 20 and gets smaller as the calculation progresses
|
|
|
|
# Iterate through blocks (two-character sections)
|
|
block = 0
|
|
while block * 2 < length:
|
|
if block % 2 == 0:
|
|
# Letters in this block
|
|
lon_cell_no = ord(grid[block * 2]) - ord('A')
|
|
lat_cell_no = ord(grid[block * 2 + 1]) - ord('A')
|
|
# Bail if the values aren't in range. Allowed values are A-R (0-17) for the first letter block, or
|
|
# A-X (0-23) thereafter.
|
|
max_cell_no = 17 if block == 0 else 23
|
|
if lat_cell_no < 0 or lat_cell_no > max_cell_no or lon_cell_no < 0 or lon_cell_no > max_cell_no:
|
|
return None, None, None, None
|
|
else:
|
|
# Numbers in this block
|
|
try:
|
|
lon_cell_no = int(grid[block * 2])
|
|
lat_cell_no = int(grid[block * 2 + 1])
|
|
except ValueError:
|
|
return None, None, None, None
|
|
# Bail if the values aren't in range 0-9
|
|
if lat_cell_no < 0 or lat_cell_no > 9 or lon_cell_no < 0 or lon_cell_no > 9:
|
|
return None, None, None, None
|
|
|
|
# Aggregate the angles
|
|
lat += lat_cell_no * lat_cell_size
|
|
lon += lon_cell_no * lon_cell_size
|
|
|
|
# Reduce the cell size for the next block, unless we are on the last cell.
|
|
if block * 2 < length - 2:
|
|
# Still have more work to do, so reduce the cell size
|
|
if block % 2 == 0:
|
|
# Just dealt with letters, next block will be numbers so cells will be 1/10 the current size
|
|
lat_cell_size = lat_cell_size / 10.0
|
|
lon_cell_size = lon_cell_size / 10.0
|
|
else:
|
|
# Just dealt with numbers, next block will be letters so cells will be 1/24 the current size
|
|
lat_cell_size = lat_cell_size / 24.0
|
|
lon_cell_size = lon_cell_size / 24.0
|
|
|
|
block += 1
|
|
|
|
# Offset back to (-180, -90) where the grid starts
|
|
lon -= 180.0
|
|
lat -= 90.0
|
|
|
|
# Return None values on maths errors
|
|
if any(x != x for x in [lat, lon, lat_cell_size, lon_cell_size]): # NaN check
|
|
return None, None, None, None
|
|
|
|
return lat, lon, lat_cell_size, lon_cell_size
|
|
|
|
|
|
def wab_wai_square_to_lat_lon(ref):
|
|
"""Convert a Worked All Britain or Worked All Ireland reference to a lat/lon point."""
|
|
|
|
# First check we have a valid grid square, and based on what it looks like, use either the Ordnance Survey, Irish,
|
|
# or UTM grid systems to perform the conversion.
|
|
if re.match(r"^[HNOST][ABCDEFGHJKLMNOPQRSTUVWXYZ][0-9]{2}$", ref):
|
|
return os_grid_square_to_lat_lon(ref)
|
|
elif re.match(r"^[ABCDEFGHJKLMNOPQRSTUVWXYZ][0-9]{2}$", ref):
|
|
return irish_grid_square_to_lat_lon(ref)
|
|
elif re.match(r"^W[AV][0-9]{2}$", ref):
|
|
return utm_grid_square_to_lat_lon(ref)
|
|
else:
|
|
logging.warning("Invalid WAB/WAI square: " + ref)
|
|
return None
|
|
|
|
|
|
def os_grid_square_to_lat_lon(ref):
|
|
"""Get a lat/lon point for the centre of an Ordnance Survey grid square"""
|
|
|
|
# Convert the letters into multipliers for the 500km squares and 100km squares
|
|
offset_500km_multiplier = ord(ref[0]) - 65
|
|
offset_100km_multiplier = ord(ref[1]) - 65
|
|
|
|
# The letter "I" is not used in the grid, so any offset of 8 or more needs to be reduced by 1.
|
|
if offset_500km_multiplier >= 8:
|
|
offset_500km_multiplier = offset_500km_multiplier - 1
|
|
if offset_100km_multiplier >= 8:
|
|
offset_100km_multiplier = offset_100km_multiplier - 1
|
|
|
|
# Convert the offsets into increments of 100km from the false origin (grid square SV):
|
|
easting_100km = ((offset_500km_multiplier - 2) % 5) * 5 + (offset_100km_multiplier % 5)
|
|
northing_100km = (19 - floor(offset_500km_multiplier / 5) * 5) - floor(offset_100km_multiplier / 5)
|
|
|
|
# Take the numeric parts of the grid square and multiply by 10000 to get metres, then combine with the 100km
|
|
# box offsets
|
|
easting = int(ref[2]) * 10000 + easting_100km * 100000
|
|
northing = int(ref[3]) * 10000 + northing_100km * 100000
|
|
|
|
# Add 5000m to each value to get the middle of the box rather than the south-west corner
|
|
easting = easting + 5000
|
|
northing = northing + 5000
|
|
|
|
# Reproject to WGS84 lat/lon
|
|
lat, lon = TRANSFORMER_OS_GRID_TO_WGS84.transform(easting, northing)
|
|
return lat, lon
|
|
|
|
|
|
def irish_grid_square_to_lat_lon(ref):
|
|
"""Get a lat/lon point for the centre of an Irish Grid square."""
|
|
|
|
# Convert the letters into multipliers for the 100km squares
|
|
offset_100km_multiplier = ord(ref[0]) - 65
|
|
|
|
# The letter "I" is not used in the grid, so any offset of 8 or more needs to be reduced by 1.
|
|
if offset_100km_multiplier >= 8:
|
|
offset_100km_multiplier = offset_100km_multiplier - 1
|
|
|
|
# Convert the offsets into increments of 100km from the false origin:
|
|
easting_100km = offset_100km_multiplier % 5
|
|
northing_100km = 4 - floor(offset_100km_multiplier / 5)
|
|
|
|
# Take the numeric parts of the grid square and multiply by 10000 to get metres, then combine with the 100km
|
|
# box offsets
|
|
easting = int(ref[1]) * 10000 + easting_100km * 100000
|
|
northing = int(ref[2]) * 10000 + northing_100km * 100000
|
|
|
|
# Add 5000m to each value to get the middle of the box rather than the south-west corner
|
|
easting = easting + 5000
|
|
northing = northing + 5000
|
|
|
|
# Reproject to WGS84 lat/lon
|
|
lat, lon = TRANSFORMER_IRISH_GRID_TO_WGS84.transform(easting, northing)
|
|
return lat, lon
|
|
|
|
|
|
def utm_grid_square_to_lat_lon(ref):
|
|
"""Get a lat/lon point for the centre of a UTM grid square (supports only squares WA & WV for the Channel Islands, nothing else implemented)"""
|
|
|
|
# Take the numeric parts of the grid square and multiply by 10000 to get metres from the corner of the letter-based grid square
|
|
easting = int(ref[2]) * 10000
|
|
northing = int(ref[3]) * 10000
|
|
|
|
# Apply the appropriate offset based on whether the square is WA or WV
|
|
easting = easting + 500000
|
|
if ref[1] == "A":
|
|
northing = northing + 5500000
|
|
else:
|
|
northing = northing + 5400000
|
|
|
|
# Add 5000m to each value to get the middle of the box rather than the south-west corner
|
|
easting = easting + 5000
|
|
northing = northing + 5000
|
|
|
|
# Reproject to WGS84 lat/lon
|
|
lat, lon = TRANSFORMER_CI_UTM_GRID_TO_WGS84.transform(easting, northing)
|
|
return lat, lon
|