Pythongdal Tutorial: reading raster data with gdal

GDAL native supports over 100 grid data types, covering all mainstream GIS and RS data formats. GDAL native supports over 100 grid 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, MPEG-4, PNG, GIF, BMP

Complete support list can refer to http://www.gdal.org/formats_list.html

Import the 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. constants in gdalconst are prefixed and try to have the least conflict with other modules. Therefore, you can import gdalconst directly as follows: from osgeo. gdalconst import *

The GDAL data driver is similar to the OGR data driver. you must first create a data driver and then create a grid dataset for the response.

Register all data drivers at one time, but only read and write: gdal. AllRegister ()

Register a data driver for a specific type. in this way, you can read or write data and create a dataset:

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)

Reads the number of records in the x direction, the number of records in the y direction, and the number of bands in the raster dataset.

Cols = ds. RasterXSize

Rows = ds. RasterYSize

Bands = ds. RasterCount

Note that there are no parentheses behind them, because they are properties and not methods)

Read the geographic coordinate reference information (georeference info)

GeoTransform is a list that stores the geographic coordinate information of the raster dataset.

AdfGeoTransform [0]/* top left x coordinate in the upper left corner */

AdfGeoTransform [1]/* w -- e pixel resolution east-west pixel resolution */

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

AdfGeoTransform [3]/* top left y coordinate in the upper left corner */

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

AdfGeoTransform [5]/* n-s pixel resolution up and down */

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

The following example extracts Geotransform from a grid data set as a list and then reads the data.

Geotransform = ds. GetGeoTransform ()

Origform = 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 coordinate, that is, the relative position of the coordinate and the pixel in the upper left corner. the formula is as follows:

XOffset = int (x-origwidth)/pixelWidth)

YOffset = int (y-originY)/pixelHeight)

Read the value of a certain pixel in two steps

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

Then use ReadAsArray ( , , , ), Read the matrix of the size (xsize, ysize) starting from (xoff, yoff. If the matrix size is set to 1X1, a 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 a whole graph at a time, set offset to 0 and size to the size of the entire graph. for example:

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

However, you must use data [yoff, xoff] to read the value of a certain pixel from data. Do not reverse it. In mathematics, the matrix is [row, col], which is the opposite! The row corresponds to the y axis, and the col corresponds to the x axis.

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

How can we read raster data more efficiently? Obviously, reading efficiency is very low one by one. it is not a good way to insert the entire raster dataset into a two-dimensional array, because there is still a lot of memory. A better way is to access data by block and put the block to the memory. This week's sample code contains an utils module that can read the block size.

For example:

Import utils

BlockSize = utils. GetBlockSize (band)

XBlockSize = blockSize [0]

YBlockSize = blockSize [1]

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

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

If it is a tiled data structure, it is the most efficient method to set the ReadAsArray () parameter value so that it can read only into a block at a time. 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

Else:

NumCols = colsnumCols = cols-j

Data = band. ReadAsArray (j, I, numCols, numRows)

This piece of code is universal and can be used frequently.

The following describes how to process a 1.2-dimension array.

Two libraries, Numeric and numpy, are used here. Numeric is old, and FWTools is used. If you install and configure it yourself, it will still have a more powerful numpy configuration.

Data type conversion:

Data = band. ReadAsArray (j, I, nCols, nRows)

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

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

Or simply write only one 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 a condition and 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

[0 4 6 0 2]

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

>>> Print mask

[0 1 1 0 1]

Array summation

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

>>> Print

[0 4 6 0 2]

>>> Print Numeric. sum ()

12

For a two-dimensional array, sum returns 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]

Therefore, the sum of two-dimensional arrays is like this.

>>> Print Numeric. sum (Numeric. sum (B ))

30

Here is a trick to count the number of pixels greater than 0. you can use the mask and sum functions together.

>>> Print

[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: use gdal to read the content of the grid data. For more information, see PHP Chinese website (www.php1.cn )!