mirror of
https://git.ianrenton.com/ian/spothole.git
synced 2025-10-27 16:59:25 +00:00
104 lines
4.4 KiB
Python
104 lines
4.4 KiB
Python
import logging
|
|
import re
|
|
from math import floor
|
|
|
|
from pyproj import Transformer
|
|
|
|
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")
|
|
|
|
|
|
# Convert a Worked All Britain or Worked All Ireland reference to a lat/lon point.
|
|
def wab_wai_square_to_lat_lon(ref):
|
|
# 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.warn("Invalid WAB/WAI square: " + ref)
|
|
return None
|
|
|
|
|
|
# Get a lat/lon point for the centre of an Ordnance Survey grid square
|
|
def os_grid_square_to_lat_lon(ref):
|
|
# 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
|
|
|
|
|
|
# Get a lat/lon point for the centre of an Irish Grid square.
|
|
def irish_grid_square_to_lat_lon(ref):
|
|
# 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
|
|
|
|
|
|
# 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)
|
|
def utm_grid_square_to_lat_lon(ref):
|
|
# 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
|