CQ/ITU zone lookups

core/geo_utils.py

@@ -2,12 +2,136 @@ import logging
 import re
 from math import floor
 
+import geopandas
 from pyproj import Transformer
+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")
 
+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"))
+
+
+# Finds out which CQ zone a lat/lon point is in.
+def lat_lon_to_cq_zone(lat, lon):
+    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"])
+    return None
+
+
+# Finds out which ITU zone a lat/lon point is in.
+def lat_lon_to_itu_zone(lat, lon):
+    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"])
+    return None
+
+
+# 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.
+def lat_lon_for_grid_centre(grid):
+    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
+
+
+# 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.
+def lat_lon_for_grid_sw_corner(grid):
+    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
+
+# 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.
+def lat_lon_for_grid_ne_corner(grid):
+    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
+
+
+# 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.
+def lat_lon_for_grid_sw_corner_plus_size(grid):
+    # 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
+
 
 # Convert a Worked All Britain or Worked All Ireland reference to a lat/lon point.
 def wab_wai_square_to_lat_lon(ref):
@@ -20,7 +144,7 @@ def wab_wai_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)
+        logging.warning("Invalid WAB/WAI square: " + ref)
         return None