1

First time posting so apologies if I miss important info. I keep getting this error but it only happens intermittently so I have had a hard time narrowing down the problem. I have an input of a geodataframe and I think the error is related to when it tries to buffer the geometry. Below is the code and the resulting error. I am reprojecting to the local UTM Zone using an epsg lookup method self.df.estimate_utm_crs(datum)name = 'NAD 83') before buffering. I have also used some geometry validation: removing empty geometries, check for invalid geometries, removing any disconnected parts.

Geometry: (-79.403895, 39.3028665, -79.3876355, 39.3174215)

Using Python 3.9.23

# Reproject back to WGS84
parcel_wgs84 = buffered_df.to_crs(epsg=4326)

# Get bbox as a tuple 
parcel_wkt = str(df['geometry'][0]) 
bbox = sh.wkt.loads(parcel_wkt)
bounds = bbox.bounds
    
# # Extract DEM 
# Using py3dep...
print(py3dep.check_3dep_availability(bounds, 4326))
dem = py3dep.get_map("DEM", bounds, 1, crs=4326)```

Error:

Cell In[8], line 1
----> 1 dem = py3dep.get_map("DEM", bounds, 10, crs=4326)

File ~\AlleteDevRepo\Anaconda\envs\contours\lib\site-packages\py3dep\py3dep.py:175, in get_map(layers, geometry, resolution, geo_crs, crs)
    173 req_crs = 5070
    174 geom_5070 = geoutils.geo2polygon(geometry, geo_crs, req_crs)
--> 175 bbox = geom_5070.buffer(20 * resolution).bounds
    176 wms = WMS(wms_url, layers=_layers, outformat="image/tiff", crs=req_crs, validation=False)
    177 r_dict = wms.getmap_bybox(bbox, resolution, box_crs=req_crs, max_px=MAX_PIXELS)

File ~\AlleteDevRepo\Anaconda\envs\contours\lib\site-packages\shapely\geometry\base.py:545, in BaseGeometry.buffer(self, distance, quad_segs, cap_style, join_style, mitre_limit, single_sided, **kwargs)
    542 elif not np.isfinite(distance).all():
    543     raise ValueError("buffer distance must be finite")
--> 545 return shapely.buffer(
    546     self,
    547     distance,
    548     quad_segs=quad_segs,
    549     cap_style=cap_style,
    550     join_style=join_style,
    551     mitre_limit=mitre_limit,
    552     single_sided=single_sided,
    553 )

File ~\AlleteDevRepo\Anaconda\envs\contours\lib\site-packages\shapely\decorators.py:77, in multithreading_enabled.<locals>.wrapped(*args, **kwargs)
     75     for arr in array_args:
     76         arr.flags.writeable = False
---> 77     return func(*args, **kwargs)
     78 finally:
     79     for arr, old_flag in zip(array_args, old_flags):

File ~\AlleteDevRepo\Anaconda\envs\contours\lib\site-packages\shapely\constructive.py:180, in buffer(geometry, distance, quad_segs, cap_style, join_style, mitre_limit, single_sided, **kwargs)
    178 if not np.isscalar(single_sided):
    179     raise TypeError("single_sided only accepts scalar values")
--> 180 return lib.buffer(
    181     geometry,
    182     distance,
    183     np.intc(quad_segs),
    184     np.intc(cap_style),
    185     np.intc(join_style),
    186     mitre_limit,
    187     np.bool_(single_sided),
    188     **kwargs
    189 )

GEOSException: Shell empty after removing invalid points```
Improve this question
5
  • could you share the geometry causing this error? Ideally you make the script included into a minimal reproducing script by removing irrelevant code and adding the geometry that causes the issue. E.g. include the geometry in the script as wkt, create a Geodataframe from it using something like df = gpd.GeodataFrame(geometry=gpd.GeoSeries.from_wkt([wkt]), crs="EPSG:4326") Commented Sep 29 at 15:28
  • I edited it, I think that should hopefully make it more clear. I believe the py3dep library requires a tuple as the argument for the bounds (I updated this as the Geometry listed towards the top of the question). Commented Sep 29 at 17:16
  • @Richelle Have you tried applying a zero-width buffer? And also check with is_valid and is_valid_reason. geopandas.org/en/stable/docs/reference/api/… Commented Sep 29 at 17:39
  • I have not tried the is_valid_reason so I will check that out Commented Sep 30 at 15:02
  • This Geometry: (-79.403895, 39.3028665, -79.3876355, 39.3174215) is not a valid WKT geometry. Can you do df["geometry"].to_wkt().iloc[0] and add the output to your question if is is the first [0] geometry that is the problem Commented Oct 28 at 6:43
Related questions
0
Error using geopandas overlay with union method
3
Why is env.render() not working on Jupyter?
1
Wrong points buffer using geopandas
Load 6 more related questions Show fewer related questions

0

Reset to default

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.