Cropping raster with GeoVector bigger than world bounds fails

[astrojuanlu]

Expected behavior and actual behavior

Cropping any raster with any GeoVector should work. Instead, if the GeoVector is bigger than the world, an error is raised.

simple_world.geojson.zip

Steps to reproduce the problem

In [1]: import telluric as tl

In [2]: rs_simple = tl.GeoRaster2.open("./tests/data/raster/creaf_gmap.tif")

In [3]: gv = tl.GeoVector.from_geojson("/tmp/simple_world.geojson")

In [4]: rs_simple.crop(gv)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-53bf33faf479> in <module>()
----> 1 rs_simple.crop(gv)

~/Development/telluric/telluric/georaster.py in crop(self, vector, resolution)
    710         :return: GeoRaster
    711         """
--> 712         bounds, window = self._vector_to_raster_bounds(vector)
    713         if resolution:
    714             xsize, ysize = self._resolution_to_output_shape(bounds, resolution)

~/Development/telluric/telluric/georaster.py in _vector_to_raster_bounds(self, vector, boundless)
    720     def _vector_to_raster_bounds(self, vector, boundless=False):
    721         # bounds = tuple(round(bb) for bb in self.to_raster(vector).bounds)
--> 722         window = self.window(*vector.get_shape(self.crs).bounds).round_offsets().round_shape(op='ceil')
    723         (ymin, ymax), (xmin, xmax) = window.toranges()
    724         bounds = (xmin, ymin, xmax, ymax)

~/.miniconda36/envs/telluric36/lib/python3.6/site-packages/rasterio/windows.py in round_offsets(self, op, pixel_precision)
    708             raise WindowError("operator must be 'ceil' or 'floor'")
    709         else:
--> 710             return Window(operator(round(self.col_off, pixel_precision)),
    711                           operator(round(self.row_off, pixel_precision)),
    712                           self.width, self.height)

ValueError: cannot convert float NaN to integer

Notice that with a simple .buffer(-1) the problem is solved in this case:

In [5]: gv.envelope
Out[5]: GeoVector(shape=POLYGON ((-180 -90, 180 -90, 180 83.70384706880299, -180 83.70384706880299, -180 -90)), crs=CRS({'init': 'epsg:4326'}))

In [6]: gv.buffer(-1).envelope  # slow
Out[6]: GeoVector(shape=POLYGON ((-179 -89, 179 -89, 179 82.70053994948313, -179 82.70053994948313, -179 -89)), crs=CRS({'init': 'epsg:4326'}))

In [7]: rs_simple.crop(gv.buffer(-1))
Out[7]: <telluric.georaster.GeoRaster2 at 0x7eff872bf630>

Operating system and installation details

Linux Mint 18.3, latest telluric.

[drnextgis]

It looks like this issue cannot be solved in the general case. For example, 180.0, 90.0 cannot be transformed from epsg:4326 into epsg:3857.

[astrojuanlu]

@drnextgis I agree. However, perhaps we could connect this to #43 in some way and produce a proper error when the conversion fails? Or, instead of an error, do a very small buffer and try to return something meaningful?

[AmitAronovitch]

We (me and @moshe742) looked at this, and it seems that the reason it does not work is because mercator is limited in latitude so cannot hold the vector. However, it is still possible to do the cropping if we represent everything in WGS84, and then convert back (Mecator must be able to hold the result, because the cropped raster is always smaller than the original one). The general algorithm we wrote is: if we fail converting the vector to raster's crs (catch the exception), then convert self to the vector.crs, crop the converted raster, then convert the result back to self.crs.

The only thing missing is a method for converting a raster to different CRS. We were not able to find an existing method, but this seems like a basic operation that we should support (see #75)