Monday, 23 September 2019

Python and GDAL: NoneType error while reading shapefile


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?



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.


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...