1

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')
Ari
  • 11
  • 4
  • 1
    Hi there! This question might be better suited for the [GIS Stack Exchange](https://gis.stackexchange.com/). But also, the code you shared in your post isn't reproducible because we don't have access to the `cis_SGRDAEA_20110201_pl_a.shp` file. You might want to consider adding a link to that somehow in your post. – Felipe D. Apr 04 '23 at 00:51
  • According to make_geocube's documentation, the documentation for the resolution argument says "Typically when using most CRSs, the first number would be negative." Have you tried using a negative value for the first component of resolution? https://corteva.github.io/geocube/html/geocube.html#make-geocube – Nick ODell Apr 04 '23 at 01:01
  • @ Nick ODell, yes, I already tried that. – Ari Apr 04 '23 at 01:07
  • You need to enable "access with the link" or else others can't access that file. – Nick ODell Apr 04 '23 at 01:08
  • @ Nick ODell Please try again now – Ari Apr 04 '23 at 01:10

1 Answers1

0

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.

Nick ODell
  • 15,465
  • 3
  • 32
  • 66
  • May I know where did you find CRS? Are you talking about the prj file? – Ari Apr 04 '23 at 02:01
  • @Ari The CRS of a geodataframe can be accessed by using `ds.crs`. It shows your CRS as WGS84. – Nick ODell Apr 04 '23 at 02:19
  • Thank you. For some reason, I need to use 5000 meters as the resolution. Can I change the CRS to Lambert first and then use the shared code? – Ari Apr 04 '23 at 02:28
  • See https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.to_crs.html – Nick ODell Apr 04 '23 at 03:24