Home Backend Development Python Tutorial python gdal tutorial: map algebra and writing of raster data

python gdal tutorial: map algebra and writing of raster data

Dec 24, 2016 pm 05:08 PM
gdal

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)!


Statement of this Website
The content of this article is voluntarily contributed by netizens, and the copyright belongs to the original author. This site does not assume corresponding legal responsibility. If you find any content suspected of plagiarism or infringement, please contact admin@php.cn

Hot AI Tools

Undresser.AI Undress

Undresser.AI Undress

AI-powered app for creating realistic nude photos

AI Clothes Remover

AI Clothes Remover

Online AI tool for removing clothes from photos.

Undress AI Tool

Undress AI Tool

Undress images for free

Clothoff.io

Clothoff.io

AI clothes remover

Video Face Swap

Video Face Swap

Swap faces in any video effortlessly with our completely free AI face swap tool!

Hot Tools

Notepad++7.3.1

Notepad++7.3.1

Easy-to-use and free code editor

SublimeText3 Chinese version

SublimeText3 Chinese version

Chinese version, very easy to use

Zend Studio 13.0.1

Zend Studio 13.0.1

Powerful PHP integrated development environment

Dreamweaver CS6

Dreamweaver CS6

Visual web development tools

SublimeText3 Mac version

SublimeText3 Mac version

God-level code editing software (SublimeText3)

How to solve the permissions problem encountered when viewing Python version in Linux terminal? How to solve the permissions problem encountered when viewing Python version in Linux terminal? Apr 01, 2025 pm 05:09 PM

Solution to permission issues when viewing Python version in Linux terminal When you try to view Python version in Linux terminal, enter python...

How to avoid being detected by the browser when using Fiddler Everywhere for man-in-the-middle reading? How to avoid being detected by the browser when using Fiddler Everywhere for man-in-the-middle reading? Apr 02, 2025 am 07:15 AM

How to avoid being detected when using FiddlerEverywhere for man-in-the-middle readings When you use FiddlerEverywhere...

How to efficiently copy the entire column of one DataFrame into another DataFrame with different structures in Python? How to efficiently copy the entire column of one DataFrame into another DataFrame with different structures in Python? Apr 01, 2025 pm 11:15 PM

When using Python's pandas library, how to copy whole columns between two DataFrames with different structures is a common problem. Suppose we have two Dats...

How to teach computer novice programming basics in project and problem-driven methods within 10 hours? How to teach computer novice programming basics in project and problem-driven methods within 10 hours? Apr 02, 2025 am 07:18 AM

How to teach computer novice programming basics within 10 hours? If you only have 10 hours to teach computer novice some programming knowledge, what would you choose to teach...

How does Uvicorn continuously listen for HTTP requests without serving_forever()? How does Uvicorn continuously listen for HTTP requests without serving_forever()? Apr 01, 2025 pm 10:51 PM

How does Uvicorn continuously listen for HTTP requests? Uvicorn is a lightweight web server based on ASGI. One of its core functions is to listen for HTTP requests and proceed...

How to solve permission issues when using python --version command in Linux terminal? How to solve permission issues when using python --version command in Linux terminal? Apr 02, 2025 am 06:36 AM

Using python in Linux terminal...

How to get news data bypassing Investing.com's anti-crawler mechanism? How to get news data bypassing Investing.com's anti-crawler mechanism? Apr 02, 2025 am 07:03 AM

Understanding the anti-crawling strategy of Investing.com Many people often try to crawl news data from Investing.com (https://cn.investing.com/news/latest-news)...

See all articles