corteva / geocube

Tool to convert geopandas vector data into rasterized xarray data.
https://corteva.github.io/geocube
BSD 3-Clause "New" or "Revised" License
344 stars 27 forks source link

Geocube doesn't load gpd Dataframe #118

Closed Stoofik closed 2 years ago

Stoofik commented 2 years ago

Code Sample, a copy-pastable example if possible

from geocube.api.core import make_geocube
import geopandas as gpd

# load
path = R"path.gpkg"
shp = gpd.read_file(path)

# create geocube
blds_cube = make_geocube(
    vector_data=shp,
    measurements=["to_tif"],
    output_crs=3857,
    resolution=(-2, 2),
    fill=0
)

Problem description

Geocube doesn't work... Ends with error: pyproj.exceptions.CRSError: Expect string or any object with .to_epsg() or .to_wkt() methods

Expected Output

Geocube works as expected. No changes to previously working code have been made

Environment Information

Installation method

snowman2 commented 2 years ago

Are you able to provide the input file to re-produce this issue?

Stoofik commented 2 years ago

Are you able to provide the input file to re-produce this issue?

I can recreate this behavior by simply drawing random polygons in QGIS into blank new file and use them as path.gpkg. I tried this on different machines and it does the same.

For example this simple file

rnp.zip

snowman2 commented 2 years ago

If you change output_crs="EPSG:3857" it works.

snowman2 commented 2 years ago

Full error when using the integer version 3857:

---------------------------------------------------------------------------
CRSError                                  Traceback (most recent call last)
Input In [3], in <cell line: 9>()
      6 shp = gpd.read_file(path)
      8 # create geocube
----> 9 blds_cube = make_geocube(
     10     vector_data=shp,
     11     measurements=["to_tif"],
     12     output_crs=3857,
     13     resolution=(-2, 2),
     14     fill=0
     15 )

File .../lib/python3.10/site-packages/geocube/api/core.py:99, in make_geocube(vector_data, measurements, datetime_measurements, output_crs, resolution, align, geom, like, fill, group_by, interpolate_na_method, categorical_enums, rasterize_function)
     35 """
     36 Rasterize vector data into an ``xarray`` object.  Each attribute will be a data
     37 variable in the :class:`xarray.Dataset`.
   (...)
     95 
     96 """
     97 geobox_maker = GeoBoxMaker(output_crs, resolution, align, geom, like)
---> 99 return VectorToCube(
    100     vector_data=vector_data,
    101     geobox_maker=geobox_maker,
    102     fill=fill,
    103     categorical_enums=categorical_enums,
    104 ).make_geocube(
    105     measurements=measurements,
    106     datetime_measurements=datetime_measurements,
    107     group_by=group_by,
    108     interpolate_na_method=interpolate_na_method,
    109     rasterize_function=rasterize_function,
    110 )

File .../lib/python3.10/site-packages/geocube/vector_to_cube.py:80, in VectorToCube.__init__(self, vector_data, geobox_maker, fill, categorical_enums)
     60 """
     61 Initialize the GeoCube class.
     62 
   (...)
     77 
     78 """
     79 self._vector_data = load_vector_data(vector_data)
---> 80 self._geobox = geobox_maker.from_vector(self._vector_data)
     81 self._grid_coords = affine_to_coords(
     82     self._geobox.affine, self._geobox.width, self._geobox.height
     83 )
     84 self._fill = fill if fill is not None else numpy.nan

File .../lib/python3.10/site-packages/geocube/geo_utils/geobox.py:172, in GeoBoxMaker.from_vector(self, vector_data)
    169     raise RuntimeError("Must specify 'resolution' if 'like' not specified.")
    171 if self.output_crs:
--> 172     crs = CRS(self.output_crs)
    173 else:
    174     crs = CRS(vector_data.crs)

File ~/miniconda/envs/midas/lib/python3.10/site-packages/odc/geo/crs.py:98, in CRS.__init__(self, crs_spec)
     95     self._crs, self._str, self._epsg = _make_crs(wkt)
     96     return
---> 98 raise CRSError(
     99     "Expect string or any object with `.to_epsg()` or `.to_wkt()` methods"
    100 )

CRSError: Expect string or any object with `.to_epsg()` or `.to_wkt()` methods
Stoofik commented 2 years ago

Can't believe it... I just ran the same script as I used before, although I am using shared scripts and machines... It still makes little sense. I am very sorry I didn't fully investigate.

snowman2 commented 2 years ago

No worries. Thanks for the report. There is a fix coming upstream. See: https://github.com/opendatacube/odc-geo/pull/64