0

This is a follow on of this post.

I have the code below in order to reproject a shapefile according to the CSR/SRS of a raster. In this case the raster SRS is EPSG:4326 and the shapefile SRS is RGF93_Lambert_93 (EPSG:2154).

I read the raster with rasterio and the shapefile with geopandas (dataframe).

The reprojection fails with the Error:

AttributeError: 'NoneType' object has no attribute 'is_empty'



import os
import warnings
import geopandas as gpd
import rasterio as ras


# Given inputs
PRODUCT = 'GRD'
YEAR = '2020'
AOI = 'Lorraine'
Index = 'No'
    
hostname = os.environ['HOSTNAME']

# setting input folders
# For GRD products we provide the address to preprocessed data stack in the folder TS2
if hostname == 'snakhost':
    root_path = os.environ['HOME']+'/Documents/LIVE/Data/'
    
    #address to S1 time-series data for the interesting product type/area/year
    ts_path = root_path+PRODUCT+'/'+AOI+'/TS2/'+YEAR+'/'
    #address to the input shapefiles  
    shp_path = root_path +'SHP/'+AOI+'/'
    #address to the input .pkl files containing glolal raster stattistics
    pkl_path = root_path + 'PKL/'+AOI+'/'+YEAR+'/'
    #address to the output files containing zonal statistics
    zst_path = root_path +'ZST/'+AOI+'/'+YEAR+'/'

    
    
raster_dir = ts_path + [f for f in os.listdir(ts_path) if f.endswith('.data')][0]
raster_lst = os.listdir(raster_dir)    
input_raster = raster_dir+ '/'+raster_lst[40] #4 May Beta0 vv

shp_name = AOI +'_'+YEAR+'_'+ 'buffer_10_validation.shp' #Validation parcels    

# Open Raster Data
raster_src = ras.open(input_raster)

# Open Vector Data
vector_gdf = gpd.read_file(shp_path+shp_name)

# Make vector data to adopt the same projection as the rasters

#Check if the CRS of Raster and Vector data are the same 
if raster_src.crs.is_projected == False or vector_gdf.crs.is_projected == False:
    warnings.warn('At least one of the sources are not projected !\n')

if raster_src.crs.is_epsg_code == True:
    raster_epsg = raster_src.crs.to_epsg()
else:
    warnings.warn('The Raster is not in EPSG format !\n')
    
if (vector_gdf.crs.srs is not None 
    and vector_gdf.crs.srs.startswith('epsg') 
    or vector_gdf.crs.srs.startswith('EPSG')):
    vector_epsg = int(vector_gdf.crs.srs[5:])
    
elif (vector_gdf.crs.srs is not None 
    and vector_gdf.crs.srs.startswith('PROJCS') 
    and 'RGF93_Lambert_93' in vector_gdf.crs.srs):
    vector_epsg = vector_gdf.crs.to_epsg()  
    
#Addapting the CRS of the vector data with the CRS of the Raster


if vector_epsg != raster_epsg:
    warnings.warn('EPSG id of the vector data is not the same as of the raster.'\
                  'The CRS of the Vector data is going to be reprojected'\
                  ' according to the Raster \n')
    vector_gdf_reproj = vector_gdf.to_crs(raster_epsg)      # Reprojection
    #vector_gdf_reproj = vector_gdf.set_crs(epsg=raster_epsg) # Setting a Projection

Below are the complete error lines:

runfile('/home/sina/Documents/LIVE/Code/test_reprojection_GPD.py', wdir='/home/sina/Documents/LIVE/Code')
/home/sina/Documents/LIVE/Code/test_reprojection_GPD.py:55: UserWarning: At least one of the sources are not projected !

  warnings.warn('At least one of the sources are not projected !\n')
/home/sina/Documents/LIVE/Code/test_reprojection_GPD.py:77: UserWarning: EPSG id of the vector data is not the same as of the raster.The CRS of the Vector data is going to be reprojected according to the Raster 

  warnings.warn('EPSG id of the vector data is not the same as of the raster.'\
Traceback (most recent call last):

  File "/home/sina/Documents/LIVE/Code/test_reprojection_GPD.py", line 80, in <module>
    vector_gdf_reproj = vector_gdf.to_crs(raster_epsg)      # Reprojection

  File "/home/sina/anaconda3/envs/py37/lib/python3.7/site-packages/geopandas/geodataframe.py", line 816, in to_crs
    geom = df.geometry.to_crs(crs=crs, epsg=epsg)

  File "/home/sina/anaconda3/envs/py37/lib/python3.7/site-packages/geopandas/geoseries.py", line 543, in to_crs
    new_data = vectorized.transform(self.values.data, transformer.transform)

  File "/home/sina/anaconda3/envs/py37/lib/python3.7/site-packages/geopandas/_vectorized.py", line 889, in transform
    result[i] = transform(func, geom)

  File "/home/sina/anaconda3/envs/py37/lib/python3.7/site-packages/shapely/ops.py", line 227, in transform
    if geom.is_empty:

AttributeError: 'NoneType' object has no attribute 'is_empty'
9
  • 1
    what fails? is there an error message?
    – Ian Turton
    Jan 27 at 15:47
  • Sorry, I included the error in the text.
    – Sinooshka
    Jan 27 at 15:52
  • 1
    there is no call to is_empty in your code - do you get a line number or a stack trace of some type
    – Ian Turton
    Jan 27 at 16:03
  • I updated the text again including the complete error message.
    – Sinooshka
    Jan 27 at 16:11
  • one or more of your geometries is missing or broken
    – Ian Turton
    Jan 27 at 16:28

1 Answer 1

1

You wrote in a comment that your code behaves well with another shapefile and it means that there is nothing fundamentally broken in your code but the issue has some connection with that certain shapefile.

Generally I think that you are doing the conversion in too complicated way. If you do not need to know the EPSG codes for some reason you can just tell geopandas to reproject vector data into the same crs that raster has without knowing anything about the codes.

vector_gdf_reproj = vector_gdf.to_crs(raster_src.crs)

EDIT

I took some random files that I found from my hard disk and made the following test:

>>> import rasterio as ras
>>> vector_gdf = gpd.read_file('shops.shp')
>>> raster_src = ras.open('MRMS_gdal.tif')
>>> vector_gdf_reproj = vector_gdf.to_crs(raster_src.crs)
>>> vector_gdf.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

>>> vector_gdf_reproj.crs
<Geographic 2D CRS: GEOGCS["Coordinate System imported from GRIB file" ...>
Name: Coordinate System imported from GRIB file
Axis Info [ellipsoidal]:
- lat[north]: Latitude (degree)
- lon[east]: Longitude (degree)
Area of Use:
- undefined
Datum: unnamed
- Ellipsoid: Spheroid imported from GRIB file
- Prime Meridian: Greenwich

I did not get any errors and the vector_gdf_reproj.crs is different from vector_gdf.crs, and is the same than raster_src.crs. Everything worked as supposed. With these details I cannot say what fails with your data and with your conda installation.

EDIT

OpenJUMP shows that there are some empty geometries (geometry='MULTIPOLYGON EMPTY') in the source data

FID SURFACE_1   AUTEUR  DATE_       REMARQUE
542 0.21119848633   TD  2020-05-27  FAUCHEE 
586 0.14784282227   YC  2020-06-10  F   
717 0.25597026367   LP  2020-06-19  F   
744 0.26447241211   LH  2020-05-25  FAUCHEE 
745 0.25622050781   LH  2020-05-25  FAUCHEE 

Ogrinfo can find the empty geometries with command

ogrinfo -dialect SQLite -sql "select * from Lorraine_2020_buffer_10_validation where IsEmpty(geometry)=1 or IsEmpty(geometry)=-1" Lorraine_2020_buffer_10_validation.shp

Notice that it is IsEmpty(geometry)=-1 that finds these geometries because as documented in https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

The return type is Integer, with a return value of 1 for TRUE, 0 for FALSE, and –1 for UNKNOWN when called with a NULL argument.

There seems to a difference between a missing (null) geometry and an empty geometry.

13
  • Well theoretically you are right. In practice I followed your suggestion and the same error is produced !
    – Sinooshka
    Jan 31 at 13:07
  • I added a Python test that I did. Works for me.
    – user30184
    Jan 31 at 13:19
  • 1
    Some failure with finding empty geometries from shapefiles I guess. Not so big surprise, they are rare and easy to leave untested. But I believe that you can find those beasts by searching the geometries with null area. QGIS can also remove them with "Remove null geometries" tool.
    – user30184
    Jan 31 at 15:00
  • 1
    I want to confirm that the issue was with 5 empty (not_valid) geometries. by removing those bad apples, the (new) shapefile does not raise error while re-projecting it.
    – Sinooshka
    Jan 31 at 20:03
  • 1
    NULL geometries are valid generally and the shapefile specification allows them. The problem seems to be that the geopandas function is_empty fails when there is no geometry at all (null, NoneType). Maybe it should handle this case somehow better. AttributeError: 'NoneType' object has no attribute 'is_empty'
    – user30184
    Jan 31 at 21:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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