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.csdn.net/download/qq_37935516/10797260