Python地理数据处理 十三:栅格数据处理(一)
目录
- 1.数据重采样
- 2. 字节序列
- 3. 子数据集
1.数据重采样
重采样是指根据一类象元的信息内插出另一类象元信息的过程。在遥感中,重采样是从高分辨率遥感影像中提取出低分辨率影像的过程。常用的重采样方法有最邻近内插法、双线性内插法和三次卷积法内插。
ReadAsArray函数可以重采样读取的数据,并且指定输出缓冲区大小或传递一个已有的缓冲区数组。
函数格式:
band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize], [buf_ysize], [buf_obj])
- xoff是开始阅读的专栏,默认值为0。
- yoff是开始阅读的行,默认值为0。
- win_xsize是要读取的列数,默认为全部读取。
- win_ysize是要读取的行数,默认为全部读取。
- buf_xsize是输出数组中的列数,默认值为使用win_xsize的值。
- buf_ysize是输出数组中的行数,默认值为使用win_ysize的值。
- buf_obj是一个NumPy数组,用于将数据放入其中,而不是创建一个新数组。如果需要,数据将被重新采样以适合此数组,对应的值也将转换为该数组的数据类型。
具有较大尺寸的数组将重采样为较小的像素,更小的尺寸的数组,将采用最近邻插值法重采样为更大的像素
重采样为更小的像素:
需要提供一个比原有数据维数更大的数组,使得重复使用像素值来填充目标数组:
import os
from osgeo import gdalos.chdir(r'D:\geodata\Landsat\Washington')in_ds = gdal.Open('p047r027_7t20000730_z10_nn10.tif')
in_band = in_ds.GetRasterBand(1)# 计算输出行列数
# 输入数翻倍,因为我将像素大小减半
out_rows = in_band.YSize * 2
out_columns = in_band.XSize * 2# 创建输出数据集
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('band1_resampled.tif',out_columns, out_rows)# 编辑地理变换
# 像素变为原来的 1/4
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] /= 2
geotransform[5] /= 2
out_ds.SetGeoTransform(geotransform)# 读取数据时,指定一个较大的缓冲
data = in_band.ReadAsArray(buf_xsize=out_columns, buf_ysize=out_rows)# 将数据写入输出光栅
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)# 为较大的图像构建合适的概视图
out_band.FlushCache()
out_band.ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16, 32, 64])del out_ds
图像展示:
属性信息:
2. 字节序列
1.ReadRaster([xoff], [yoff], [xsize], [ysize], [buf_xsize], [buf_ysize], [buf_type], [band_list],[buf_pixel_space], [buf_line_space], [buf_band_space])
- xoff是列读取起点,默认值为0。
- yoff是行读取起点,默认值为0。
- xsize是读取的列数,默认为全部读取。
- ysize是读取的行数,默认为全部读取。
- buf_xsize是输出数组里的列数,默认值为使用 win_xsize 值。如果此值不同于 win_xsize,则数据将重采样。
- buf_ysize是输出数组里的行数,默认值为使用 win_ysize 值。如果此值不同于 win_ysize,则数据将重采样。
- buf_type 是返回序列的GDAL目标数据类型,默认值与源数据相同。
- band_list是要读取的波段列表,默认读取所有波段。
- buf_pixel_space是序列中像素之间的字节偏移,默认值为buf_type的大小。
- buf_line_space是序列中行间的字节偏移,默认值是 buf_type 乘以 xsize 。
- buf_band_space是序列中列间的字节偏移,默认值是 buf_line_space 乘以 ysize 。
import os
import numpy as np
from osgeo import gdaldata_dir = r'D:\geodata'os.chdir(os.path.join(data_dir, 'Landsat', 'Washington'))
ds = gdal.Open('Landsat_color.tif')
data = ds.ReadRaster(1400, 6000, 2, 2, band_list=[1])
print(data)# 取出第一个值,将从字节转换为数字
print(data[0])# 尝试更改第一个像素的值。
# 输出失败,因为不能更改字节字符串
data[0] = 50# 将字节字符串转换为字节数组
# 更改第一个值,输出成功
bytearray_data = bytearray(data) # bytearray是字节数组
bytearray_data[0] = 50
print(bytearray_data[0])
b'\x1c\x1d\x1c\x1e'
28
50
2. 将字节字符串转换为像素值的元组:
import struct
tuple_data = struct.unpack('B' * 4, data) # 指定4个字节
print(tuple_data)
(28, 29, 28, 30)
3. 将元组转换为numpy数组:
numpy_data1 = np.array(tuple_data)
print(numpy_data1)
[28 29 28 30]
4. 将字节字符串转换为numpy数组:
# 将字节字符串转换为numpy数组
# 重构一个numpy数组,使其具有2行2列,就像读入的原始数据一样
numpy_data2 = np.fromstring(data, np.int8)
reshaped_data = np.reshape(numpy_data2, (2,2))
print(reshaped_data)
[[28 29][28 30]]
5. 利用字节序列将图像重采样为更大的像素尺寸:
import os
import numpy as np
from osgeo import gdalos.chdir(r'D:\geodata\Landsat\Washington')# 打开输入栅格
in_ds = gdal.Open('Landsat_color.tif')# 计算输出的行和列的数量
# 输入数字的一半,因为要使像素的两倍大
out_rows = int(in_ds.RasterYSize / 2)
out_columns = int(in_ds.RasterXSize / 2)
num_bands = in_ds.RasterCount# 创建输出数据集
gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create('Landsat_color_resampled.tif',out_columns, out_rows, num_bands)# 编辑地理变换,让像素尺寸变大
# 像素宽度、像素高度加倍
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())
geotransform[1] *= 2
geotransform[5] *= 2
out_ds.SetGeoTransform(geotransform)# 利用较小的缓冲来读写数据
data = in_ds.ReadRaster(buf_xsize=out_columns, buf_ysize=out_rows)# 将数据写入栅格中
out_ds.WriteRaster(0, 0, out_columns, out_rows, data)# 为较小的图像创建合适数量的概视图
out_ds.FlushCache()
for i in range(num_bands):out_ds.GetRasterBand(i + 1).ComputeStatistics(False)
out_ds.BuildOverviews('average', [2, 4, 8, 16])del out_ds
结果显示:
3. 子数据集
MODIS图像(中分辨率成像光谱仪)是一个分层数据格式(HDF)文件。如果数据包含子数据集,可用 GetSubDatasets() 函数得到他们的子数据集列表,再利用这个列表信息打开所需要的子数据集。
数据集结构:
import os
import numpy as np
from osgeo import gdaldata_dir = r'D:\geodata'# 从MODIS文件中获取子数据集的名称和描述
# 并打印,输出NDVI(归一化植被指数)
os.chdir(os.path.join(data_dir, 'Modis'))
ds = gdal.Open('MYD13Q1.A2014313.h20v11.005.2014330092746.hdf')
subdatasets = ds.GetSubDatasets() # GetSubDatasets函数返回一个元组列表
print('Number of subdatasets: {}'.format(len(subdatasets)))
for sd in subdatasets:print('Name: {0}\nDescription:{1}\n'.format(*sd))# 打开Modis文件中的第一个子数据集:[0][0]
# 第五个数据集为: [4][0]
ndvi_ds = gdal.Open(subdatasets[0][0])# 通过打印尺寸来确保正常工作
print('Dataset dimensions: {} {}'.format(ndvi_ds.RasterXSize, ndvi_ds.RasterYSize))# 读取数据之前先获取波段
ndvi_band = ndvi_ds.GetRasterBand(1)
print('Band dimensions: {} {}'.format(ndvi_band.XSize, ndvi_band.YSize))
结果显示:
Python地理数据处理 十三:栅格数据处理(一)相关推荐
- Python读写矢量数据(2)矢量数据写入(属性数据)——Python地理数据处理学习分享
这一节主要介绍矢量数据的写入(只有属性数据,无几何),如果有读者没有读取的基础建议先看一下上一篇文章,需要对矢量数据读取有一定的了解才能继续学习本节.在这里我们用到的数据仍为goble文件夹下的数据, ...
- Python Gdal 栅格数据处理之hgt转tif数据
Python Gdal 栅格数据处理之.hgt转.tif数据 1 介绍 2 Python代码 3 结果展示 1 介绍 SRTM是地形数据,其存储形式为.hgt格式,根据精度的不同可以分为两种:SR ...
- 《Python 地理数据处理》by Chris Garrard
科研小白从头开始学习用Python处理栅格数据.矢量数据等. 下面就跟我一起开启<Python 地理数据处理>的学习之旅吧!!! 本书主要利用Python + GDAL 等相关库进行地理空 ...
- inner join 重复数据_Ramp;Python Data Science 系列:数据处理(2)
承接R&Python Data Science 系列:数据处理(1)继续介绍剩余的函数. 1 衍生字段函数 主要有两个函数,mutate()和transmute(),两个函数在Python和R ...
- 【优秀课设】武汉光迅科技22校招笔试Python题改进(增加GUI)——基于Python的125温度传感器模块数据处理
武汉光迅科技22校招笔试Python题改进(增加GUI) 基于Python的125温度传感器模块数据处理 原本的基础代码: blog.csdn.net/weixin_53403301/article/ ...
- python txt文件处理软件,对python .txt文件读取及数据处理方法总结
1.处理包含数据的文件 最近利用python读取txt文件时遇到了一个小问题,就是在计算两个np.narray()类型的数组时,出现了以下错误: 作为一个python新手,遇到这个问题后花费了挺多时间 ...
- 【笔记】《Python地理空间分析指南(第2版)》
转载地址:https://blog.csdn.net/jianbinzheng/article/details/80215228 概述部分 地理空间数据 地理空间技术概览 Python地理空间分析工具 ...
- python空间分析_读书笔记——《python地理空间分析指南》
本文为<Python地理空间分析指南(第2版)>的读书摘录,顺便挖个坑,进一步对python的几个包做学习整理. 本笔记的用途:了解python地理空间处理的技术框架和实现途径. 第三章 ...
- python学习[第十三篇] 条件和循环
python学习[第十三篇] 条件和循环 if语句 单一if 语句 if语句有三个部分构成,关键字if本身,判断结果真假的条件表达式,以及表达式为真或非0是执行的代码 if expression: e ...
- python地理空间分析指南pdf邓世超_Python地理空间分析指南(第2版)源代码.zip
[实例简介] Python地理空间分析指南(第2版)的随书源代码,需要的朋友可以下载一下~~ [实例截图] [核心代码] Python地理空间分析指南(第2版)源代码 └── Python地理空间分析 ...
最新文章
- 关于Java为什么配置好环境变量但是不能在命令行cmd运行javac的问题
- 微软陆续更新Win8应用 否认靠金钱争取开发者
- awgn信道中的噪声功率谱密度_从OFC2020看高级算法在光通信中的应用
- 2019年最流行的10个前端框架
- 深入浅出:Microsoft分布式事务处理协调器
- selenium-python中文文档
- goodFeaturesToTrack——Shi-Tomasi角点检测
- 电脑怎么隐藏文件夹?6个步骤完成!
- 好看的皮囊 · 也是大自然的杰作 · 全球高质量 · 美图 · 集中营 · 美女 · 2017-08-20期...
- 8cm等于多少像素_一寸照片像素是多少
- UEBA架构设计之路1
- 转-基于OpenGL的3D天空仿真
- 一本通 1335:【例2-4】连通块
- 计算机毕设之餐厅选座订餐系统的设计与实践
- “SEO是什么意思?”Kyw的通俗回答
- 查询linux系统中文件名颜色分别代表什么
- 微信怎样查绑定的服务器地址,你的微信绑定了哪些网站和应用?这个方法可以一键查看......
- 执行力在ERP系统中发挥的作用
- XidianOJ 1081 锘爷与矩阵
- spring整合quartz框架定时任务实战
热门文章
- 【R】 Error in is.data.frame(x) : (list) object cannot be coerced to type 'double'
- php+mysql+jquery瀑布流
- JavaIO流——文件的读取与传输
- 走进龙芯3A3000(二)安装Gentoo N64
- wincc变量数据归档(案例)
- 制造主数据集成开发心得
- AKULAKU笔试题(还有1题未答)
- 7.1 认识Access报表
- STM32F401CCU6 核心板的功能描述(针对采集数据)
- centos有道linux安装,centos7安装有道词典(不能发音和取词)