python gdal tutorial: map algebra and writing of raster data

Take the calculation of NDVI as an example:

NDVI=(NIR-RED)/(NIR+RED)

where NIR is band 3 and RED is band 2

The programming points are as follows:

1. Read band 3 into the array data3, read band 2 into the array data2

2. The calculation formula is:

3. When data3 and data2 are both 0 (for example, 0 represents NODATA), a division by 0 error will occur, causing the program to crash. You need to use mask with choose to remove the 0 value

The code is as follows, there are 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)))

Create a new raster dataset

Write the data just calculated into the new raster dataset

First, copy a data driver:

driver = inDataset.GetDriver()

Create a new data set afterwards

Create(, , , [], [])

The default value of bands is 1, and the default value of GDALDataType is 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 also need to introduce the band object first

outBand = outDataset.GetRasterBand(1)

The band object supports direct writing to the matrix, the two parameters are x-direction offset and y-direction offset

outBand.WriteArray(ndvi, 0 , 0)

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

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

                                                                                                                                                               

            outBand.WriteArray(outData, j, i)

The band object can set the NoData value

outBand.SetNoDataValue(-99)

You can also read the NoData value

ND = outBand.GetNoDataValue()

Calculate the statistics of the band

First use FlushCache( )Write cached data to disk

and then use GetStatistics(, ) to calculate statistics. If approx_ok=1 then the calculation is based on pyramid, if force=0 then the statistics are not calculated when the entire graph has to be re-read.

outBand.FlushCache()

outBand.GetStatistics(0, 1)

Set the geographical reference point of the new image

If the new image has exactly the same geographical reference information as another image, it is very simple

geoTransform = inDataset.GetGeoTransform()

outDataset.SetGeoTransform(geoTransform)

proj = inDataset.GetProjection()

outDataset.SetProjection(proj)

Create pyramids

Set Imagine style pyramids

gdal.SetConfigOption ('HFA_USE_RRD', 'YES')

Force to build pyramids

outDataset.BuildOverviews(overviewlist=[2,4, 8,16,32,64,128])

Stitching of images

1.​​​​​​​​​​​​​​​​​​​​ 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 image to be stitched: Calculate the offset value

xOffset1 = int((minX1 - minX) / pixelWidth)

yOffset1 = int((maxY1 - maxY) / pixelHeight)

Read the data and write it according to the offset calculated above

6. For the output image: calculate statistics, set geotransform: [minX, pixelWidth, 0, maxY, 0, pixelHeight], set projection, and create pyramids

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