Clipping Shapefile to Raster Extent using GDAL with Python

Question

I am trying to clip a shapefile (containing multiple polygons) to a raster's extent. A similar question was asked as Clipping polygon precisely to raster (extent) using GDAL?, with an answer that I have not been able to effectively use. In my case, I have a large raster file that contains one shapefile which holds over 6000 polygons. I split my raster up into equal segments (600x600). I want to also split the polygons up into their respective new .tif file.
The code below does run without any errors though it saves over each shapefile that is saved. So my variable i is not increasing. When I import the shapefile into QGIS it loads correctly, but has no data. So it is a blank shapefile. If I zoom to layer, it does not zoom. So I am not sure what is going on.
Any suggestions?
import gdal
import subprocess
import os

# define file paths
Ras_path = 'D:/images/'
i = 0
for Raster in os.listdir('D:/images/'):
    src = gdal.Open(Ras_path + Raster)
    ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform()
    sizeX = src.RasterXSize * xres
    sizeY = src.RasterYSize * yres
    lrx = ulx + sizeX
    lry = uly + sizeY
    src = None

extent = '{0} {1} {2} {3}'.format(ulx, lry, lrx, uly);

#i = 0
    inVectorPath = 'D:/data/shape_merge/compiled.shp'
    outVectorPath = 'D:/data/split_shape/' + str(i) + '.shp'
    i = i + 1

# make clip command with ogr2ogr - default to shapefile format
    cmd = 'ogr2ogr ' + outVectorPath + ' ' + inVectorPath + ' -clipsrc ' + extent

# call the command
    subprocess.call(cmd, shell=True)

As @mikewatt commented, I did have my i = 0 inside my loop so it wasn't incrementing correctly. After the fix I am running into an error:
ERROR 1: Failed to create file D:/data/split_shape103.shp: No such file or directory
ERROR 1: Terminating translation prematurely after failed
translation of layer compiled (use -skipfailures to skip errors)

Zman3 · Accepted Answer

I ended up taking a pretty round about way, but I did manage to achieve my goal. Though, I had to use QGIS and I would like to just use purely python for efficiency and automation.
From the .tif I used, I extracted the extent and used the xmin, ymin, xmax, ymax values to create basically a grid.
import shapefile
import os
import shutil
import gdal
import glob
import geopandas
import pandas

save_name = 0

Ras_path = 'D:/images/'
shp_path = 'D:/data/extent/'
for Raster in os.listdir('D:/images/'):
    src = gdal.Open(Ras_path + Raster)
    xmin, xres, xskew, ymax, yskew, yres = src.GetGeoTransform()
    sizeX = src.RasterXSize * xres
    sizeY = src.RasterYSize * yres
    xmax = xmin + sizeX
    ymin = ymax + sizeY

extent = '{0} {1} {2} {3}'.format(xmin, ymin, xmax, ymax)

# get vertices of the box
    # geocoordinates from pixel coordinates
    p1 = (xmin, ymin)
    p2 = (xmax, ymin)
    p3 = (xmax, ymax)
    p4 = (xmin, ymax)

# save shapefile containing one bounding box shape
    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('bbox')
    w.close()

# generate .PRJ file
    crs_wkt = src.GetSpatialRef()
    prj = open(shp_path + str(save_name) + '.prj', "w")
    prj.write(crs_wkt.ExportToWkt())
    prj.close()

save_name = save_name + 1

root = 'D:/data/extent/'
path = 'D:/data/extent_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 + '/' + 'compiled.shp')

I then took the .shp identified as the red boxes in the question asked pictures and performed QGIS's clip method, iterating over each extent. Some of the extents did not contain any polygons in them, but would still save, so I created a script to remove those who did not hold any data. I found that the size of the data was 100 bytes:
import os

path = 'D:/data/clipped/'
for clipped_shape in os.listdir(path):
    b = os.path.getsize(path + clipped_shape)
    if b <= 100:
        os.remove(path + clipped_shape)
    else:
        pass

Clipping Shapefile to Raster Extent using GDAL with Python

One Answer

Add your own answers!

Ask a Question