Gdalwarp -crop_to_cutline results in shifted raster

[mjigmond]

When running gdalwarp using -cutline and -crop_to_cutline options, the resulting raster is slightly shifted. Both cutline shapefile and raster are using the same projection. On the sample raster dataset (pixel=28.5m) the shift is about 5m to the east and 3.5m to the north. Pixel values remain the same (so resampling doesn't take place, as expected).

I should add that not adding -crop_to_cutline results in no shift. Commands I ran:

(1) gdalwarp -cutline cluj_co.shp -crop_to_cutline sample.tif cropped.tif

result is shifted

(2) gdalwarp -cutline cluj_co.shp sample.tif out.tif

result is not shifted

Sample dataset was too large for attachment. Temporary link provided below: ​https://www.twdb.state.tx.us/filetxfr/emaillinkdownload.aspx?id=&FileTransferID=774740096

-marius

[Even Rouault]

I don't think this is really a bug, but rather an enhancement.

First, -crop_to_cutline is just an automated way of computing the xmin ymin xmax ymax values that you would pass with the -te option of gdalwarp. So it will behave exactly like if you used -te directly.

When specifying -te, the warping algorithm must recompute the appropriate pixel size, then the image height/width and the actual raster extents. So the origin of the new raster doesn't necessarily match with the origin of a pixel of the source raster, hence this apparent "shift". But notice it remains under the pixel size of the source raster, which is the precision of the data, so I think this is acceptable from a theoretical point of view (yes in practice if you overlay the 2 rasters, you'll notice a visual shift...)

Perhaps in that particular case (no reprojection, just one input raster), a specific code could be made to alter the "raw" bounding box given by the cropping polygons so that it matches the coordinates of the border of source pixels. In all other cases, there's no way of guaranteing a 0 shift.

[mjigmond]

Replying to rouault:

I don't think this is really a bug, but rather an enhancement.

First, -crop_to_cutline is just an automated way of computing the xmin ymin xmax ymax values that you would pass with the -te option of gdalwarp. So it will behave exactly like if you used -te directly.

When specifying -te, the warping algorithm must recompute the appropriate pixel size, then the image height/width and the actual raster extents. So the origin of the new raster doesn't necessarily match with the origin of a pixel of the source raster, hence this apparent "shift". But notice it remains under the pixel size of the source raster, which is the precision of the data, so I think this is acceptable from a theoretical point of view (yes in practice if you overlay the 2 rasters, you'll notice a visual shift...)

Perhaps in that particular case (no reprojection, just one input raster), a specific code could be made to alter the "raw" bounding box given by the cropping polygons so that it matches the coordinates of the border of source pixels. In all other cases, there's no way of guaranteing a 0 shift.

I agree that if actual warping was involved you can no longer guarantee the origin. But when input datasets use identical projections (and a different output projection is not being requested) I believe -crop_to_cutline should behave more like -projwin from gdal_translate which forces the computation of -srcwin by maintaining the relative origin.

A secondary issue which I discovered after the submission (possibly another enhancement request) is that the output raster extent is not greater than or equal to the shapefile extent. This results in vector areas not being covered by the output raster.

I disagree that this should be marked as only an enhancement. Because rasters cannot be re-projected, any operations applied to them should be minimally disruptive. Subsequent -crop_to_cutline operations on an already shifted raster only compounds the problem. Or even warping to a new projection will amplify the otherwise normal warping errors.

-marius

[Even Rouault]

To avoid "polluting" this ticket with other things, could you open a ticket for the secondary issue you mention ? I agree that the output extent should be the intersection of the extent of the input raster files and the shape extent.

As far as the main issue of the ticket is concerned, I'd note that currently if you do "gdalwarp in.tif out.tif" and don't specify -tr, you can have resampling to occur if the pixel width is different from the pixel height (in your case, they are equal), so it would also imply altering this behaviour, but perhaps some people actually prefer it...

[mjigmond]

Replying to rouault:

To avoid "polluting" this ticket with other things, could you open a ticket for the secondary issue you mention ? I agree that the output extent should be the intersection of the extent of the input raster files and the shape extent.

Will do (pending the conclusion of this discussion).

As far as the main issue of the ticket is concerned, I'd note that currently if you do "gdalwarp in.tif out.tif" and don't specify -tr, you can have resampling to occur if the pixel width is different from the pixel height (in your case, they are equal), so it would also imply altering this behaviour, but perhaps some people actually prefer it...

Is there a reason square pixels are treated differently than non-square pixels? Should my output had been resampled, then? I'm not looking to open a can of worms :), just looking for consistent behavior.

Our discussion got me thinking. Ideally overlay operations should occur only on identically projected datasets to minimize the errors (obviously the user's responsibility). Particularly, -crop_to_cutline (extract by mask) is usually performed to downsize a large raster to a specific area of interest and then, later, extract the raster statistics for that area (aside from the visual aspect).

Wouldn't -crop_to_cutline be better suited for gdal_translate? After all, when executing -crop_to_cutline I would think most users are not trying to resample the raster, but simply extract the raster pixels overlaying the area of interest.

-marius

[Even Rouault]

Is there a reason square pixels are treated differently than non-square pixels? Should my output had been resampled, then? I'm not looking to open a can of worms :), just looking for consistent behavior.

square and non square pixels are not treated differently in fact. This is just the consequence of how the algorithm that computes the spatial extent and the resolution of the output raster after reprojection work (the no reprojection case is not a particular case, but just a transform that does nothing). It is written such that the computed output pixel width = pixel height, and that the area of the output pixel is close enough to the area of the input pixel. So, when your input raster has square pixel and you don't do any reprojection, then the output raster will have the same pixel size. But if the input is a non square pixel, you should get a pixel size = sqrt(src_pixel_width * src_pixel_width + src_pixel_height * src_pixel_height)

Our discussion got me thinking. Ideally overlay operations should occur only on identically projected datasets to minimize the errors (obviously the user's responsibility). Particularly, -crop_to_cutline (extract by mask) is usually performed to downsize a large raster to a specific area of interest and then, later, extract the raster statistics for that area (aside from the visual aspect).

Wouldn't -crop_to_cutline be better suited for gdal_translate? After all, when executing -crop_to_cutline I would think most users are not trying to resample the raster, but simply extract the raster pixels overlaying the area of interest.

Perhaps, but I reiterate that -crop_to_cutline is no more than a stupid (~ 20 lines of code) and convenient way of automatically filling the target extent of the raster, which you could do manually with the -te option. The hard work of masking raster by feature geometries is done by the -cutline option, which is an optionnal step of the warping processing. And there's currently no infrastructure to do cutline in gdal_translate. Would require work in the VRT driver for this.

[mjigmond]

Replying to rouault:

Is there a reason square pixels are treated differently than non-square pixels? Should my output had been resampled, then? I'm not looking to open a can of worms :), just looking for consistent behavior.

square and non square pixels are not treated differently in fact. This is just the consequence of how the algorithm that computes the spatial extent and the resolution of the output raster after reprojection work (the no reprojection case is not a particular case, but just a transform that does nothing). It is written such that the computed output pixel width = pixel height, and that the area of the output pixel is close enough to the area of the input pixel. So, when your input raster has square pixel and you don't do any reprojection, then the output raster will have the same pixel size. But if the input is a non square pixel, you should get a pixel size = sqrt(src_pixel_width * src_pixel_width + src_pixel_height * src_pixel_height)

If I understand correctly, the pixel essentially becomes square and resampling takes place.

Our discussion got me thinking. Ideally overlay operations should occur only on identically projected datasets to minimize the errors (obviously the user's responsibility). Particularly, -crop_to_cutline (extract by mask) is usually performed to downsize a large raster to a specific area of interest and then, later, extract the raster statistics for that area (aside from the visual aspect).

Wouldn't -crop_to_cutline be better suited for gdal_translate? After all, when executing -crop_to_cutline I would think most users are not trying to resample the raster, but simply extract the raster pixels overlaying the area of interest.

Perhaps, but I reiterate that -crop_to_cutline is no more than a stupid (~ 20 lines of code) and convenient way of automatically filling the target extent of the raster, which you could do manually with the -te option. The hard work of masking raster by feature geometries is done by the -cutline option, which is an optionnal step of the warping processing. And there's currently no infrastructure to do cutline in gdal_translate. Would require work in the VRT driver for this.

In light of this, I'm inclined to request closing this ticket and perhaps suggest a feature request. I feel that -cutline and -crop_to_cutline are more appropriate in the gdal_translate world where unexpected warping/resampling and shifting are not likely to occur. What are your thoughts on this?

-marius

[antoniorocha]

Was there any feature request or even any update?

[mjigmond]

Antonio, I had not had a chance to request the enhancement as my life took a busy turn. I do plan to make the feature request where "-cutline" and "-crop_to_cutline" be added (or moved) to gdal_translate. There was a side issue that came out of the testing I had done and I will need to test trunk to see if a fix slipped in. Based on my discussion with Even (and I agree with him) I don't think the current behavior of gdalwarp can be changed.

[mjigmond]

I propose we close this ticket and reference the request for enhancement, #4875.

[Jukka Rahkonen]

Suggestion made by the creator of the ticket accepted.