2

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?

  • 1
    You could try to open the shapefile via the driver function. You can use : ogr.GetDriverByName("ESRI Shapefile").Open(vector_fn) – Matte Nov 21 '16 at 10:34
  • Does not work. Returns AttributeError: open. – FaCoffee Nov 21 '16 at 10:39
  • The problem is that the shapefile is long/lat projected and you are assigning a pixel_size = 25. – xunilk Nov 21 '16 at 11:35
4

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.

enter image description here

  • It worked! Can I visualize this raster from within Python itself, without any GIS software? – FaCoffee Nov 21 '16 at 11:48
  • 1
    I think that the answer is not but, it is not very difficult to do a new raster with array values. – xunilk Nov 21 '16 at 11:53

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.