Test validity of CRS for a single coordinate in python

Question

Does anyone know if there is a simple way to check if a CRS (given by EPSG code) is valid for a single coordinate (lon, lat)?

Answer 1

Using pyproj and shapely

from shapely.geometry import box, Point
from pyproj import CRS

def within_crs_bounds(epsg_code, lon, lat):
    """ Returns `True` if point (`lon`, `lat`) is within the 
        bounds of the given projection (`epsg_code`).
        Otherwise returns `False`. """

    crs = CRS.from_user_input(epsg_code)

    bounding_box = box(*crs.area_of_use.bounds)

    return Point(lon, lat).within(bounding_box)

---

Using PyQGIS

from qgis.core import QgsCoordinateReferenceSystem, QgsPointXY

def within_crs_bounds(epsg_code, lon, lat):
    """ Returns `True` if point (`lon`, `lat`) is within the 
        bounds of the given projection (`epsg_code`).
        Otherwise returns `False`. """

    spatial_reference = QgsCoordinateReferenceSystem.fromEpsgId(epsg_code)
    
    bounding_box = spatial_reference.bounds()
    
    return bounding_box.contains(QgsPointXY(lon, lat))

---

Using GDAL/OGR

from osgeo import osr, ogr

def within_crs_bounds(epsg_code, lon, lat):
    """ Returns `True` if point (`lon`, `lat`) is within the 
        bounds of the given projection (`epsg_code`).
        Otherwise returns `False`. """

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg_code)

    area = srs.GetAreaOfUse()
    
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(lon, lat)

    min_x = area.west_lon_degree
    min_y = area.south_lat_degree
    max_x = area.east_lon_degree
    max_y = area.north_lat_degree

    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(min_x, min_y)
    ring.AddPoint(max_x, min_y)
    ring.AddPoint(max_x, max_y)
    ring.AddPoint(min_x, max_y)
    ring.AddPoint(min_x, min_y)

    box = ogr.Geometry(ogr.wkbPolygon)
    box.AddGeometry(ring)

    return point.Within(box)

Since creating geometries in GDAL leads to very verbose code, I would probably also check the spatial relation myself:

from osgeo import osr

def within_crs_bounds(epsg_code, lon, lat):
    """ Returns `True` if point (`lon`, `lat`) is within the 
        bounds of the given projection (`epsg_code`).
        Otherwise returns `False`. """

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(epsg_code)

    area = srs.GetAreaOfUse()

    # oddly formatted for readability
    return (  
            area. west_lon_degree <= lon <= area. east_lon_degree
        and area.south_lat_degree <= lat <= area.north_lat_degree
    )

---

Test Code

epsg_code = 2056  # Projection: CH1903+ / LV95

# Aarhus (Denmark)
assert not within_crs_bounds(epsg_code, lon=10.2107, lat=56.1572)

# Lausanne (Switzerland)
assert within_crs_bounds(epsg_code, lon=6.6335, lat=46.519833)

Answer 2

In the end I wrote the following function which seems to do the job:

def verify_crs_bounds(lon, lat, epsg_code):
    crs = pyproj.CRS.from_user_input(epsg_code)
    if (crs.area_of_use.west <= lon <= crs.area_of_use.east) and (
        crs.area_of_use.south <= lat <= crs.area_of_use.north
    ):
        result = True
    else:
        print(f"EPSG:{epsg_code} is invalid as the coordinate given is not within the permitted area of use.")
        result = False
    return result