【GDAL學習】更多柵格資料處理函式——滑動視窗與過濾器
阿新 • • 發佈：2018-12-08
例如設計一個3 x 3的滑動視窗，寫演算法執行就有兩種方式：
1.pixel by piexl每個進行逐畫素運算，效率太低，速度慢
2.使用 slice切片形式迴圈，效率高，速度快
兩個作業就是分別用pixel和slice方式完成高通濾波操作進行對比
1.Assignment 6a
- Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img
- The output data type will be Float
- Use pixel notation (that’s why you’re doing it on smallaster.img instead of aster.img)
1）自己寫的程式碼，對輸入影象的三個通道進行高通濾波，輸出影象有三個影象
# Assignment 6a:Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img # writen by myself from osgeo import gdal from gdalconst import * import numpy as np import os import sys import time # record start time startTime = time.time() # set the directory os.chdir('E:/data/GDAL/ospy_data6') # register all driver gdal.AllRegister() # open file ds = gdal.Open('smallaster.img', GA_ReadOnly) if ds is None: print('Cannot open smallaster.img') sys.exit(1) # get the size of smallaster.img rows = ds.RasterYSize cols = ds.RasterXSize bands = ds.RasterCount # create a data array to store result outdata = np.zeros((bands, rows, cols)) # filt = np.full((3, 3), 0.111).astype(np.float32)# low pass filt = np.array([[-0.7, -1.0, -0.7], [-1.0, 6.8, -1.0], [-0.7, -1.0, -0.7]]).astype(np.float32) # high pass # get current driver driver = ds.GetDriver() # create output file outds = driver.Create('highpass_smallaster.img', cols, rows, bands, GDT_Float32) if outds is None: print('Cannot open highpass_smallaster.img') sys.exit(1) # cycle bands to calculate for i in range(bands): # get each band band = ds.GetRasterBand(i + 1) # get each band data data = band.ReadAsArray(0, 0, cols, rows).astype(np.float32) # cycle row direction for j in range(1, rows - 1): # cycle col direction for k in range(1, cols - 1): # high pass filter calculate outdata[i, j, k] = ( data[j - 1, k - 1] * filt[0, 0] + data[j - 1, k] * filt[0, 1] + data[j - 1, k + 1] * filt[0, 2] + data[j, k - 1] * filt[1, 0] + data[j, k] * filt[1, 1] + data[j, k + 1] * filt[1, 2] + data[j + 1, k - 1] * filt[2, 0] + data[j + 1, k] * filt[2, 1] + data[j + 1, k + 1] * filt[2, 2]) # get current band to store result outband = outds.GetRasterBand(i + 1) outband.WriteArray(outdata[i, :, :], 0, 0) # write to disk outband.FlushCache() stats = outband.GetStatistics(0, 1) # set the Geotransform and projection outds.SetGeoTransform(ds.GetGeoTransform()) outds.SetProjection(ds.GetProjection()) del ds del outds # record the end time endTime = time.time() # print the time for run last print('Last ', endTime - startTime, ' seconds')
程式執行時間：
Last 104.1204240322113 seconds2）參考答案給的程式碼，對輸入影象的第一個通道進行高通濾波，輸出影象只有一個通道
# homework 6a:high-pass filter using pixel notation # reference answer import os, sys, numpy, gdal, time from gdalconst import * start = time.time() os.chdir('E:/data/GDAL/ospy_data6') # register all of the GDAL drivers gdal.AllRegister() # open the image inDs = gdal.Open('smallaster.img', GA_ReadOnly) if inDs is None: print('Could not open smallaster.img') sys.exit(1) # get image size rows = inDs.RasterYSize cols = inDs.RasterXSize # read the input data inBand = inDs.GetRasterBand(1) inData = inBand.ReadAsArray(0, 0, cols, rows).astype(numpy.float) # do the calculation outData = numpy.zeros((rows, cols), numpy.float) for i in range(1, rows - 1): for j in range(1, cols - 1): outData[i, j] = ((-0.7 * inData[i - 1, j - 1]) + (-1.0 * inData[i - 1, j]) + (-0.7 * inData[i - 1, j + 1]) + (-1.0 * inData[i, j - 1]) + (6.8 * inData[i, j]) + (-1.0 * inData[i, j + 1]) + (-0.7 * inData[i + 1, j - 1]) + (-1.0 * inData[i + 1, j]) + (-0.7 * inData[i + 1, j + 1])) # create the output image driver = inDs.GetDriver() outDs = driver.Create('highpass1.img', cols, rows, 1, GDT_Float32) if outDs is None: print('Could not create highpass1.img') sys.exit(1) outBand = outDs.GetRasterBand(1) # write the output data outBand.WriteArray(outData, 0, 0) # flush data to disk, set the NoData value and calculate stats outBand.FlushCache() stats = outBand.GetStatistics(0, 1) # georeference the image and set the projection outDs.SetGeoTransform(inDs.GetGeoTransform()) outDs.SetProjection(inDs.GetProjection()) inDs = None outDs = None print(time.time() - start, 'seconds')
程式執行時間：
34.041293144226074 seconds2.Assignment 6b
- Use a 3x3 high pass filter to detect edges in band 1 of aster.img (good idea to test on smallaster.img first)
- The output data type will be Float
- Use slice notation
1）自己寫的程式碼，輸出影象有三個通道
# Assignment 6b:Use a 3x3 high pass filter to detect edges in band 1 of smallaster.img
# writen by myself
from osgeo import gdal
from gdalconst import *
import numpy as np
import os
import sys
import time
# record start time
startTime = time.time()
# set the directory
os.chdir('E:/data/GDAL/ospy_data6')
# register all driver
gdal.AllRegister()
# open file
ds = gdal.Open('smallaster.img', GA_ReadOnly)
if ds is None:
print('Cannot open smallaster.img')
sys.exit(1)
# get the size of smallaster.img
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount
# create a data array to store result
outdata = np.zeros((bands, rows, cols))
# filt = np.full((3, 3), 0.111).astype(np.float32)# low pass
filt = np.array([[-0.7, -1.0, -0.7], [-1.0, 6.8, -1.0], [-0.7, -1.0, -0.7]]).astype(np.float32) # high pass
# get current driver
driver = ds.GetDriver()
# create output file
outds = driver.Create('highpass_smallaster_slice.img', cols, rows, bands, GDT_Float32)
if outds is None:
print('Cannot open highpass_smallaster_slice.img')
sys.exit(1)
# cycle bands to calculate
for i in range(bands):
# get each band
band = ds.GetRasterBand(i + 1)
# get each band data
data = band.ReadAsArray(0, 0, cols, rows).astype(np.float32)
outdata[i, 1:rows - 1, 1:cols - 1] = (
data[0:rows - 2, 0:cols - 2] * filt[0, 0] + data[0:rows - 2, 1:cols - 1] * filt[0, 1] +
data[0:rows - 2, 2:cols] * filt[0, 2] + data[1:rows - 1, 0:cols - 2] * filt[1, 0] +
data[1:rows - 1, 1:cols - 1] * filt[1, 1] +
data[1:rows - 1, 2:cols] * filt[1, 2] + data[2:rows, 0:cols - 2] * filt[2, 0] +
data[2:rows, 1:cols - 1] * filt[2, 1] + data[2:rows, 2:cols] * filt[2, 2])
# get current band to store result
outband = outds.GetRasterBand(i + 1)
outband.WriteArray(outdata[i, :, :], 0, 0)
# write to disk
outband.FlushCache()
stats = outband.GetStatistics(0, 1)
# set the Geotransform and projection
outds.SetGeoTransform(ds.GetGeoTransform())
outds.SetProjection(ds.GetProjection())
del ds
del outds
# record the end time
endTime = time.time()
# print the time for run last
print('Last ', endTime - startTime, ' seconds')
# last 3s
程式執行時間，可以看到速度由104s提升到了3s
Last 3.429971933364868 seconds2）參考答案給的程式碼，輸出圖形只有一個通道
# homework 6b：high-pass filter using slice notation
# reference answer
import os, sys, numpy, gdal, time
from gdalconst import *
start = time.time()
os.chdir('E:/data/GDAL/ospy_data6')
# register all of the GDAL drivers
gdal.AllRegister()
# open the image
inDs = gdal.Open('aster.img', GA_ReadOnly)
if inDs is None:
print('Could not open aster.img')
sys.exit(1)
# get image size
rows = inDs.RasterYSize
cols = inDs.RasterXSize
# read the input data
inBand = inDs.GetRasterBand(1)
inData = inBand.ReadAsArray(0, 0, cols, rows).astype(numpy.float)
# do the calculation
outData = numpy.zeros((rows, cols), numpy.float)
outData[1:rows - 1, 1:cols - 1] = (
(-0.7 * inData[0:rows - 2, 0:cols - 2]) +
(-1.0 * inData[0:rows - 2, 1:cols - 1]) + (-0.7 * inData[0:rows - 2, 2:cols]) +
(-1.0 * inData[1:rows - 1, 0:cols - 2]) + (6.8 * inData[1:rows - 1, 1:cols - 1]) +
(-1.0 * inData[1:rows - 1, 2:cols]) + (-0.7 * inData[2:rows, 0:cols - 2]) +
(-1.0 * inData[2:rows, 1:cols - 1]) + (-0.7 * inData[2:rows, 2:cols]))
# create the output image
driver = inDs.GetDriver()
outDs = driver.Create('highpass3.img', cols, rows, 1, GDT_Float32)
if outDs is None:
print('Could not create highpass3.img')
sys.exit(1)
outBand = outDs.GetRasterBand(1)
# write the output data
outBand.WriteArray(outData, 0, 0)
# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
stats = outBand.GetStatistics(0, 1)
# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())
inDs = None
outDs = None
print(time.time() - start, 'seconds')
資料下載：https://download.csdn.net/download/qq_37935516/10798222