0

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)

enter image description here

enter image description here


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_shape\103.shp: No such file or directory
ERROR 1: Terminating translation prematurely after failed
translation of layer compiled (use -skipfailures to skip errors)
5
  • 2
    You're resetting i to zero on every iteration, hence why it's not incrementing on every iteration – mikewatt Mar 4 at 23:29
  • Whoops, that would be a problem, let me update that. I am running into another error. Please see edits. – Zman3 Mar 5 at 0:19
  • It looks like I need to add -skipfailures to my code? :gis.stackexchange.com/questions/103161/… Though, I am not sure where to input this??? – Zman3 Mar 5 at 0:37
  • Okay I was able to get it to work after changing my cmd line to: cmd = 'ogr2ogr ' + outVectorPath + ' ' + inVectorPath + ' -clipsrc ' + extent + ' ' + ' -skipfailures ' – Zman3 Mar 5 at 0:55
  • Though, still when imported into QGIS, nothing is shown. Could that be due to not having anything in my attribute table? So I might need to add just a column heading with a name? – Zman3 Mar 5 at 0:56
0

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

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.