Python gdal Tutorial: Map algebra and the writing of raster data

Take the calculation of NDVI as an example:

Ndvi= (nir-red)/(nir+red)

Where Nir is band 3,red is band 2

The programming points are as follows:

1. Read Band 3 into the array data3 and read Band 2 into the array data2

2. The calculation formula is:

3. When both Data3 and data2 are 0 o'clock (for example, with 0 for NoData), an error that is removed by 0 causes the program to crash. You need to use mask with choose to remove the 0 value

The code is as follows, with only 4 lines

Data2 = Band2. Readasarray (0, 0, cols,rows). Astype (numeric.float16)

Data3 = Band3. Readasarray (0, 0, cols,rows). Astype (numeric.float16)

Mask = numeric.greater (data3 + data2, 0)

Ndvi = Numeric.choose (Mask, ( -99, (DATA3-DATA2)/(Data3 + data2)))

New Raster DataSet

Write the data you just calculated into a new raster dataset

The first thing to do is to copy a data-driven:

Driver = Indataset.getdriver ()

After you create a new data set

Create (<FILENAME>, <xsize>, <ysize>, [<bands>], [<gdaldatatype>])

Where the default value for bands is 1,gdaldatatype, the default type is Gdt_byte, for example

Outdataset = driver. Create (filename, cols, rows, 1, gdt_float32)

During the execution of this statement, the storage space has been allocated to the hard disk.

Before writing, you need to introduce a band object first

Outband = Outdataset.getrasterband (1)

Band object supports direct write matrix, two parameters are X-offset and y-offset, respectively

Outband.writearray (NDVI, 0, 0)

The following example summarizes this and the last block-by-write method

Xblocksize = 64

Yblocksize = 64

For I in range (0, rows, yblocksize):

If i + yblocksize < rows:

NumRows = Yblocksize

Else

NumRows = Rowsnumrows = Rows––ii

For j in range (0, cols, xblocksize):

If J + xblocksize < cols:

Numcols = Xblocksize

Else

Numcols = Cols–j

data = Band. Readasarray (J, I, Numcols, NumRows)

# do calculations here to create Outdata array

Outband.writearray (Outdata, J, i)

Band object can set NoData value

Outband.setnodatavalue (-99)

You can also read NoData values

ND = Outband.getnodatavalue ()

Calculate the band statistics

First write the cached data to disk with Flushcache ()

The statistics are then calculated using Getstatistics (<APPROX_OK> <force>). If Approx_ok=1 's calculations are based on pyramid, if force=0 then the statistics are not counted when the entire picture is to be reread again.

Outband.flushcache ()

Outband.getstatistics (0, 1)

Set a geographic reference point for a new diagram

If the new diagram is exactly the same as the geographic reference information for another diagram, it's easy.

Geotransform = Indataset.getgeotransform ()

Outdataset.setgeotransform (Geotransform)

proj = Indataset.getprojection ()

Outdataset.setprojection (proj)

Establish pyramids

Set the pyramids of Imagine style

Gdal. Setconfigoption (' hfa_use_rrd ', ' YES ')

Forcing the establishment of pyramids

Outdataset.buildoverviews (overviewlist=[2,4, 8,16,32,64,128])

Splicing of images

1. For each picture: Read the number of rows and columns, origin (minx,maxy), pixel length, pixel width, and calculate the coordinate range

MaxX1 = minX1 + (COLS1 * pixelwidth)

MinY1 = maxY1 + (ROWS1 * pixelheight)

2. Calculate the coordinate range of the output image:

MinX = min (minX1, minX2, ...) MaxX = max (maxX1, maxX2, ...)

miny = min (minY1, minY2, ...) Maxy = max (maxY1, maxY2, ...)

3. Calculate the number of rows and columns of the output image:

cols = Int ((maxx–minx)/pixelwidth)

rows = Int ((maxy–miny)/abs (Pixelheight)

4. Create and initialize the output image

5. For each picture to be spliced: Calculate the value of offset

XOffset1 = Int ((minx1-minx)/pixelwidth)

YOffset1 = Int ((maxy1-maxy)/pixelheight)

Read the data and write the offset as calculated above

6. For output Image: Calculate statistics, set geotransform: [MinX, Pixelwidth, 0, Maxy, 0, Pixelheight], set projection, build pyramids

The above is the Python gdal tutorial: Map algebra and raster data write content, more relevant content please pay attention to topic.alibabacloud.com (www.php.cn)!

