6

I'm currently trying to cutline crop a tif with a shapefile and no matter what I do I have pixel either overhanging the cutline or being cropped too early.

So here is an example picture of what is happening when I use gdalwarp: enter image description here

When I use arcpy to crop I get this result:

enter image description here

The code I am currently using with gdal warp is:

 gdalwarp.exe -ot "float32" -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 
 +no_defs " -t_srs "+proj=utm
 +zone=48 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "-of "GTiFF" -cutline 
 mask.shp -crop_to_cutline infile.asc outfile.tif -wm 5000 -tap -tr 100 100 -overwrite

I am quite keen to use gdal if possible to do this cropping as it is considerably faster


I have been asked to provide gdalinfo info, see below. This is from the gdal code: enter image description here

This is from the arcpy code: enter image description here

SOLUTION was found in comments.

Exactly what I wanted:

enter image description here

My new problem to ponder next is: [enter image description here]

11
  • 1
    What is the output from a gdalinfo on the arcpy result? Commented Feb 12, 2018 at 14:02
  • 1
    And is the resolution of the input data also 100m? Commented Feb 12, 2018 at 14:42
  • 1
    The resolution of your arcpy image seems to be much higher along the cutline. Or is the cutline covering part of pixels? Are the visible blocks in the raster 100 by 100 meters in size?
    – user30184
    Commented Feb 12, 2018 at 14:53
  • 2
    You could try setting -wo CUTLINE_ALL_TOUCHED=TRUE, so that all pixels overlapping the cutline polygon will be selected, not just those whose center point falls within the polygon. see gdal.org/…
    – dr_times
    Commented Feb 12, 2018 at 15:22
  • 1
    Please consider writing the solution as an answer, and accepting it. This will mark the question as answered for future reference. Commented Feb 12, 2018 at 16:46

2 Answers 2

8

The solutiuon was to add the option: -wo CUTLINE_ALL_TOUCHED=TRUE

0

In trying to spatially upscale a raster and ensure the upscaled perimeter completely overlapped the source raster I developed the following solution that could also be adapted and used here.

# Define input and output files
INPUT_FILE="source_raster.tif"
OUTPUT_FILE="final_raster.tif"  # Final output after masking
WARPED_FILE="temporary_raster.tif"  # Temporary output before masking
INVERTED_MASK_FILE="inverted_mask.tif"
MASK_FILE="mask.tif"
CLIP_MASK="clip_mask.shp"

# Step 1: Create an inverted mask to define the clipping area
# - Invert the mask: Set data pixels to 0 (keep), NODATA to 1 (clip)

gdal_calc.py --overwrite -A $INPUT_FILE --outfile=$INVERTED_MASK_FILE
--calc="1*(A>0)" --NoDataValue=0

#  - A: refers to the input raster in the calculation
#  - --calc: expression to create the mask (data > 0 = keep, set to 1)
#  - --NoDataValue: set NODATA pixels in the output to 0

# Step 2: Generate a shapefile from the inverted mask for clipping 
gdal_polygonize.py $INVERTED_MASK_FILE -f "ESRI Shapefile" $CLIP_MASK

#  - Convert inverted mask to a shapefile defining the clipping boundary

# Step 3: Perform the upscaling with clipping

gdalwarp -cutline "$CLIP_MASK" -crop_to_cutline -wo CUTLINE_ALL_TOUCHED=TRUE -overwrite -tr 1 1 -r average -of GTiff "$INPUT_FILE" "$WARPED_FILE"

#  - -cutline: specify the shapefile for clipping
#  - -crop_to_cutline: clip the output to the shapefile boundary
#  - -wo CUTLINE_ALL_TOUCHED: ensure all features in the shapefile are considered
#  - -tr 1 1: set output pixel size to 1 degree
#  - -r average: resampling method (averaging)

# Step 4: Create a mask of valid data in the upscaled raster

gdal_calc.py --overwrite -A "$WARPED_FILE" --outfile="$MASK_FILE" --calc="A*(A!=0)" --NoDataValue=0

#  - Create a mask where valid data pixels are set to 1 (keep)

# Step 5: Apply the mask to set zero values (nodata in upscaled) as NODATA in the final output

gdal_calc.py --overwrite -A "$WARPED_FILE" -B "$MASK_FILE" --outfile="$OUTPUT_FILE" --calc="A*(B>0)" --NoDataValue=0

#  - A: upscaled raster with potential NODATA values
#  - B: mask with valid data identified
#  - The calculation sets pixels in A to NODATA if the corresponding pixel in B is 0 (area to clip)

# Clean up temporary file
rm $WARPED_FILE

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.