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.
No comments:
Post a Comment