4

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:

1

1 Answer 1

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.

2
  • This is awesome, thanks for the tip!!! For some reason, though, it doesn't work properly with my second example (EPSG:3081). I've added some notes to the bottom of my answer to explain the issues. Any idea of how to get around those limitations?
    – Felipe D.
    Commented Jul 8, 2021 at 17:20
  • I tried to find an answer to EPSG:3081. I could not find a method where those coordinates are obtained automatically. Good that you opened an issue in Github! I suggest you enter the bounds manually to solve your problem until a solution is implemented in Pyproj library. Commented Jul 9, 2021 at 5:26

Your Answer

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

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