Python – Fixing NoneType Error While Reading Shapefile with GDAL

gdalgdal-rasterizepythonshapefile

I am a newbie to GDAL. I have a shapefile of Britain which I want to rasterize. This file can be obtained from here selecting United Kingdom as Country and Administrative areas as Subject.

Following the GDAL tutorial ("Convert vector layer to array"), I have this script:

import pyproj
import osgeo
import numpy as np
from osgeo import gdal
import os 
import ogr, gdal
from gdalconst import *

#The dir and file to work with
mydir=r'C:\Users\GBR_adm'
os.chdir(mydir)
vector_fn = 'GBR_adm0.shp'

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

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
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('MEM').Create('', x_res, y_res, 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=[1])

# Read as array
array = band.ReadAsArray()
print array

This script raises an error: AttributeError: 'NoneType' object has no attribute 'SetGeoTransform', and this relates to type(target_ds) being NoneType. I have read somewhere that this means that GDAL has failed to load the shapefile. Why does this happen, and how to avoid it?

Best Answer

The problem is that the shapefile is long/lat projected and you are assigning a pixel_size = 25. I tried your code (slightly modified for adapting to my system) by using (e.g.) a pixel_size = 0.0025 and, it ran without any errors.

import pyproj
import osgeo
import numpy as np
from osgeo import gdal
from osgeo import ogr
import os 
from gdalconst import *

#The dir and file to work with
mydir=r'/home/zeito/pyqgis_data/GBR_adm'
os.chdir(mydir)
vector_fn = 'GBR_adm0.shp'

# Define pixel_size and NoData value of new raster
pixel_size = 0.0025
NoData_value = 255

# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
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('MEM').Create('', x_res, y_res, 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=[1])

source_ds = None

# Read as array
array = band.ReadAsArray()
print "min:", array.min(), "max;", array.max()
print "rows:", len(array), "columns:", len(array[0])

After running it, it can be observed at the Python Console of QGIS the min and max values and rows and columns of this array.

Related Solutions

Python GDAL/OGR – Handling NoneType Feature in Python GDAL/OGR

It is a pure Python problem and not an GDAL/OGR problem. But if you want to program in Python,you need to understand the difference between an iterator and an iterable (many, many references on the Web, but one of the best is Loop Like A Native by Ned Batchelder). These are fundamental concepts in Python.

"The important operation on an iterable is iter(), which will return an iterator. And the only operation available on an iterator is next(), which either returns the next value, or raises StopIteration, a special exception that means the iteration is finished." (from Ned Batchelder)

An iterable produces an iterator and an iterator produces a stream of values.

So what's the difference between the for loop and pointsLayer.GetNextFeature() ?

the for loops = iterator and pointsLayer= iterable

for feature in iterable:
   statements

so with your example:

from osgeo import ogr
source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
for feature in pointsLayer:
    geom =feature.GetGeometryRef()
    xy = geom.GetPoint()
    print xy

(272022.68669955182, 155404.12013724342, 0.0)
(272904.99338241993, 152881.6706538822, 0.0)
.....
(272796.14718989708, 152075.00336062044, 0.0)

But we can create another type of iterator with iter() and next():

source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
iterator = iter(pointsLayer)
# first feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272022.68669955182, 155404.12013724342, 0.0)
# second feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272904.99338241993, 152881.6706538822, 0.0)
....
# end raises StopIteration error to signal that iteration is complete
feature = iterator.next()
Traceback (most recent call last):
...
StopIteration

Unlike the for loop that handles the StopIteration error, you need here another control structure (while loop,if,...)

Thus, what is pointsLayer.GetNextFeature() ?

It is an iterator and you can replace iterator = iter(pointsLayer) and iterator.next() by

feature = pointsLayers.GetNextFeature()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
.....
feature = pointsLayers.GetNextFeature()
....

I hope that I have been able to explain why you should not mix the for loops (iterator) with .GetNextFeature() (another iterator).

GDAL Python – Suppress GDAL Error Messages in Python

There are some python 'only' rules stated at the 'Python Gotchas' page.

At the beginning they state that:

By default, the GDAL and OGR Python bindings do not raise exceptions when errors occur. Instead they return an error value such as None and write an error message to sys.stdout.

So first enable the use of exceptions, by issuing gdal.UseExceptions() somewhere in the beginning of your script. And secondly catch if any exceptions and do whatever you want with them (including nothing)

gdal.UseExceptions()  # Enable errors
...
try:
    band.GetStatisticsts()  
except RuntimeError:  # <- Check first what exception is being thrown
    pass

If you're curious here's a link from python's wiki where they describe the methods of exception handling.

As far as I can tell, using gdal.UseExceptions() will normalize the behaviour of all the methods in the gdal lib, by making those methods to use the Python Exceptions. If you really want to ignore all exceptions put the problematic part of your code inside a try/except block eg.:

>>>
try:
    print "1"
    print "2"
    a = 0
    b = 1
    c = b/a   # <- I just divided by zero.
except:
    print "Everything is ok" # <-  but Everything is Ok. 
                             # If you really want to be silent replace with pass
>>> 1
    2
    'Nothing to see here'

Just be careful, because that way when an exception happens, your script will silently break from the procedure without any indication of what went wrong.

Also I want to point a couple of things as well:

If you're opt to go that way, you're responsible if your program does. You're choosing to turn a blind eye at any errors that come in your way.
If you catch ALL the exceptions, you will also catch the 'KeyboardInterrupt' exception as well. No more ctrl-c (if you're inside a try/except block).
Your code at some point is encountering a processing error. Since you're not using gdal.UseExceptions() you're experiencing that error as a sys.stdout response. Use a try/except block at that point to plug that hole.

Best Answer

Related Solutions

Python GDAL/OGR – Handling NoneType Feature in Python GDAL/OGR

GDAL Python – Suppress GDAL Error Messages in Python

Related Question