0

I am having trouble clipping a raster with GDAL. I have the coordinates of the bounding box I would like to clip to. There is no error but also no output file is being saved. Could it have to do with my coordinate system? My coordinates are in a UTM zone (same as raster I am clipping)

import arcpy
from arcpy import env
from arcpy.sa import *
import os
import subprocess
import gdal

path = "R:\\Sentinel-1\\kevf\\"
image = "S1A_IW_GRDH_1SDH_20190303T100649_20190303T100714_026172_02EBE6_AD79_out"
segmentedRasterSave = "R:\\Sentinel-1\\kevf\\seg\\S1A_IW_GRDH_1SDH_20190303T100649_20190303T100714_026172_02EBE6_AD79_out_seg.tif"

def split(segmentedRasterSave, image):
    arcpy.CheckOutExtension("Spatial")

    elevRaster = arcpy.sa.Raster(segmentedRasterSave)
    Extent = elevRaster.extent
    #set values
    xL = float(Extent.XMin)
    xR = float(Extent.XMax)

    yB = float(Extent.YMin)
    yT = float(Extent.YMax)

    x_delta = xL - xR
    x1 = xL + float((x_delta/3))
    x2 = xR - float((x_delta/3))

    leftSave = path + "split\\" + "_left.tif"

    ds = gdal.Translate(leftSave, gdal.Open(segmentedRasterSave), projWin = [xL, yB, x1, yT])
    ds = None

split(segmentedRasterSave, image)

There is no error message, just no output file is being saved.

  • 1
    Use os.path.join instead of combining paths like this: path + "split\\" + "_left.tif" – BERA Mar 24 at 19:13
  • 2
    Why arent you using Clip? – BERA Mar 24 at 19:17
  • clip will write out a file that is not clipped to my dimensions. – Kevin Ferris Mar 24 at 19:36
  • clip won't even work on my raster in arcPRO. just returns full extent of input – Kevin Ferris Mar 24 at 19:55
1

You can use the Extract By Mask tool which:

Extracts the cells of a raster that correspond to the areas defined by a mask.

The arcpy function is:

ExtractByMask(in_raster, in_mask_data)

You can convert your bounding box to a shapefile (if you don't have one already) and pass that file as the in_mask_data parameter.

| improve this answer | |
  • that's what my plan is now. Any suggestions for easily converting a bounding box to a shapefile? – Kevin Ferris Mar 25 at 13:42
  • @Kevin Ferris At the end of the arcpy.Polygon documentation, there is a code sample of how to create a shapefile from a set of coordinates. You would have to change the points to match your bounding box and set the spatial reference to the UTM zone you are using. – Marcelo Villa-Piñeros Mar 25 at 15:16
  • I keep getting an empty shapefile for whatever reason – Kevin Ferris Mar 25 at 16:49
  • Post a new question with the code you are running and link it here so I can take a look. – Marcelo Villa-Piñeros Mar 25 at 20:25

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.