Why rasterizing a shapefile in python using make_geocube gives an empty geotiff file?

513 Views Asked by At

When I try to rasterize a shapefile (ice chart) for attribute CA, it does not work. It gives me a GeoTIFF file with a size of 1 KB. When I open it on QGIS, nothing shows up because it is empty.

Here is the Warning: C:\Users\ccw\anaconda3\envs\geo_env\lib\site-packages\rasterio_init_.py:230: NotGeoreferencedWarning: The given matrix is equal to Affine.identity or its flipped counterpart. GDAL may ignore this matrix and save no geotransform without raising an error. This behavior is somewhat driver-specific. s = writer(path, mode, driver=driver,

Here is the link to data: https://drive.google.com/file/d/12gavTn0rrIXnzVKXbW-dcKhRVzDnVhX1/view?usp=share_link

I've also used this prj file: https://drive.google.com/file/d/1JQqcdz-W3OuGbqs0jNi2jLu-RcZTarSc/view?usp=share_link

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

ds = gpd.read_file('cis_SGRDAEA_20110201_pl_a.shp')
ds['CA'] = ds['CA'].astype(float)
grid = make_geocube(vector_data=ds, measurements=['CA'], resolution=(5000,-5000))
grid.CA.rio.to_raster('test.tif')
1

There are 1 best solutions below

4
Nick ODell On

This CRS is in units of degrees. The code is asking for tiles 5000 degrees by 5000 degrees. Since it only needs one of those to cover the whole map, that's all it uses. It then produces a 1x1 geotiff.

Changing it to this:

grid = make_geocube(vector_data=ds, measurements=['CA'], resolution=(1,-1))

causes it to create a larger 61x16 raster output.