python 黑白tif提取边界像素坐标_python+gdal实现MODIS_HDF文件转TIF（可批量）

最新推荐文章于 2022-04-11 22:16:39 发布
最新推荐文章于 2022-04-11 22:16:39 发布
阅读量421 收藏
点赞数
文章标签： python 黑白tif提取边界像素坐标
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39751871/article/details/111635437
版权

def06e58546fa13eae5d7f6193496197.png

准备

安装python以及gdal

馨意：基于python的遥感图像处理(1.1)--Anaconda安装步骤​zhuanlan.zhihu.com
c1eee02ff78128aa9e5cb50591ad4aad.png

正文

本文以提取MODIS植被指数产品MOD13A3的NDVI数据为例：

首先打开hdf数据集并查看子数据集和元数据：

#  gdal打开hdf数据集
datasets = gdal.Open(r"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf")

#  获取hdf中的子数据集
SubDatasets = datasets.GetSubDatasets()
#  获取子数据集的个数
SubDatasetsNum =  len(datasets.GetSubDatasets())
#  输出各子数据集的信息
print("子数据集一共有{0}个: ".format(SubDatasetsNum))
for i in range(SubDatasetsNum):
    print(datasets.GetSubDatasets()[i])

#  获取hdf中的元数据
Metadata = datasets.GetMetadata()
#  获取元数据的个数
MetadataNum = len(Metadata)
#  输出各子数据集的信息
print("元数据一共有{0}个: ".format(MetadataNum))
for key,value in Metadata.items():
    print('{key}:{value}'.format(key = key, value = value))

输出为：

子数据集一共有11个: 
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly NDVI"', '[1200x1200] 1 km monthly NDVI MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly EVI"', '[1200x1200] 1 km monthly EVI MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly VI Quality"', '[1200x1200] 1 km monthly VI Quality MOD_Grid_monthly_1km_VI (16-bit unsigned integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly red reflectance"', '[1200x1200] 1 km monthly red reflectance MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly NIR reflectance"', '[1200x1200] 1 km monthly NIR reflectance MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly blue reflectance"', '[1200x1200] 1 km monthly blue reflectance MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly MIR reflectance"', '[1200x1200] 1 km monthly MIR reflectance MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly view zenith angle"', '[1200x1200] 1 km monthly view zenith angle MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly sun zenith angle"', '[1200x1200] 1 km monthly sun zenith angle MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly relative azimuth angle"', '[1200x1200] 1 km monthly relative azimuth angle MOD_Grid_monthly_1km_VI (16-bit integer)')
('HDF4_EOS:EOS_GRID:"E:Remote_Sensing_DataTVDIMOD13A32012MOD13A3.A2012001.h26v04.005.hdf":MOD_Grid_monthly_1km_VI:"1 km monthly pixel raliability"', '[1200x1200] 1 km monthly pixel raliability MOD_Grid_monthly_1km_VI (8-bit integer)')

元数据一共有163个: 
ALGORITHMPACKAGEACCEPTANCEDATE:102004
ALGORITHMPACKAGEMATURITYCODE:Normal
......元数据太多，只展示这几个为例
DATACOLUMNS:1200
DATAROWS:1200
......
GRINGPOINTLATITUDE.1:39.8624390029821, 50.1154198026055, 49.9369825110409, 39.7325278324118
GRINGPOINTLONGITUDE.1:104.235558676363, 124.607874257256, 140.825259076411, 117.770602913265
......
VERTICALTILENUMBER:04
WESTBOUNDINGCOORDINATE:104.43258313171

我们可以发现子数据集中的第一个数据为NDVI数据，我们将数值提取出来：

#  第一个子数据集合,也就是NDVI数据
DatasetNDVI = datasets.GetSubDatasets()[0][0]
RasterNDVI = gdal.Open(DatasetNDVI)
NDVI = RasterNDVI.ReadAsArray()

如果我们转成tif格式数据，需要知道图像的仿射变换矩阵，仿射变换矩阵可以通过图像的四个角的经纬度计算得出(下面链接可以帮助了解一下仿射变换矩阵)：

馨意：python+GDAL遥感图像像素坐标与地理坐标转换（仿射变换）​zhuanlan.zhihu.com 
#  获取四个角的维度
Latitudes = Metadata["GRINGPOINTLATITUDE.1"]
#  采用", "进行分割
LatitudesList = Latitudes.split(", ")
#  获取四个角的经度
Longitude = Metadata["GRINGPOINTLONGITUDE.1"]
#  采用", "进行分割
LongitudeList = Longitude.split(", ")

# 图像四个角的地理坐标
GeoCoordinates = np.zeros((4, 2), dtype = "float32")
GeoCoordinates[0] = np.array([float(LongitudeList[0]),float(LatitudesList[0])])
GeoCoordinates[1] = np.array([float(LongitudeList[1]),float(LatitudesList[1])])
GeoCoordinates[2] = np.array([float(LongitudeList[2]),float(LatitudesList[2])])
GeoCoordinates[3] = np.array([float(LongitudeList[3]),float(LatitudesList[3])])

#  列数
Columns = float(Metadata["DATACOLUMNS"])
#  行数
Rows = float(Metadata["DATAROWS"])
#  图像四个角的图像坐标
PixelCoordinates = np.array([[0, 0],
                [Columns - 1, 0],
                [Columns - 1, Rows - 1],
                [0, Rows - 1]], dtype = "float32")
    
#  计算仿射变换矩阵
from scipy.optimize import leastsq
def func(i):
    Transform0, Transform1, Transform2, Transform3, Transform4, Transform5 = i[0], i[1], i[2], i[3], i[4], i[5]
    return [Transform0 + PixelCoordinates[0][0] * Transform1 + PixelCoordinates[0][1] * Transform2 - GeoCoordinates[0][0],
            Transform3 + PixelCoordinates[0][0] * Transform4 + PixelCoordinates[0][1] * Transform5 - GeoCoordinates[0][1],
            Transform0 + PixelCoordinates[1][0] * Transform1 + PixelCoordinates[1][1] * Transform2 - GeoCoordinates[1][0],
            Transform3 + PixelCoordinates[1][0] * Transform4 + PixelCoordinates[1][1] * Transform5 - GeoCoordinates[1][1],
            Transform0 + PixelCoordinates[2][0] * Transform1 + PixelCoordinates[2][1] * Transform2 - GeoCoordinates[2][0],
            Transform3 + PixelCoordinates[2][0] * Transform4 + PixelCoordinates[2][1] * Transform5 - GeoCoordinates[2][1],
            Transform0 + PixelCoordinates[3][0] * Transform1 + PixelCoordinates[3][1] * Transform2 - GeoCoordinates[3][0],
            Transform3 + PixelCoordinates[3][0] * Transform4 + PixelCoordinates[3][1] * Transform5 - GeoCoordinates[3][1]]
#  最小二乘法求解
GeoTransform = leastsq(func,np.asarray((1,1,1,1,1,1)))
print(GeoTransform)

为了方便处理数据，我们获取数据的时间并以此为文件名保存：

#  获取数据时间
date = Metadata["RANGEBEGINNINGDATE"]

#  保存为tif
def array2raster(TifName, GeoTransform, array):
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(TifName, cols, rows, 1, gdal.GDT_Float32)
    # 括号中两个0表示起始像元的行列号从(0,0)开始
    outRaster.SetGeoTransform(tuple(GeoTransform))
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

TifName = date + ".tif"
array2raster(TifName, GeoTransform[0], NDVI)

至此，大功告成。

批量转换代码：

import gdal, osr
import numpy as np
import os

#  数组保存为tif
def array2raster(TifName, GeoTransform, array):
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(TifName, cols, rows, 1, gdal.GDT_Float32)
    # 括号中两个0表示起始像元的行列号从(0,0)开始
    outRaster.SetGeoTransform(tuple(GeoTransform))
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
    
#  hdf批量转tif
def hdf2tif_batch(hdfFolder):
    #  获取文件夹内的文件名
    hdfNameList = os.listdir(hdfFolder)
    for i in range(len(hdfNameList)):
        #  判断当前文件是否为HDF文件
        if(os.path.splitext(hdfNameList[i])[1] == ".hdf"):
            hdfPath = hdfFolder+"/"+hdfNameList[i]
            #  gdal打开hdf数据集
            datasets = gdal.Open(hdfPath)
            
            #  获取hdf中的元数据
            Metadata = datasets.GetMetadata()
            #  获取四个角的维度
            Latitudes = Metadata["GRINGPOINTLATITUDE.1"]
            #  采用", "进行分割
            LatitudesList = Latitudes.split(", ")
            #  获取四个角的经度
            Longitude = Metadata["GRINGPOINTLONGITUDE.1"]
            #  采用", "进行分割
            LongitudeList = Longitude.split(", ")

            # 图像四个角的地理坐标
            GeoCoordinates = np.zeros((4, 2), dtype = "float32")
            GeoCoordinates[0] = np.array([float(LongitudeList[0]),float(LatitudesList[0])])
            GeoCoordinates[1] = np.array([float(LongitudeList[1]),float(LatitudesList[1])])
            GeoCoordinates[2] = np.array([float(LongitudeList[2]),float(LatitudesList[2])])
            GeoCoordinates[3] = np.array([float(LongitudeList[3]),float(LatitudesList[3])])

            #  列数
            Columns = float(Metadata["DATACOLUMNS"])
            #  行数
            Rows = float(Metadata["DATAROWS"])
            #  图像四个角的图像坐标
            PixelCoordinates = np.array([[0, 0],
                                         [Columns - 1, 0],
                                         [Columns - 1, Rows - 1],
                                         [0, Rows - 1]], dtype = "float32")
    
            #  计算仿射变换矩阵
            from scipy.optimize import leastsq
            def func(i):
                Transform0, Transform1, Transform2, Transform3, Transform4, Transform5 = i[0], i[1], i[2], i[3], i[4], i[5]
                return [Transform0 + PixelCoordinates[0][0] * Transform1 + PixelCoordinates[0][1] * Transform2 - GeoCoordinates[0][0],
                        Transform3 + PixelCoordinates[0][0] * Transform4 + PixelCoordinates[0][1] * Transform5 - GeoCoordinates[0][1],
                        Transform0 + PixelCoordinates[1][0] * Transform1 + PixelCoordinates[1][1] * Transform2 - GeoCoordinates[1][0],
                        Transform3 + PixelCoordinates[1][0] * Transform4 + PixelCoordinates[1][1] * Transform5 - GeoCoordinates[1][1],
                        Transform0 + PixelCoordinates[2][0] * Transform1 + PixelCoordinates[2][1] * Transform2 - GeoCoordinates[2][0],
                        Transform3 + PixelCoordinates[2][0] * Transform4 + PixelCoordinates[2][1] * Transform5 - GeoCoordinates[2][1],
                        Transform0 + PixelCoordinates[3][0] * Transform1 + PixelCoordinates[3][1] * Transform2 - GeoCoordinates[3][0],
                        Transform3 + PixelCoordinates[3][0] * Transform4 + PixelCoordinates[3][1] * Transform5 - GeoCoordinates[3][1]]
            #  最小二乘法求解
            GeoTransform = leastsq(func,np.asarray((1,1,1,1,1,1)))

            #  获取数据时间
            date = Metadata["RANGEBEGINNINGDATE"]

            #  第一个子数据集合,也就是NDVI数据
            DatasetNDVI = datasets.GetSubDatasets()[0][0]
            RasterNDVI = gdal.Open(DatasetNDVI)
            NDVI = RasterNDVI.ReadAsArray()
            
            TifName = date + ".tif"
            array2raster(TifName, GeoTransform[0], NDVI)
            print(TifName,"Saved successfully!")
            
hdf2tif_batch(r"E:Remote_Sensing_DataTVDIMOD13A32012")

输出：

2012-01-01.tif Saved successfully!
2012-02-01.tif Saved successfully!
2012-03-01.tif Saved successfully!
2012-04-01.tif Saved successfully!
2012-05-01.tif Saved successfully!
2012-06-01.tif Saved successfully!
2012-07-01.tif Saved successfully!
2012-08-01.tif Saved successfully!
2012-09-01.tif Saved successfully!
2012-10-01.tif Saved successfully!
2012-11-01.tif Saved successfully!
2012-12-01.tif Saved successfully!
weixin_39751871
关注
大数据之5V特点
无心⚡的博客
12-08 3139
5V特点是什么？Volume(大量)Variety(多样)Value(价值)Velocity(高速)Veracity(真实)**总结：大、多、值、快、信** Volume(大量) 数据的采集、计算、存储量都非常的庞大。数据的大小决定所考虑的数据价值和潜在的信息。 Variety(多样) 种类的来源多，数据类型的多样性。 Value(价值) 数据价值密度相对较低，犹如浪里淘金，合理运用大数据，以低成本创造高价值。 Velocity(高速) 数据增长速度快，处理速度也快，获取数据的速度也要快，时效性要求高。
大数据的5V特征
wyqwilliam的博客
06-08 2万+
IBM提出了大数据”5V”特点： 一、Volume：数据量大，包括采集、存储和计算的量都非常大。大数据的起始计量单位至少是P（1000个T）、E（100万个T）或Z（10亿个T）。 二、Variety：种类和来源多样化。包括结构化、半结构化和非结构化数据，具体表现为网络日志、音频、视频、图片、地理位置信息等等，多类型的数据对数据的处理能力提出了更高的要求。 三、Value：数据价值密度相对较低，或者说是浪里淘沙却又弥足珍贵。随着互联网以及物联网的广泛应用，信息感知无处不在，信息海量，但价值密度较低，如
大数据的 5V 特征
weixin_42511070的博客
07-02 631
这里写自定义目录标题欢迎使用Markdown编辑器新的改变功能快捷键合理的创建标题，有助于目录的生成如何改变文本的样式插入链接与图片如何插入一段漂亮的代码片生成一个适合你的列表创建一个表格设定内容居中、居左、居右SmartyPants创建一个自定义列表如何创建一个注脚注释也是必不可少的KaTeX数学公式新的甘特图功能，丰富你的文章UML 图表FLowchart流程图导出与导入导出导入 欢迎使用Markdown编辑器 你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Mar
大数据特点——5V
热门推荐
qigemantou的博客
06-14 4万+
大数据5V特点 我拆分理解了下，可以用五个字来概括 一、Volume： 大 数据的采集，计算，存储量都非常的庞大。 二、Variety： 多 种类和来源多样化。种类有：结构化、半结构化和非结构化数据等，常见的来源有：网络日志、音频、视频、图片等等。 三、Value： 值 数据价值密度相对较低，犹如浪里淘金，百炼成钢般才能获取到大量信息中的部分有价值的信息 四、Velocity： 快 数据...
大数据的5V特点
黄道婆的专栏
11-09 8792
大数据的5V特点（IBM提出）：Volume（大量）、Velocity（高速）、Variety（多样）、Value（低价值密度）、Veracity（真实性）。 特征 编辑 容量（Volume）：数据的大小决定所考虑的数据的价值和潜在的信息；[5] 种类（Variety）：数据类型的多样性；[5] 速度（Velocity）：指获得数据的速度；[5] 可变性（Variability）：妨碍了处理和有效地管理数据的过程。[5] 真实性（Veracity）：数据的质量...
GDAL.zip_GDAL_HDF GDAL_gdal modis_hdf4_hdf5
09-14
这扩展了GDAL的兼容性，使得处理MODIS（Moderate Resolution Imaging Spectroradiometer）等卫星遥感数据变得更加容易，因为MODIS数据通常以HDF格式发布。通过GDAL，用户可以进行数据的切片、重采样、投影转换等操作...
gdal.zip_GDAL_HDF View_gdal modis_hdf_hdf5
09-24
在这个名为“gdal.zip_GDAL_HDF View_gdal modis_hdf_hdf5”的压缩包中，包含了特别针对HDF（Hierarchical Data Format）格式的增强版GDAL库。HDF是一种用于存储和组织大量科学数据的文件格式，分为HDF4和HDF5两个...
gdal裁剪tif_tif遥感图像_GDAL裁剪_遥感影像_GDAL_gdal批量裁剪_
10-01
GDAL为遥感图像的处理提供了强大工具，通过Python脚本实现的批量裁剪极大地提高了工作效率。理解GDAL的基本操作和API，以及遥感图像的特性，是成功进行TIFF图像裁剪的关键。通过实践和优化，你可以定制化脚本来满足...
nctotiff.rar_dancexs1_gdal nc to tif_nc数据转tif_nc转txt_python .nc4
07-15
利用python与GDAL实现nc数据转成TIFF
hdf.rar_C hdf_GDAL_HDF GDAL_hdf_波段提取
09-19
在本项目中，我们将探讨如何使用C++和GDAL库来实现HDF文件的波段提取。 首先，理解GDAL与HDF的交互是至关重要的。GDAL提供了对HDF4和HDF5两种版本的支持，它们是两种不同的文件格式，但都可以通过GDAL进行读取和...
通过GDAL\OGR获取栅格矢量边界Python代码
09-11
Python代码，获取栅格矢量边界，文件说明：https://xcosmic.blog.csdn.net/article/details/108536771
批量删除GDAL裁剪是边缘产出的纯白或者黑白无用图片
07-13
专门批量删除GDAL裁剪是边缘产出的纯白或者黑白无用图片；代码中路径需要自己修改，有任何问题可以反馈给我，我会及时修改并联系您！第一上传资源，谢谢合作！
大数据的5v特征
Tiechitan的博客
07-08 1万+
1，数据体量大2，数据类型多样3，价值密度低，商业价值高4，要求输入输出的书读快
已知一个GeoTiff，得到其边界矢量的方法
yuezhilanyi的博客
11-28 1204
在制作GeoTiff的样本时，经常需要知道它的边界矢量，本文给出了一个简单、准确的方法。
利用GDAL库读取tif格式的高光谱波段的波长
weixin_38536617的博客
11-13 2424
利用Python中的GDAL库读取tif格式的高光谱波段的波长 因为文件比较大，都是以读字节或者是读头文件的方式读取高光谱各波段波长，怕自己忘了，所以记录一下 直接贴代码吧 dataset = gdal.Open(path)#path为高光谱文件所在路径 if dataset == None: print(path + "文件无法打开") im_width = dataset.RasterXSize # 列数 im_height = dataset.RasterYSize # 行数 im_bands
用python画花瓣-用python画花瓣
weixin_37988176的博客
11-01 5363
广告关闭2017年12月，云+社区对外发布，从最开始的技术博客到现在拥有多个社区产品。未来，我们一起乘风破浪，创造无限可能。p_line(t, n, step_length, step_angle) t.rt(step_angle 2) defpetal(t, r, angle): draws a 花瓣 using two arcs. for i in range(2):arc(t, r, ang...
Python-gdal读取卫星数据，并提取三个波段写入新的tif
weixin_39317090的博客
03-15 2600
#导入gdal包 from osgeo import gdal #读取遥感影像数据 path=r'E:\oil.tif' ds=gdal.Open(path) num=ds.RasterCount #获取影像波段数 #栅格矩阵的列数 im_width = ds.RasterXSize #栅格矩阵的行数 im_height = ds.RasterYSize #获取仿射矩阵，左上...
多波段TIF图像中选择指定几个波段导出 (ENVI)
Vanvan's Blog
04-11 7166
1、打开所要处理的tif图片； 2、如下图，在右边栏搜索Layer Stacking并双击打开 3、点击Import File，选择要处理的tiff 4、 选择Spectral Subset 5、选择指定范围波段，这里选择1-3，按ok结束 6、返回该界面，填写导出文件名称 7、按下OK导出 ......
python读取tif_在python中利用GDAL对tif文件进行读写的方法
weixin_39622419的博客
12-04 1126
利用GDAL库对tif影像进行读取示例代码默认波段为[B、G、R、NIR的顺序，且为四个波段]import gdaldef readTif(fileName):dataset = gdal.Open(fileName)if dataset == None:print(fileName+"文件无法打开")returnim_width = dataset.RasterXSize #栅格矩阵的列数im_...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

抵扣说明：

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

余额充值