2

I have the following polygons shapefile consisting of polygons of the five boroughs of New York City from "NYC Open Data", linked here: https://data.cityofnewyork.us/City-Government/Borough-Boundaries/tqmj-j8zm

I want to take this shapefile and rasterize it to a .tif file using gdal in Python.

I am using the following python script example to accomplish this task, cited from "Python GDAL/OGR Cookbook 1.0 documentation", linked here: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html

This is what I am using:

from osgeo import gdal, ogr

gdal.UseExceptions()

# Define pixel_size and NoData value of new raster
pixel_size = 25
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp'

# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

However, this returns the following error message:

RuntimeError                              Traceback (most recent call last)
<ipython-input-32-0e6d6c7404d4> in <module>
     21 x_res = int((x_max - x_min) / pixel_size)
     22 y_res = int((y_max - y_min) / pixel_size)
---> 23 target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)
     24 target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
     25 band = target_ds.GetRasterBand(1)

~\Anaconda3\lib\site-packages\osgeo\gdal.py in Create(self, *args, **kwargs)
   1892     def Create(self, *args, **kwargs):
   1893         """Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None) -> Dataset"""
-> 1894         return _gdal.Driver_Create(self, *args, **kwargs)
   1895 
   1896 

RuntimeError: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

How do I address this error so that the raster .tif file is correctly created and sent to my directory folder? I am not sure why this message is telling me I am creating a 0x0 raster file or what is wrong with the raster I am creating, since it seems like all I have to do to run this code is simply enter my input shapefile, which it seems like I did correctly.

UPDATE (with second attempt at code, incorporating pixel size and setting the projection):

from osgeo import gdal, ogr

# Define NoData value of new raster
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp'

# Filename of the raster Tiff that will be created
raster_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/test2.tif'

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
proj = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, 100, 100, 1, gdal.GDT_Byte)
target_ds.SetProjection(proj.ExportToWkt())
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

This code runs and outputs a 0. I am not sure what this means. I then go to open test2.tif in QGIS, and all I see is a black square with all 0 values, not the rasterized polygon shapes I expect.

4
  • Try enabling exceptions and seeing if you get a more helpful error
    – mikewatt
    Nov 24 '21 at 20:57
  • Good suggestion, thanks. I used gdal.UseExceptions() and pasted the new error traceback message in my post. The error is telling my my 0x0 dataset is illegal, and I am not sure what this means, but perhaps it could be an error with my projection. I am trying to figure out what this means, or which step resulted in this problem raster. Nov 24 '21 at 21:07
  • Good point there. I had tried changing the burn value there to 1, but this resulted in a raster that just a black square with all 0 value. I tried your gdal.Rasterize() suggestion and it works, and produces the .tif raster like I want, but when I open the file in QGIS, it is as well just a black square, with all 0 value pixels. I tried leaving xRes and yRes blank, but that didn't work. I then tried setting them both to 100, but then the black square was far beyond the extent of the NYC polygons. I am just trying to get the output raster to match the size/extent of the input shapefile. Nov 25 '21 at 19:42
  • Success! This worked! Thank you! I see that the out-of-scale issue is caused by a discrepancy between degrees and lon/lat, but how did you know to choose exactly 0.0001? Sorry I have always been a bit confused between pixel size and resolution. Is 0.0001 the one "correct" size value or just the best approximation? I am just trying to understand so I know how to set the most appropriate pixel size next time. Nov 25 '21 at 23:25
2

https://github.com/corteva/geocube

from geocube.api.core import make_geocube

vector_fn = 'C:/Users/MyName/Desktop/NYC_GIS/Boroughs_Test/geo_export_d92842a7-2631-4fe5-ab83-ee070743980c.shp'

out_grid = make_geocube(
    vector_data=vector_fn,
    measurements=["column_name"],
    resolution=(-25, 25),
    fill=-9999,
)
out_grid["column_name"].rio.to_raster("my_rasterized_column.tif")
10
  • I really really want to use Geocube as my main tool for raster processing! However, I just cannot seem to get it to download properly. When I try to simply run pip install geocube from Anaconda Prompt (Anaconda3), I get a a long list of error messages (much too long to post the complete traceback here) including: WARNING: Ignoring invalid distribution -umpy (c:\users\MyName\anaconda3\lib\site-packages), Building wheel for gdal (setup.py) ... error ERROR: Command errored out with exit status 1:, andERROR: Failed building wheel for gdal. I am not sure how to address this. Nov 25 '21 at 0:37
  • For geocube, you should use conda to install from conda-forge.
    – snowman2
    Nov 25 '21 at 0:41
  • I tried using conda-forge via those instructions and received the messages: Solving environment: failed with initial frozen solve. Retrying with flexible solve. and Solving environment: failed with repodata from current_repodata.json, will retry with next repodata source., ultimately being stuck on `Solving environment: `. I am not sure how to fix this, but maybe it is something to do with how I set up my system. Nov 25 '21 at 1:15
  • Are you creating a new environment? See the conda install section here: pyproj4.github.io/pyproj/stable/installation.html
    – snowman2
    Nov 25 '21 at 2:11
0

The parameters you have entered into the Create() function are incorrect. You can find them using the python help function:

from osgeo import gdal
driver = gdal.GetDriverByName('GTiff')
help(driver.Create)
Create(Driver self, char const * utf8_path, int xsize, int ysize, int bands=1, GDALDataType eType, char ** options=None)

In your case, you have entered the pixel resolution (xres, yres) inplace of the raster dimensions (xsize, ysize):

target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, 1, gdal.GDT_Byte)

To avoid having to manually calculate the raster dimensions, you could rasterise using a reference image as illustrated here.

You have also failed to define a projection system. Consequently, the output image will not be georeferenced. You can do this using the SetProjection() function as follows:

vector = ogr.Open(vector_fn)
layer = vector.GetLayer()
proj = layer.GetSpatialRef()
target_ds.SetProjection(proj.ExportToWkt())

Your burn value is set to zero. You should probably change this to 1:

gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])
7
  • Ok, I see that I need actual values for xres and yres. I just replaced them with 100 and 100, and that part of the code runs now. Though I am having some trouble with the SetProject() component. Nov 24 '21 at 22:47
  • Particularly I am not sure where in the script to place the function to set the projection. The code runs, but the output .tif file is just a blank 0 value square when I open in QGIS. Nov 24 '21 at 22:55
  • To reiterate, the Create() function requires the number of pixels along each dimension (xsize, ysize) NOT the resolution (xres, yres). The SetProjection() function should be applied before you rasterize.
    – Osian
    Nov 24 '21 at 22:58
  • Understood, I got rid of xres and yres and just placed 100 and 100 for the xsize and ysize to address the number of pixels. I then placed target_ds.SetProjection(prok.ExportToWkt()) after proj = layer.GetSpatialRef() and before I use .RasterizeLayer(), but still the output .tif is all blank 0 value. I pasted my latest attempt in my post to show where I used your suggestions. Nov 24 '21 at 23:13
  • I see the problem. Your burn value is zero. Change this to 1 and it should be OK.
    – Osian
    Nov 24 '21 at 23:20
0

You are telling GDAL to burn all 0 values with burn_values=[0].

Additionally, setting xRes and yRes to 100 won't work as your data uses lon, lat coords (1 degree is 111,000m at the equator so 100 degrees is huuuuge).

By the way, the newer gdal.Rasterize is easier to use than the old gdal.RasterizeLayer.

Particularly there's no need to pre-create the output dataset, just pass an output filename. For example:

from osgeo import gdal

# Define NoData value of new raster
NoData_value = -9999

# Filename of input OGR file
vector_fn = 'Boroughs.geojson'

# Filename of the raster Tiff that will be created
raster_fn = 'Boroughs.tif'

# Open the data source and read in the extent
source_ds = gdal.OpenEx(vector_fn)
pixel_size = 0.001  # about 11 metres(ish) use 0.001 if you want roughly 100m 

gdal.Rasterize(raster_fn, source_ds, format='GTIFF', outputType=gdal.GDT_Byte, creationOptions=["COMPRESS=DEFLATE"], noData=NoData_value, initValues=NoData_value, xRes=pixel_size, yRes=-pixel_size, allTouched=True, burnValues=1)
1
  • Apologies this took so long, but I tried this out, and it is working, but it is just taking so incredibly long that I had to end the code run, and I am not sure if that means there is something wrong with my file. A raster file is produced still, but it is 453 megabytes big, which seems a bit much for my simple NYC example, and it crashes QGIS when I try to open it. Jan 10 at 21:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.