Python Array from Shapefile – How to Read Using GDAL and QGIS

gdalosgeopythonqgisshapefile

I am trying to read a shapefile as an array with Python but I keep getting the same error:

Warning 1: Failed to find field ID on layer shapefile_maido_tipe, skipping.

ERROR 4: `temp.tif' not recognized as a supported file format.

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 714, in runfile
    execfile(filename, namespace)
  File "C:\Anaconda2\lib\site-packages\spyderlib\widgets\externalshell\sitecustomize.py", line 74, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)
  File "C:/Users/python/Desktop/programme_python_tipe/test.py", line 38, in <module>
    print layer('shapefile_maido_tipe.shx')
  File "C:/Users/python/Desktop/programme_python_tipe/test.py", line 36, in layer
    return gdal.Open('temp.tif').ReadAsArray()
AttributeError: 'NoneType' object has no attribute 'ReadAsArray'

Here is the code that I got from Convert polygons in shapefile to a NumPy array? :

import gdal
from osgeo import osr
from osgeo import ogr

def layer(shapefile):

    # 1) opening the shapefile
    source_ds = ogr.Open(shapefile)
    source_layer = source_ds.GetLayer()

    # 2) Creating 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()

    # 4) 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'])


    # 5) Adding a spatial reference

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


    return gdal.Open('temp.tif').ReadAsArray()

print layer('shapefile_maido_tipe.shx')

# maido is the name of a mountain
# tipe is the name of a french school project

Can you tell what am I doing wrong please ?

In this link you will find the shapefile I use: https://drive.google.com/drive/folders/0B7P95aWmH4DUYUZfcTVvbzF4bG8?usp=sharing

PS: The Reunion Island National Park is the owner and also the creator of the shapefile. I have to mention them every time I communicate this file.

Best Answer

Your script was generally correct, but you didn't change the name of the attribute field you wanted to rasterize.

In the example you posted, you set ['ATTRIBUTE=ID'] as field, but it doesn't exists in your shapefile. You only have "Habitats" and "surface" as fields, so you need properly edit the code.

Therefore, you needed to edit the folders for both shapefile and rasterized layers, and the crs.

I slightly edited the code in this way:

import gdal
from osgeo import osr
from osgeo import ogr

def layer(shapefile):

    # 1) opening the shapefile
    source_ds = ogr.Open(shapefile)
    source_layer = source_ds.GetLayer()

    # 2) Creating 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(raster_path, 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()

    # 4) 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=surface'])

    # 5) Adding a spatial reference
    target_dsSRS = osr.SpatialReference()
    target_dsSRS.ImportFromEPSG(2975)
    target_ds.SetProjection(target_dsSRS.ExportToWkt())
    return gdal.Open(raster_path).ReadAsArray()


raster_path = 'C:/Users/path_to_the_rasterized_output/temp.tif'

shapefile = 'C:/Users/path_to_the_shapefile/shapefile_maido_tipe.shp'

print layer(shapefile)

and I think it is working by now because I obtain this rasterized layer (which overlaps the shapefile):

and this return from the print layer(shapefile) line (you see only '255' value because you set it as nodata value):

[[255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 ..., 
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]
 [255 255 255 ..., 255 255 255]]

Related Solutions

Python GDAL – Splitting Shapefile Per Feature

OK so a second attempt to answer your question with a pure GDAL solution.

Firstly, GDAL (Geospatial Data Abstraction Library) was originally just a library for working with raster geo-spatial data, while the separate OGR library was intended to work with vector data. However, the two libraries are now partially merged, and are generally downloaded and installed together under the combined name of GDAL. So the solution really falls under OGR. You have this in your initial code so I guess you knew this but it is an important distinction to remember when searching for tips and hints.

To read data from a vector layer, you're initial code is fine:

from osgeo import ogr
shapefile = ogr.Open(shapefile)
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    name = feature.GetField("NAME")
    geometry = feature.GetGeometryRef()
    print i, name, geometry.GetGeometryName()

We need to create a new feature before we can write it to a shapefile (or any other vector data set). To create a new feature, we first need: - A geometry - A feature definition, which will probably include field definitions Use the Geometry constructor ogr.Geometry() to create an empty Geometry object. Define what the geometry is in a different way for each type (point, line, polygon, etc.). So for example:

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10,20)

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10,10)
line.AddPoint(20,20)
line.SetPoint(0,30,30) #(10,10) -> (30,30)

For a field definition

fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)

Now you can create your vector layer. In this instance, a square polygon:

#create simple square polygon shapefile:
from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')

datasource = driver.CreateDataSource('YOUR_PATH')
layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbPolygon)

#create polygon object:
myRing = ogr.Geometry(type=ogr.wkbLinearRing)
myRing.AddPoint(0.0, 0.0)     #LowerLeft
myRing.AddPoint(0.0, 10.0)    #UpperLeft
myRing.AddPoint(10.0, 10.0)   #UpperRight
myRing.AddPoint(10.0, 0.0)    #Lower Right
myRing.AddPoint(0.0, 0.0)     #close ring
myPoly = ogr.Geometry(type=ogr.wkbPolygon)
myPoly.AddGeometry(myRing)
print ('Polygon area =',myPoly.GetArea())  #returns correct area of 100.0

#create feature object with point geometry type from layer object:
feature = ogr.Feature( layer.GetLayerDefn())
feature.SetGeometry(myPoly)
layer.CreateFeature(feature)

#flush memory - very important
feature.Destroy()
datasource.Destroy()

[GIS] Writing array to vrt

You can't write data to a VRT. It's a "Virtual Raster", and is actually just an XML file. Try something like:

gdal_translate -of VRT -b1 -b1 -b1 in.tif out.vrt

Then open the out.vrt in python and set the colour interpretation.

You can write the vrt xml in python instead of using the commandline, but it's easier to use the subprocess module to run the gdal_translate command.

Best Answer

Related Solutions

Python GDAL – Splitting Shapefile Per Feature

[GIS] Writing array to vrt

Related Question