Python gdal Tutorial: Reading raster data with Gdal

Gdal native supports more than 100 raster data types, covering all major 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

The complete support list can be consulted http://www.gdal.org/formats_list.html

Import Gdal Support Library

Older version (1.5 ago): Import Gdal, Gdalconst

New version (after 1.6): From OSGeo import Gdal, Gdalconst

Gdal and Gdalconst are best to import, where the constants in Gdalconst are prefixed, trying to minimize collisions with other module. So for gdalconst you can import this directly: from osgeo.gdalconst import *

Gdal data driven, similar to OGR data driven, you need to create a type of data-driven, and then create a raster dataset of the response.

One-time registration of all data-driven, but can only read can not write: Gdal. Allregister ()

Register a type of data driver separately so that it can be read or write, and you can create a new dataset:

Driver = Gdal. Getdriverbyname (' HFA ')

Driver. Register ()

To 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 X-directional pixels of the raster dataset, the number of pixels in the Y-direction, and how many bands

cols = ds. Rasterxsize

rows = ds. Rasterysize

Bands = ds. Rastercount

Note that there are no parentheses behind because they are attributes (properties) that are not methods (methods)

Read geo-coordinate reference information (georeference info)

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

ADFGEOTRANSFORM[0]/* Top left x x coordinate */

ADFGEOTRANSFORM[1]/* W--E pixel resolution East-West pixel resolution */

ADFGEOTRANSFORM[2]/* rotation, 0 if image is ' north up ' if the northern side is upward, the map's rotation angle */

ADFGEOTRANSFORM[3]/* Top left y y coordinate */

ADFGEOTRANSFORM[4]/* rotation, 0 if image is ' north up ' if the northern side is upward, the map's rotation angle */

ADFGEOTRANSFORM[5]/* N-S pixel resolution pixel resolution in North/south direction */

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

The following example removes Geotransform from a raster dataset as a list and then reads the data from it

Geotransform = ds. Getgeotransform ()

Originx = geotransform[0]

Originy = Geotransform[3]originy = Geotransform[3]

Pixelwidth = geotransform[1]

Pixelheight = geotransform[5]

Calculates the relative position (pixel offset) of the corresponding pixel of a coordinate, which is the relative position of the pixel in the upper-left corner, calculated by the number of pixels, as follows:

Xoffset = Int ((x–originx)/pixelwidth)

Yoffset = Int ((y–originy)/pixelheight)

It takes two steps to read the value of a pixel.

First read a band (band): Getrasterband (<index>), whose parameter is the index number of the band

Then using Readasarray (<XOFF>, <yoff>, <xsize>, <ysize>), read the matrix from (Xoff,yoff) to the size (xsize,ysize). If you set the matrix size to 1x1, you are reading one pixel. 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 picture at once, set offset to 0,size to the size of the entire frame, for example:

data = Band. Readasarray (0, 0, cols, rows)

Note, however, that reading the value of a pixel from data must be done with Data[yoff, Xoff]. Be careful not to go against the wrong. The matrix in mathematics is [Row,col], and here is the opposite! The row corresponds to the y-axis, and col corresponds to the x-axis.

Note free memory at the appropriate time, such as Band = None or DataSet = None. Especially when it's a big picture.

How can I read raster data more efficiently? Obviously one of the reading efficiency is very low, it is not a good idea to plug the entire raster dataset into a two-dimensional array, because this occupies a lot of memory. The better way is to access the data by block (blocks), just put the piece that you want to use in 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]

Tiling (tiled), where raster data is stored by block. Some formats, such as GeoTiff, are not tiled, and one row is a block. The Erdas imagine format is paved with 64x64 like Suping.

If a row is a block, reading by rows is more resource-efficient.

If it is a tiled data structure, then setting the parameter value of Readasarray () allows it to read into only one block at a time, which is the most efficient method. For example:

rows = cols = one, 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 at all times.

Here are the techniques for handling 1.2-D arrays

There are two libraries, numeric and numpy. Numeric older, Fwtools used it. If you install the configuration, you will be equipped with a more powerful numpy.

Data type conversions:

data = Band. Readasarray (J, I, NCols, nRows)

data = Data.astype (numeric.float) # Numeric

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

Or simply write one sentence.

data = Band. Readasarray (J, I, NCols, nRows). Astype (Numeric.float)

Masking mask

This is the function of the numeric and numpy libraries, enter an array and 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 a

[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 a>>> Print a

[0 4 6 0 2]

>>> print Numeric.sum (a)

12

If it is a two-dimensional array, then 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]

So, the summation of the two-dimensional array is going to be

>>> Print numeric.sum (numeric.sum (b))

30

Here's a tip to count the number of pixels greater than 0, and you can use mask and sum two functions together

>>> 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: Read the contents of raster data with Gdal, more relevant content please pay attention to topic.alibabacloud.com (www.php.cn)!

