【GDAL学习】更多栅格数据处理函数——滑动窗口与过滤器

最新推荐文章于 2024-03-13 12:09:36 发布
最新推荐文章于 2024-03-13 12:09:36 发布
阅读量1.4k 收藏 3
点赞数
分类专栏： GDAL Python 文章标签： GDAL Python
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_37935516/article/details/84317032
版权
Python 同时被 2 个专栏收录
24 篇文章 1 订阅
订阅专栏
GDAL
10 篇文章 3 订阅
订阅专栏

例如设计一个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  seconds

2）参考答案给的代码，对输入图像的第一个通道进行高通滤波，输出图像只有一个通道

# 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 seconds

 2.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  seconds

2）参考答案给的代码，输出图形只有一个通道

# 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

 

众里寻你千百度
关注
专栏目录
python实现栅格滑动窗口统计值
Orange_GISer的博客
06-01 583
任务目标：在25km的栅格内统计90m栅格的均值。任务转化：构建一个滑动窗口统计均值。
遥感深度学习二维图像样本制作—异常值处理—零值（均值）填充—滑动窗口分割—数据增强（多波段遥感图像），python
最新发布
qq_63277082的博客
05-11 1073
最近在做多波段遥感语义分割的实验时，发现样本的制作是一个比较棘手的问题，因此想要整理下代码，以供学习交流。本文提供了python代码，可直接利用，具体步骤包括：异常值处理—零值（均值）填充—滑动窗口分割—数据增强（多波段遥感图像）
IDL语言实现滑动窗口
12-10
使用ENVIIDL语言实现9*9窗口的某一类值像元所占比例
[自用代码]ENVI输出的Sentinel-1图像预处理(翻转，滑动窗口裁剪)
weixin_38757163的博客
03-29 1502
import cv2 import os from osgeo import gdal import numpy as np def tif_jpg(rasterfile): in_ds = gdal.Open(rasterfile) # 打开样本文件 xsize = in_ds.RasterXSize # 获取行列数 ysize = in_ds.RasterYSize block_data = in_ds.ReadAsArray(0, 0, xsize, ysize
栅格影像高斯滤波平滑代码分享
weixin_43539146的博客
06-20 244
记得更改一下命名变量，同时这个是单波段不是多波段哦，多波段代码可以看参考链接里的小姐姐的分享。
【GDAL学习】地图代数与栅格数据的写入
qq_37935516的博客
11-20 732
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 statist...
在Python中用GDAL实现矢量对栅格的切割实例
09-17
在Python中，GDAL（Geospatial Data Abstraction Library）是一个强大的开源库，用于处理地理空间数据，包括矢量和栅格数据。本实例主要关注如何利用GDAL库来实现矢量对栅格数据的切割操作，即根据矢量边界（如形状...
基于GDAL实现多shp多影像的栅格数据批量裁剪
05-26
GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间数据转换库。...本资源使用Python的GDAL库，实现一个或多个矢量文件，批量对一个或多个栅格数据进行裁剪，大大减少了工作量。
通过GDAL\OGR获取栅格矢量边界Python代码
09-11
GDAL主要用于处理栅格数据，而OGR则专注于矢量数据操作。这两个库通常一起使用，因为它们都是GDAL项目的一部分，提供了广泛的数据读取、写入和处理功能。 标题提到的“通过GDAL\OGR获取栅格矢量边界”是一个常见的...
gdal处理栅格数据
08-25
利用python,安装GDAL库之后，来读取栅格数据，是一个简要的例子
gdal开源栅格空间数据转换库包括gdal-3.6.4源码，编译64位库文件
05-12
GDAL（Geospatial Data Abstraction Library）是一个强大的开源栅格空间数据转换库，它提供了对多种地理数据格式的支持，并且具有丰富的数据处理和转换功能。GDAL不仅仅是一个库，还包括一系列的命令行工具，使得...
gdal库读取tif影像坐标
pylyy的博客
12-29 1970
# coding=utf-8 from osgeo import gdal from gdalconst import * def readImage(img_path): data = [] # 以只读方式打开遥感影像 dataset = gdal.Open(img_path, GA_ReadOnly) if dataset is None: print("Unable to open image file.") return data
import gdal 报错
软考鸭
03-13 675
我的版本Python 3.7.9。（2）放到Scripts文件夹下。（1）查看python版本。（3）查看安装是否成功。进入python环境。
GDAL库初学
yezizp的专栏
06-04 1146
打开文件 在打开一个GDAL支持的栅格资料之前，必需要注册驱动。 每个驱动对应各自支持的格式。 通常这个会被GDALAllRegister（）函数完成。   C++ Python     import gdal     from gdalconst import *       dataset = gdal.Open( filename, GA_ReadOnly )
python reshape shape astype 理解
as472780551的博客
07-20 9421
x_train = x_train.reshape(x_train.shape[0], 1, img_rows, img_cols) x_test = x_test.reshape(x_test.shape[0], 1, img_rows, img_cols) from numpy import * l = zeros((5,4))#构建一个5*4的零矩阵 for i in range(5)...
C#调用GDAL矢量裁剪栅格
三千思的博客
01-05 814
在 GDAL 中，地理变换（transform）参数是一个长度为 6 的 double 数组，其中包含了栅格数据的投影坐标系下的坐标范围和像元大小。
基于OSGEO的遥感影像窗口滑动裁剪
John的博客
09-23 374
对遥感影像滑动窗口裁剪.
滑动窗口滤波的C语言实现（简单易移植）
qq_37662088的博客
06-01 1万+
滑动窗口滤波，含matlab代码与c语言实现。简单易懂易移植
Python地学分析 — GDAL进行遥感影像拼接
热门推荐
一叶知秋的博客
12-25 1万+
欢迎关注博主的微信公众号：“智能遥感”。 该公众号将为您奉上Python地学分析、爬虫、数据分析、Web开发、机器学习、深度学习等热门源代码。 本人的GitHub代码资料主页（持续更新中，多给Star，多Fork）： https://github.com/xbr2017 CSDN也在同步更新： https://blog.csdn.net/XBR_2014 “ 很多时候，我们需要研究...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值