Reproject and crop are not commutative

[jmarintur]

In the first case I perform crop and then reproject:

1) raster.crop(roi.envelope).reproject(dst_crs=WEB_MERCATOR_CRS).save(output_raster)

In the second case I reproject and then I crop:

2) raster.reproject(dst_crs=WEB_MERCATOR_CRS).crop(roi.envelope).save(output_raster)

The resulting images when I visualize them using Qgis 2.14.11-Essen look like this:

1)

And

2)

I am using a Sentinel-2 image Band 02 and tile T18HXB and the following roi:

coord = [-8083749.056641963, -4737939.582734095, -8079749.056641963, -4733939.582734095]
roi = tl.GeoVector.from_bounds(xmin=coordindates[0], ymin=coordindates[1], xmax=coordindates[2], ymax=coordindates[3], crs=WEB_MERCATOR_CRS)

Thanks for your support.

[astrojuanlu]

Thanks a lot for the good bug report @jmarintur! Could you please give the output of crs.raster for completeness?

[jmarintur]

This is it: CRS({'init': 'epsg:32718'})

Thanks @Juanlu001 !

[drnextgis]

It is expected that shapes in these cases are different (at least from code perspective). Let's see what is happening.

First case (crop first):

1 roi.envelope.get_bounds(self.crs) is called inside crop here 2 roi.envelope in original CRS epsg:3857 is not the same as bounds of the roi.envelope in epsg:32718 3 image after cropping in epsg:32718 and result image reprojected to epsg:3857, red polygon - roi.envelope:

Result image in the second case (reproject first):

[drnextgis]

If reprojecting first is not a solution you can use mask method of GeoRaster:

buffer = GeoVector(roi.get_shape(WEB_MERCATOR_CRS).buffer(100), WEB_MERCATOR_CRS)
raster.crop(buffer).reproject(dst_crs=WEB_MERCATOR_CRS).mask(roi).save(output_raster)

Example of output:

[astrojuanlu]

Thanks @drnextgis! I have a couple of questions:

1. I understand that crop and reproject are not meant to be commutative. Reprojecting a big image can be extremely expensive, especially when in the end you only want one crop. The buffer solution is nice, but a bit manual. Do you think it's worth it that telluric includes such a function?
2. Why mask instead of crop as the last step?

[drnextgis]

1. I think that it would be handy to have crop parameter for mask: raster.mask(vector, crop=True)
2. Let's see what happening in case of using mask (first image) and crop (second one) as the last step: