Convert polygons in shapefile to a NumPy array?

Question

I have a shapefile containing a number of polygons and all attributes being text. I'm trying to import the polygons into NumPy as an array where each polygon is represented as unique values.

I approach this by using gdal_rasterize to generate a GeoTIFF, which I then can convert to an array:

gdal_rasterize -a provcode -l longhurst PG:'host=localhost dbname=biomes' -tr 0.25 0.25 out.tif

and

tif = gdal.Open('out.tif')
tifArray = tif.ReadAsArray()

My problem is that all polygons get the same values since provcode is of type string. How can I make gdal_rasterize burn different values to represent different polygons?

Alternatively, is there a better way to do this conversion?

Answer 1

There is an option in GDAL to rasterize polygons based on their attribute. But as far as I know it can not be string. But you can just add an attribute to your features and then give each feature a unique id. Let's say we call this field ID.

Open your shapefile

source_ds = ogr.Open("Longhurst_world_v4_2010.shp")
source_layer = source_ds.GetLayer()

Create the destination raster data source

pixelWidth = pixelHeight = 1 # depending how fine you want your raster
x_min, x_max, y_min, y_max = source_layer.GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)
target_ds = gdal.GetDriverByName('GTiff').Create('temp.tif', cols, rows, 1, gdal.GDT_Byte) 
target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds.GetRasterBand(1)
NoData_value = 255
band.SetNoDataValue(NoData_value)
band.FlushCache()

Here is the important part. Instead of setting a general burn_value, use optionsand set it to the attribute that contains the relevant unique value ["ATTRIBUTE=ID"]

gdal.RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=ID"])  

Add a spatial reference

target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())
target_ds = None 

Now you have your shapefile as a raster and can read it with gdal.Open('temp.tif').ReadAsArray()