【GDAL学习】地图代数与栅格数据的写入

1.Assignment 5a：

• Create an NDVI image
• Read in data from aster.img
• Create an NDVI image
• Write out NDVI to new file
• Can do entire image at once or block by block 
• Don't forget to calculate statistics, set projection and georeferencing information, and build pyramids

# Assignment 5a: Create an NDVI image
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys

os.chdir('E:/data/GDAL/ospy_data5')

gdal.AllRegister()

ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:
    print('Cannot open aster.img')
    sys.exit(1)

rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

band2 = ds.GetRasterBand(2)
band3 = ds.GetRasterBand(3)

xBlockSize = 64
yBlockSize = 64

driver = ds.GetDriver()
ndviDS = driver.Create('ndvi.img', cols, rows, 1, GDT_Float32)
if ndviDS is None:
    print('Cannot open ndvi.img')
    sys.exit(1)
ndviBand = ndviDS.GetRasterBand(1)

for i in range(0, rows, yBlockSize):
    if (i + yBlockSize) < rows:
        ySize = yBlockSize
    else:
        ySize = rows - i

    for j in range(0, cols, xBlockSize):
        if (j + xBlockSize) < cols:
            xSize = xBlockSize
        else:
            xSize = cols - j

        data2 = band2.ReadAsArray(j, i, xSize, ySize).astype(np.float32)
        data3 = band3.ReadAsArray(j, i, xSize, ySize).astype(np.float32)

        mask = np.greater((data2

                           + data3), 0)
        ndvi = np.choose(mask, (-99, (data3 - data2) / (data2 + data3)))

        ndviBand.WriteArray(ndvi, j, i)

ndviBand.FlushCache()
ndviBand.SetNoDataValue(-99)
stats = ndviBand.GetStatistics(0, 1)

ndviDS.SetGeoTransform(ds.GetGeoTransform())
ndviDS.SetProjection(ds.GetProjection())

gdal.SetConfigOption('HFA_USE_RRD', 'YES')
ndviDS.BuildOverviews(overviewlist=[2, 4, 8, 16, 32, 64, 128])

del ds
del ndviDS

2. Assignment 5b

• Mosaic doq1.img and doq2.img together
• The pixel sizes are the same for both images
• Read in each image all at once – that will make it easier
• If you display it in ArcMap, change the symbology so it doesn’t stretch the data and it will look better
• Because it has different pyramid levels than the originals it might look offset when zoomed out, but zoom in and you’ll see no difference

# Assignment 5B: mosaic two images together
import os
import sys
from osgeo import gdal
from gdalconst import *

os.chdir('E:/data/GDAL/ospy_data5')

# register all of the GDAL drivers
gdal.AllRegister()

# read in doq1 and get info about it
ds1 = gdal.Open('doq1.img', GA_ReadOnly)
if ds1 is None:
    print('Cannot open doq1.img')
    sys.exit(1)
band1 = ds1.GetRasterBand(1)
rows1 = ds1.RasterYSize
cols1 = ds1.RasterXSize

# get the corner coordinates for doq1
transform1 = ds1.GetGeoTransform()
minX1 = transform1[0]
maxY1 = transform1[3]
pixelWidth1 = transform1[1]
pixelHeight1 = transform1[5]
maxX1 = minX1 + (cols1 * pixelWidth1)
minY1 = maxY1 + (rows1 * pixelHeight1)

# read in doq2 and get info about it
ds2 = gdal.Open('doq2.img', GA_ReadOnly)
if ds2 is None:
    print('Cannot open doq2.img')
    sys.exit(1)
band2 = ds2.GetRasterBand(1)
rows2 = ds2.RasterYSize
cols2 = ds2.RasterXSize

# get the corner coordinates for doq2
transform2 = ds2.GetGeoTransform()
minX2 = transform2[0]
maxY2 = transform2[3]
pixelWidth2 = transform2[1]
pixelHeight2 = transform2[5]
maxX2 = minX2 + (cols2 * pixelWidth2)
minY2 = maxY2 + (rows2 * pixelHeight2)

# get the corner coordinates for the output
minX = min(minX1, minX2)
maxX = max(maxX1, maxX2)
minY = min(minY1, minY2)
maxY = max(maxY1, maxY2)

# get the number of rows and columns for the output
cols = int((maxX - minX) / pixelWidth1)
rows = int((maxY - minY) / abs(pixelHeight1))

# compute the origin (upper left) offset for doq1
xOffset1 = int((minX1 - minX) / pixelWidth1)
yOffset1 = int((maxY1 - maxY) / pixelHeight1)

# compute the origin (upper left) offset for doq2
xOffset2 = int((minX2 - minX) / pixelWidth1)
yOffset2 = int((maxY2 - maxY) / pixelHeight1)

# create the output image
driver = ds1.GetDriver()
dsOut = driver.Create('mosiac.img', cols, rows, 1, band1.DataType)
bandOut = dsOut.GetRasterBand(1)

# read in doq1 and write it to the output
data1 = band1.ReadAsArray(0, 0, cols1, rows1)
bandOut.WriteArray(data1, xOffset1, yOffset1)

# read in doq2 and write it to the output
data2 = band2.ReadAsArray(0, 0, cols2, rows2)
bandOut.WriteArray(data2, xOffset2, yOffset2)

# compute statistics for the output
bandOut.FlushCache()
stats = bandOut.GetStatistics(0, 1)

# set the geotransform and projection on the output
geotransform = [minX, pixelWidth1, 0, maxY, 0, pixelHeight1]
dsOut.SetGeoTransform(geotransform)
dsOut.SetProjection(ds1.GetProjection())

# build pyramids for the output
gdal.SetConfigOption('HFA_USE_RRD', 'YES')
dsOut.BuildOverviews(overviewlist=[2, 4, 8, 16])

 

数据下载：https://download.youkuaiyun.com/download/qq_37935516/10797260