Reprojection from RGF93_Lambert_93 fails

Question

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'

Answer 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.