3

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)?

3
  • 1
    What do you mean by "check validity"?
    – Erik
    Aug 24, 2022 at 14:41
  • Sorry wasn't very clear of me. I mean check if the coordinate lies inside or outside the valid lon/lat ranges for the CRS. Aug 24, 2022 at 15:14
  • Look up the geographic bounding box and do a point in polygon test.
    – Ian Turton
    Aug 24, 2022 at 16:10

2 Answers 2

4

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)
1
  • Amazing! Thank you for such a thorough answer Thomas, I've learnt a lot from it :) Aug 25, 2022 at 12:48
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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.