【GDAL學習】地圖代數與柵格資料的寫入
阿新 • • 發佈:2018-12-08
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