1.根據座標讀取遙感影像的單個畫素值

# week 4: get pixel values at a set of coordinates by reading in one pixel at a time
import os, sys, time
from osgeo import gdal
from gdalconst import *

# start timing
startTime = time.time()

# coordinates to get pixel values for
xValues = [447520.0, 432524.0, 451503.0]
yValues = [4631976.0, 4608827.0, 4648114.0]

# set directory
os.chdir('E:/data/GDAL/ospy_data4')
# register all of the drivers
gdal.AllRegister()

# open the image
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:
    print('cannot open this file.')
    sys.exit(1)

# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# loop through the coordinates
for i in range(3):
    # get x,y
    x = xValues[i]
    y = yValues[i]

    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    # create a string to print out
    s = str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '
    # loop through the bands
    for j in range(bands):
        band = ds.GetRasterBand(j + 1)
        # read data and add the value to the string
        data = band.ReadAsArray(xOffset, yOffset, 1, 1)
        value = data[0, 0]
        s = s + str(value) + ' '
    # print out the data string
    print(s)
# figure out how long the script took to run
endTime = time.time()
print('The script took ' + str(endTime - startTime) + ' seconds.')
 

2.根據座標讀取遙感影像的單個畫素值，通過獲取影象所有畫素值獲得

# week 4:script to get pixel values at a set of coordinates by reading in entire bands
from osgeo import gdal
from gdalconst import *
import os, sys, time

# start timing
startTime = time.time()

# coordinates to get pixel values for
xValues = [447520.0, 432524.0, 451503.0]
yValues = [4631976.0, 4608827.0, 4648114.0]

# set directory
os.chdir('E:/data/GDAL/ospy_data4')

# register all of the drivers
gdal.AllRegister()

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

# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# create a list to store band data in
bandList = []

# read in bands and store all the data in bandList
for i in range(bands):
    band = ds.GetRasterBand(i + 1)
    data = band.ReadAsArray(0, 0, cols, rows)
    bandList.append(data)

# loop through the coordinates
for i in range(3):
    x = xValues[i]
    y = yValues[i]

    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    # create a string to print out
    s = str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '

    # loop through the bands and get the pixel value
    for j in range(bands):
        data = bandList[j]
        value = data[xOffset-1, yOffset-1]
        s = s + str(value)+' '

    # print out the data string
    print(s)

# figure out how long the script took to run
endTime = time.time()
print('The script took ' + str(endTime - startTime) + ' seconds.')
 

3.Assignment 4a：Read pixel values from an image

• Print out the pixel values for all three bands of aster.img at the points contained in sites.shp •
• Use any method of reading the raster data that you want, but I would suggest one pixel at a time (fastest in this case since we don't need much data)

1）自己寫的程式碼，現獲取所有sities的座標後再獲取全部的pixel values

# Assignment 4a: Read pixel values from an image
# write by myself
from osgeo import gdal
from gdalconst import *
import ogr, os, sys

# set the directory
os.chdir('E:/data/GDAL/ospy_data4')

# set the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# open the file
siteDS = driver.Open('sites.shp', 0)
if siteDS is None:
    print('Cannot open sites.shp')
    sys.exit(1)

# get the layer
siteLy = siteDS.GetLayer()
# get the feature
siteFe = siteLy.GetNextFeature()

# create lists to store the x,y coordination of the feature
sitePosX = []
sitePosY = []

# cycle the feature to get the x,y coordination
while siteFe:
    geom = siteFe.GetGeometryRef()
    fex = geom.GetX()
    fey = geom.GetY()
    sitePosX.append(fex)
    sitePosY.append(fey)
    siteFe.Destroy()
    siteFe = siteLy.GetNextFeature()

# create a txt file to store result
txt = open('output.txt', 'w')

# register all of the drivers
gdal.AllRegister()
# open the image
asterDS = gdal.Open('aster.img', GA_ReadOnly)
if asterDS is None:
    print('Cannot open aster.img')
    sys.exit(1)

# get the image size
rows = asterDS.RasterYSize
cols = asterDS.RasterXSize
bands = asterDS.RasterCount

# get the georeference info
transform = asterDS.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# get the pixel value of the coordination
for i in range(len(sitePosX)):
    x = sitePosX[i]
    y = sitePosY[i]
    # get the offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    s = str(i + 1) + ' ' + str(x) + ' ' + str(y) + ' ' + str(xOffset) + ' ' + str(yOffset) + ' '
    # get all bands of the pixel value
    for j in range(bands):
        band = asterDS.GetRasterBand(j + 1)
        data = band.ReadAsArray(xOffset, yOffset, 1, 1)
        value = data[0, 0]
        s = s + str(value) + ' '
    s = s + '\n'
    # write to a text file
    txt.write(s)
    print(s)

# destroy datasource
siteDS.Destroy()

2）參考答案給的程式碼，每獲取一個site座標就獲取該座標的pixel value

import os, sys
try:
    from osgeo import ogr, gdal
    from osgeo.gdalconst import *

    os.chdir('E:/data/GDAL/ospy_data4')
except ImportError:
    import ogr, gdal
    from gdalconst import *

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

# open the shapefile and get the layer
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open('sites.shp')
if shp is None:
    print('Could not open sites.shp')
    sys.exit(1)
shpLayer = shp.GetLayer()

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

# open the image
img = gdal.Open('aster.img', GA_ReadOnly)
if img is None:
    print('Could not open aster.img')
    sys.exit(1)

# get image size
rows = img.RasterYSize
cols = img.RasterXSize
bands = img.RasterCount

# get georeference info
transform = img.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# loop through the features in the shapefile
feat = shpLayer.GetNextFeature()
while feat:
    geom = feat.GetGeometryRef()
    x = geom.GetX()
    y = geom.GetY()

    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)

    # create a string to print out
    s = feat.GetFieldAsString('ID') + ' '

    # loop through the bands
    for j in range(bands):
        band = img.GetRasterBand(j + 1)  # 1-based index

        # read data and add the value to the string
        data = band.ReadAsArray(xOffset, yOffset, 1, 1)
        value = data[0, 0]
        s = s + str(value) + ' '

    # print out the data string
    print(s)

    # get the next feature
    feat.Destroy()
    feat = shpLayer.GetNextFeature()

# close the shapefile
shp.Destroy()

4.為了更有效率地讀取柵格，一種更好的方法是分塊讀取資料。資料檔案中給出了一個utils.py檔案中的GetBlockSize()函式可以實現獲取影象中的xBlockSize和yBlockSize，但是在我的電腦裡不能執行。於是嘗試手動定義xBlockSize和yBlockSize，更改了幾次數值發現程式執行結果都是一樣的，因此我認為手動賦值給xBlockSize和yBlockSize也是可行的。下面程式碼是迴圈分塊讀取柵格資料的樣例，具體原理可參見PPT

from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys
import utils

os.chdir('E:/data/GDAL/ospy_data4')
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

band = ds.GetRasterBand(1)
# 由於utils.py檔案中找不到_gdal模組，因此以下三行程式碼不能使用
# xblocksize和yblocksize改為手動設定
# 更改了幾個值發現程式執行最終結果不變
# blockSizes = utils.GetBlockSize(band)
# xBlockSize = blockSizes[0]
# yBlockSize = blockSizes[1]

xBlockSize = 20
yBlockSize = 20

count = 0

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

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

        data = band.ReadAsArray(j, i, numCols, numRows).astype(np.float)
        mask = np.greater(data, 0)
        count = count + np.sum(np.sum(mask))

print('Number of non-zero pixels: ', count)

5. Assignment 4b：

• Write a script to calculate the average pixel value for the first band in aster.img
• Read in the data one block at a time 
• Do the calculation two ways ：
  • Including zeros in the calculation
  • Ignoring zeros in the calculation

1）自己寫的程式碼

# Assignment 4b:Write a script to calculate the average pixel value for the first band in aster.img
# write by myself
# import models
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys

# set the directory
os.chdir('E:/data/GDAL/ospy_data4')

# register all drivers
gdal.AllRegister()

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

# get the raster size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# get the first band and set blocksize
band = ds.GetRasterBand(1)
xBlockSize = 10
yBlockSize = 10

sumIncludeZero = 0.0
zeroCount = 0

# cycle row blocksize
for i in range(0, rows, yBlockSize):
    # set yblocksize
    if (i + yBlockSize) < rows:
        ySize = yBlockSize
    else:
        ySize = rows - i

    # cycle col blocksize
    for j in range(0, cols, xBlockSize):
        # set xblocksize
        if (j + xBlockSize) < cols:
            xSize = xBlockSize
        else:
            xSize = cols - j

        # read block data
        data = band.ReadAsArray(j, i, xSize, ySize).astype(np.float)
        sumIncludeZero = sumIncludeZero + np.sum(np.sum(data))
        mask = np.greater(data, 0)
        zeroCount = zeroCount + np.sum(np.sum(mask))

# calculate
avg_includeZero = sumIncludeZero / (rows * cols)
avg_ignoreZero = sumIncludeZero / zeroCount
print('Average include zero: ', avg_includeZero)
print('Average ignore zero: ', avg_ignoreZero)

2）參考答案給的程式碼

# Assignment 4b:Write a script to calculate the average pixel value for the first band in aster.img
#reference answer
import os, sys
# import utils

use_numeric = True

try:
    from osgeo import ogr, gdal
    from osgeo.gdalconst import *
    import numpy

    os.chdir('E:/data/GDAL/ospy_data4')
    use_numeric = False
except ImportError:
    import ogr, gdal
    from gdalconst import *
    import Numeric

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

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

# open the image
ds = gdal.Open('aster.img', GA_ReadOnly)
if ds is None:
    print('Could not open aster.img')
    sys.exit(1)

# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# # get the band and block sizes
# band = ds.GetRasterBand(1)
# blockSizes = utils.GetBlockSize(band)
# xBlockSize = blockSizes[0]
# yBlockSize = blockSizes[1]
band = ds.GetRasterBand(1)
xBlockSize = 20
yBlockSize = 20

# initialize variables
count = 0
total = 0

# loop through the rows
for i in range(0, rows, yBlockSize):
    if i + yBlockSize < rows:
        numRows = yBlockSize
    else:
        numRows = rows - i

    # loop through the columns
    for j in range(0, cols, xBlockSize):
        if j + xBlockSize < cols:
            numCols = xBlockSize
        else:
            numCols = cols - j

        # read the data and do the calculations
        if use_numeric:
            data = band.ReadAsArray(j, i, numCols, numRows).astype(Numeric.Float)
            mask = Numeric.greater(data, 0)
            count = count + Numeric.sum(Numeric.sum(mask))
            total = total + Numeric.sum(Numeric.sum(data))
        else:
            data = band.ReadAsArray(j, i, numCols, numRows).astype(numpy.float)
            mask = numpy.greater(data, 0)
            count = count + numpy.sum(numpy.sum(mask))
            total = total + numpy.sum(numpy.sum(data))

# print results
print('Ignoring 0:', total / count)
print('Including 0:', total / (rows * cols))

資料下載：https://download.csdn.net/download/qq_37935516/10795361