How to check which geometries in a GeoDataFrame are within the bounds/limits of the projection/coordinate system?

Question

Is there a pre-built function that I can use to identify which of the geometries in a GeoDataFrame fall within the bounds of a CRS?

I'm asking because, currently, GeoPandas allows me to create points that are wildly outside of a CRS' bounds.

For example, WGS84 (EPSG:4326) has the following bounds:

• x_min=-180, x_max=180, y_min=-90, y_max=90

However, I can create POINT (-500 900) and GeoPandas won't throw any kind of error.

Similarly, the Texas State Mapping System (EPSG:3081) has the following bounds:

• x_min=334235.2323, x_max=1651722.6917, y_min=428079.3639, y_max=1608994.9702

However, I can create POINT (-50000000 70000000) and GeoPandas won't throw any kind of error.

Here's a quick example (in code) of what I'm talking about.

import geopandas as gpd

test_1 = gpd.GeoDataFrame({'my_id':[1,2,3]},
                          geometry=gpd.points_from_xy([-180,180,-500], 
                                                      [-90,90,900]),
                          crs='epsg:4326')


test_2 = gpd.GeoDataFrame({'my_id':[11,12,13]},
                          geometry=gpd.points_from_xy([334235.2323, 
                                                       1651722.6917, 
                                                       -50000000], 
                                                      [428079.3639, 
                                                       1608994.9702, 
                                                       70000000]),
                          crs='epsg:3081')

Notice how in GeoPandas allowed the creation of both GeoDataFrames without any hiccups.

So, my question is: is there an easy way to identify geometries that fall outside of a CRS' "legal" bounds?

Edit

I've built on top of what @Babak Fi Foo suggested and created the following function:

import pyproj
def within_crs_bounds(gdf):
    ref_crs = pyproj.CRS.from_epsg(gdf.crs.to_epsg())
    
    x_min, y_min, x_max, y_max = ref_crs.area_of_use.bounds
    
    proj = pyproj.Transformer.from_crs(4326, ref_crs.to_epsg(), always_xy=ref_crs.is_projected)
    
    crs_x_min, crs_y_min = proj.transform(x_min, y_min)
    crs_x_max, crs_y_max = proj.transform(x_max, y_max)
    
    gdf_bounds = gdf.geometry.bounds
    result = ((gdf_bounds['minx']>=crs_x_min)&
              (gdf_bounds['miny']>=crs_y_min)&
              (gdf_bounds['maxx']<=crs_x_max)&
              (gdf_bounds['maxy']<=crs_y_max))
    return result

For the test_1 GeoDataFrame, this implementation works just fine. However, for the test_2 GeoDataFrame, it does not work. The crs_x_min, crs_y_min, crs_x_max and crs_y_max variables get values far away from the true bounds (as seen on the EPSG:3081 Spatial Reference Site). I expect them to hold (334235.2323, 428079.3639, 1651722.6917, 1608994.9702), but they actually hold (331226.35401596106, 429093.9412772637, 1583305.6741846378, 1608302.6853412518) instead.

Note

Related GitHub issues and pull requests:

• Issue: ENH: Method that checks which geometries in a GeoDataFrame are within the bounds/limits of the projection/coordinate system Issue #1915
• Issue: would it be possible to have the CRS bounds in its original unit ? Issue #809
• PR: ENH: Add Transformer.transform_bounds

Answer 1

First you need to get the bounds of the CRS. you can find it here.

from pyproj import CRS
import geopandas as gpd

crs= CRS.from_epsg(4326)
x_min, y_min, x_max, y_max = crs.area_of_use.bounds

Then you can get the bounds of each geometry.

gdf[['minx','miny','maxx','maxy']] = gdf.geometry.bounds

What remains is filtering them out:

    gdf = gdf[
             (gdf['minx']>x_min)&
             (gdf['miny']>y_min)&
             (gdf['maxx']<x_max)&
             (gdf['maxy']<y_max)&
             ]

This way all shapes that are within the bounds of CRS will remain in your geo data frame.