python gdal tutorial: reading raster data with gdal

GDAL natively supports more than 100 raster data types, covering all mainstream GIS and RS data formats, including

ArcInfo grids, ArcSDE raster, Imagine, Idrisi, ENVI, GRASS, GeoTIFF

HDF4, HDF5

USGS DOQ, USGS DEM

ECW, MrSID

TIFF, JPEG, JPEG2000, PNG, GIF, BMP

For a complete support list, please refer to http://www.gdal.org/formats_list.html

Import GDAL support library

Old version (Before 1.5): import gdal, gdalconst

New version (after 1.6): from osgeo import gdal, gdalconst

It is best to import both gdal and gdalconst, where the constants in gdalconst are prefixed, trying to be consistent with other modules Conflict is minimal. So you can directly import gdalconst like this: from osgeo.gdalconst import *

GDAL data driver, similar to OGR data driver, needs to create a certain type of data driver first, and then create the corresponding raster dataset.

Register all data drivers at once, but can only read but not write: gdal.AllRegister()

Register a certain type of data driver separately, so that you can read and write, and you can create a new data set:

driver = gdal.GetDriverByName('HFA')

driver.Register()

Open an existing raster dataset:

fn = 'aster.img'

ds = gdal.Open(fn, GA_ReadOnly)

if ds is None:

      print 'Could not open ' + fn

    sys.exit(1)

Read the number of pixels in the x direction, the number of pixels in the y direction, and the number of bands in the raster dataset

cols = ds. Raster A list that stores the geographical coordinate information of the raster dataset

adfGeoTransform[0] /* top left x x coordinate of the upper left corner*/

adfGeoTransform[1] /* w--e pixel resolution pixel resolution in the east-west direction Rate*/

adfGeoTransform[2] /* rotation, 0 if image is "north up" If the north is facing up, the rotation angle of the map*/

adfGeoTransform[3] /* top left y y coordinate of the upper left corner*/

adfGeoTransform[4] /* rotation, 0 if image is "north up" If the north is facing up, the rotation angle of the map*/

adfGeoTransform[5] /* n-s pixel resolution pixel resolution in the north-south direction*/

Note that the coordinates of raster datasets are generally based on the upper left corner.

The following example is to extract Geotransform from a raster dataset as a list, and then read the data in it

geotransform = ds.GetGeoTransform()

originX = geotransform[0]

originY = geotransform[3]originY = geotransform[3]

pixelWidth = geotransform[1]

pixelHeight = geotransform[5]

Calculate the relative position (pixel offset) of the pixel corresponding to a certain coordinate, which is the relative position of the coordinate and the pixel in the upper left corner. Calculated by the number of pixels, the calculation formula is as follows:

xOffset = int((x – originX) / pixelWidth)

yOffset = int((y – originY) / pixelHeight)

Reading the value of a certain pixel requires dividing Two steps

First read a band: GetRasterBand(), whose parameter is the index number of the band

Then use ReadAsArray(, , , < ;ysize>), read out the matrix starting from (xoff,yoff) and having size (xsize,ysize). If the matrix size is set to 1X1, one pixel is read. However, this method can only put the read data into the matrix, even if only one pixel is read. For example:

band = ds.GetRasterBand(1)

data = band.ReadAsArray(xOffset, yOffset, 1, 1)

If you want to read an entire image at once, then set offset to 0, size Then set it to the size of the entire image, for example:

data = band.ReadAsArray(0, 0, cols, rows)

But please note that to read the value of a certain pixel from data, you must use data[ yoff, xoff]. Be careful not to do it backwards. A matrix in mathematics is [row,col], but here it's just the opposite! Here row corresponds to the y-axis, and col corresponds to the x-axis.

Pay attention to releasing memory when appropriate, such as band = None or dataset = None. Especially when the image is very large

How to read raster data more efficiently? Obviously reading one by one is very inefficient. It is not a good idea to stuff the entire raster dataset into a two-dimensional array, because it still takes up a lot of memory. A better approach is to access data in blocks and put only the block you need into memory. This week's sample code has a utils module that can read the block size.

For example:

import utils

blockSize = utils.GetBlockSize(band)

xBlockSize = blockSize[0]

yBlockSize = blockSize[1]

Tile (tiled), that is, raster data is stored in blocks. Some formats, such as GeoTiff, are not tiled, and one line is a block. The Erdas imagine format is tiled by 64x64 pixels.

If a row is a block, then reading by row is more resource-saving.

If it is a tiled data structure, then setting the parameter value of ReadAsArray() so that it can only read one block at a time is the most efficient method. For example:

rows = 13, cols = 11, xBSize = 5, yBSize = 5

for i in range(0, rows, yBSize):

if i + yBSize < rows:

numRows = yBSize

else:

numRows = rows – i

for j in range(0, cols, xBSize):

if j + xBSize < cols:

numCols = xBSize

                                                                                                                  j

           data = band.ReadAsArray(j, i, numCols, numRows)

This code is universal and can be used from time to time.

The following introduces some two-dimensional array processing skills

Two libraries are used here, Numeric and numpy. Numeric is older and FWTools uses it. If you install and configure it yourself, you should use the more powerful numpy.

Data type conversion:

data = band.ReadAsArray(j, i, nCols, nRows)

data = data.astype(Numeric.Float) # Numeric

data = data.astype(numpy.float) # numpy

Or just write a simple sentence

data = band.ReadAsArray(j, i, nCols, nRows).astype(Numeric.Float)

mask

This is the function of Numeric and numpy libraries, input an array and condition, output a binary array. For example

mask = Numeric.greater(data, 0)mask = Numeric.greater(data, 0)

>>> a = Numeric.array([0, 4, 6, 0, 2])

>>> print a

[0 4 6 0 2]

>>> mask = Numeric.greater(a, 0)

>>> 1 1 0 1]

Array sum

>>> a = Numeric.array([0, 4, 6, 0, 2])

>>> print a>>> ; print a

[0 4 6 0 2]

>>> print Numeric.sum(a)

12

If it is a two-dimensional array, then sum will return a one-dimensional array

> ;>> b = Numeric.array([a, [5, 10, 0, 3, 0]])

>>> print b

[[ 0 4 6 0 2]

[ 5 10 0 3 0]]

>>> print Numeric.sum(b)>>> print Numeric.sum(b)

[ 5 14 6 3 2]

So, two-dimensional The sum of arrays is like this

>>> print Numeric.sum(Numeric.sum(b))

30

Here is a little trick, count the number of pixels greater than 0, you can use mask in combination and sum two functions

>>> print a

[0 4 6 0 2]

>>> mask = Numeric.greater(a, 0)

>>> print mask

[0 1 1 0 1]

>>> print Numeric.sum(mask)

3

The above is the python gdal tutorial: using gdal to read the content of raster data, more For related content, please pay attention to the PHP Chinese website (www.php.cn)!