CRS mismatch - Geopandas.clip

Question

I am working on a script that will iterate through a grid and clip out polygons that lay within their respective grid cell. I am running into a CRS mismatch error exclaiming that my left and right CRS's do not match. I tired setting them with to_crs(), but it changes my data type to None. My second issue that I run into is a No geometry data set yet (expected in column 'geometry'.) Well that should be because I no longer have a dataframe (class Nonetype).

Any advice?

import geopandas
import shapefile
import rasterio
import glob
import pandas

save_name = 0

little_polygon_path = 'D:/data/shape_merge/'
little_polygon_extent = geopandas.read_file(little_polygon_path + 'little_polygon_merge.shp')
little_polygon_extent.to_crs(epsg=32610, inplace=True)
little_polygon = geopandas.GeoDataFrame(little_polygon)

big_shape_extent_path = 'D:/data/extents/'
big_shape = geopandas.read_file('D:/data/extents/big_polygon.shp')
big_shape.to_crs(epsg=32610, inplace=True)
big_shape = geopandas.GeoDataFrame(big_shape)

for each_shape in little_polygon_extent.iterrows():
    df = geopandas.GeoDataFrame(each_shape[1])

    clip_shape = geopandas.clip(df, big_shape)
    xmin = clip_shape['geometry'].bounds['minx']
    xmax = clip_shape['geometry'].bounds['maxx']
    ymin = clip_shape['geometry'].bounds['miny']
    ymax = clip_shape['geometry'].bounds['maxy']

    p1 = ((xmin.iat[0]), (ymin.iat[0]))
    p2 = ((xmax.iat[0]), (ymin.iat[0]))
    p3 = ((xmax.iat[0]), (ymax.iat[0]))
    p4 = ((xmin.iat[0]), (ymax.iat[0]))

    shp_path = 'D:/data/clipped/'

    w = shapefile.Writer(shp_path + str(save_name) + '.shp')
    w.field("name", "C")  # pyshp needs at least one field
    w.poly([[p4, p3, p2, p1]])  # generate bbox polygon
    w.record('box')
    w.close()

    # generate .PRJ file
    tif_file = 'D:/farmer_32610.tif'
    image = rasterio.open(tif_file)
    crs_wkt = image.crs.to_wkt()
    prj = open(shp_path + str(save_name) + '.prj', "w")
    prj.write(crs_wkt)
    prj.close()

    save_name = save_name + 1



root = 'D:/data/clipped/'
path = 'D:/data/clipped_merge/'
shapefiles = glob.glob(root + '/' + '/*.shp')
gdf = pandas.concat([
geopandas.read_file(shp)
for shp in shapefiles
]).pipe(geopandas.GeoDataFrame)
# gdf.crs = {'init': 'epsg:32610'}
gdf.to_file(path + '/' + 'clipped_merge.shp')

Error Output before using to_crs(). After using to_crs() the CRS issue is cleared, but the "geometry" column issue occurs.

D:\PycharmProjects\one.py:23: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.

Left CRS: None
Right CRS: EPSG:32610

  clip_shape = geopandas.clip(df, big_shape)
Traceback (most recent call last):
  File "D:\PycharmProjects\one.py", line 23, in <module>
    clip_shape = geopandas.clip(df, big_shape)
  File "C:\Users\zlesl\anaconda3\envs\testing\lib\site-packages\geopandas\tools\clip.py", line 140, in clip
    box_gdf = gdf.total_bounds
  File "C:\Users\zlesl\anaconda3\envs\testing\lib\site-packages\pandas\core\generic.py", line 5462, in __getattr__
    return object.__getattribute__(self, name)
  File "C:\Users\zlesl\anaconda3\envs\testing\lib\site-packages\geopandas\base.py", line 2577, in total_bounds
    return GeometryArray(self.geometry.values).total_bounds
  File "C:\Users\zlesl\anaconda3\envs\testing\lib\site-packages\pandas\core\generic.py", line 5462, in __getattr__
    return object.__getattribute__(self, name)
  File "C:\Users\zlesl\anaconda3\envs\testing\lib\site-packages\geopandas\geodataframe.py", line 186, in _get_geometry
    raise AttributeError(
AttributeError: No geometry data set yet (expected in column 'geometry'.)

Edits. Needed to add these two lines because the to_crs() changes the data type to None. When I print out little_polygon and big_shape, I get a geopandas data frame that is empty.

little_polygon = geopandas.GeoDataFrame(little_polygon)
big_shape = geopandas.GeoDataFrame(big_shape)

Edits 2. So it seems that the CRS mismatch is actually a bug and can just be ignored. So I don't need to include the to_crs() anymore. With the updated code below I still however get the geometries error. little_polygon and big_shape are GeoDataFrames and calling head() to them outputs a column containing the geometry. I'll keep trouble shooting. I also changed iterrows to itertuples because there are some problems with the geometry data being currupted.

little_polygon_path = 'D:/data/shape_merge/'
little_polygon_extent = geopandas.read_file(little_polygon_path + 'little_polygon_merge.shp')


big_shape_extent_path = 'D:/data/extents/'
big_shape = geopandas.read_file('D:/data/extents/big_polygon.shp')


for each_shape in little_polygon_extent.itertuples():
    df = geopandas.GeoDataFrame(each_shape)

Edits 3. Using itertuples() changes the data type to a pandas data frame. I was able to index out the geometry and tried creating my own GeoDataFrame though am running into a different error: tuple object has no attribute 'is_empty'.

for each_shape in little_polygon_extent.itertuples():
    new = Polygon(each_shape[2])
    df = geopandas.GeoDataFrame({'geometry': [new.bounds]})

    shape_clip = geopandas.clip(df, big_shape, keep_geom_type=True)

Answer 1

Okay, I have finally figured it out. Along with my edits, I needed to pass a try block becuase some little_polygon were completely outside of my big_shape polygon. This will throw and error due to the each_shape variable will contain no data. So when p1 trys to grab xmin and ymin, there is no data to grab. The try block allows me to see if there is any data there. If there isn't go to the next polygon. If there is, complete the rest of the code.

import geopandas
import shapefile
import rasterio
import glob
import pandas
from shapely.geometry import Polygon

save_name = 0

little_polygon_path = 'D:/data/shape_merge/'
little_polygon_extent = geopandas.read_file(little_polygon_path + 'little_polygon_merge.shp')


big_shape_extent_path = 'D:/data/extents/'
big_shape = geopandas.read_file('D:/data/extents/big_polygon.shp')

for each_shape in little_polygon_extent.itertuples():
    new = Polygon(each_shape[2])
    df = geopandas.GeoDataFrame({'geometry': [new.bounds]})

    shape_clip = geopandas.clip(df, big_shape, keep_geom_type=True))

    xmin = bbox_clip['geometry'].bounds['minx']
    xmax = bbox_clip['geometry'].bounds['maxx']
    ymin = bbox_clip['geometry'].bounds['miny']
    ymax = bbox_clip['geometry'].bounds['maxy']

    try:
        p1 = ((xmin.iat[0]), (ymin.iat[0]))
    except:
        pass
    else:
        p1 = ((xmin.iat[0]), (ymin.iat[0]))
        p2 = ((xmax.iat[0]), (ymin.iat[0]))
        p3 = ((xmax.iat[0]), (ymax.iat[0]))
        p4 = ((xmin.iat[0]), (ymax.iat[0]))

        shp_path = 'D:/data/clipped/'

        w = shapefile.Writer(shp_path + str(save_name) + '.shp')
        w.field("name", "C")  # pyshp needs at least one field
        w.poly([[p4, p3, p2, p1]])  # generate bbox polygon
        w.record('box')
        w.close()

        # generate .PRJ file
        tif_file = 'D:/farmer_32610.tif'
        image = rasterio.open(tif_file)
        crs_wkt = image.crs.to_wkt()
        prj = open(shp_path + str(save_name) + '.prj', "w")
        prj.write(crs_wkt)
        prj.close()

        save_name = save_name + 1



root = 'D:/data/clipped/'
path = 'D:/data/clipped_merge/'
shapefiles = glob.glob(root + '/' + '/*.shp')
gdf = pandas.concat([
geopandas.read_file(shp)
for shp in shapefiles
]).pipe(geopandas.GeoDataFrame)
# gdf.crs = {'init': 'epsg:32610'}
gdf.to_file(path + '/' + 'clipped_merge.shp')