Tuesday, 30 January 2018

python - Convert polygons in shapefile to a NumPy array?


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



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())


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


enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...