python实现hdf/h5文件转tif（可批量） |
您所在的位置：网站首页 › h5转tif › python实现hdf/h5文件转tif（可批量） |
|
所用数据：AMSR_U2_L3_DailySnow_B02_20121230.he5 雪水当量数据产品 环境：python2.7 并安装h5py模块 H5数据介绍： HDF(Hierarchical Data Format)，设计用于存储和组织大量数据的文件格式。 h5文件中有两个核心的概念：组“group”和数据集“dataset”。组可以理解为文件夹，数据集理解为文件，文件夹中可以递归地包含文件夹、文件等。 dataset ：简单来讲类似数组组织形式的数据集合，像 numpy 数组一样工作，一个dataset即一个numpy.ndarray。具体的dataset可以是图像、表格，甚至是pdf文件和excel。group：包含了其它 dataset(数组) 和 其它 group ，像字典一样工作。一、查看数据信息 用python查看没搞明白，我是用matlab看的 二、转tif 查阅了一些文章和自己摸索，发现以下两种方法都可以 方法一： 我用之前nc转tif的代码改写的 先看函数 def dayextract_hdf(hdf_file, output_dir): fp = h5py.File(hdf_file, 'r') Lat = fp['HDFEOS']['GRIDS']['Northern Hemisphere']["YDim"][:] Lon = fp['HDFEOS']['GRIDS']['Northern Hemisphere']["XDim"][:] # 影像的左上角和右下角坐标 # LonMin, LatMax, LonMax, LatMin = [Lon.min(), Lat.max(), Lon.max(), Lat.min()] # 分辨率计算 N_Lat = len(Lat) N_Lon = len(Lon) Lon_Res = (LonMax - LonMin) / (float(N_Lon) - 1) Lat_Res = (LatMax - LatMin) / (float(N_Lat) - 1) cols = N_Lon rows = N_Lat # 设置影像的显示范围 # -Lat_Res一定要是-的 # geo = (LonMin, Lon_Res, 0, LatMax, 0, -Lat_Res) geo = (-9036843.07384, 25067.52586, 0, 9036843.07384, 0, -25067.52586) # 构造projection src_srs = osr.SpatialReference() src_srs.ImportFromEPSG(3408) #查看文件投影代码为3408 src_srs_wkt = src_srs.ExportToWkt() # 给新建图层赋予投影信息 dayarray = fp['HDFEOS']['GRIDS']['Northern Hemisphere']['Data Fields']['SWE_NorthernDaily'][:] # print(len(dayarray)) out_file =os.path.join(output_dir, hdf_file[11:-4] + ".tif") print(out_file) arcpy.env.overwriteOutput = 1 # 输出文件夹里面已经有内容的，就覆盖掉 arcpy.CheckOutExtension("Spatial") # 生成tif，写入 driver = gdal.GetDriverByName("Gtiff") outdataset = driver.Create(out_file, cols, rows, 1, gdal.GDT_Float32) outdataset.SetGeoTransform(geo) outdataset.SetProjection(src_srs_wkt) band = outdataset.GetRasterBand(1) band.WriteArray(data_array) fp.close()之后再批量操作 # 读取所有HDF数据 hdf_file = r"D:\test\nc" tiff_dir = r"D:\test\tif" hdfFileList = getFileName(hdf_file, {'.he5'}) for hdfFile in hdfFileList: # nc数据转换为tif数据 print(hdfFile) hdf_to_tif(hdfFile, tiff_dir)这里如果xy或是经纬度读取有问题的话，可以直接把h5文件拖进arcmap点击图层属性查看图像的范围坐标，可以手动输入geo。同理，rows和cols也可以在arcmap中查看直接赋值 方法二： #-*- coding: UTF-8 -*- # import arcpy import glob import os import h5py import xlrd arcpy.env.overwriteOutput = 1 #输出文件夹里面已经有内容的，就覆盖掉 arcpy.CheckOutExtension("Spatial") #检查ArcGIS扩展模块是否可用 inPath = 'D:\\test\\nc\\' outDir = 'D:\\test\\tiff\\' def h5_to_tiff(inPath, outDir): file = os.listdir(inPath) # 读取路径中所有文件，以ls来表示，常用file for i in file: # 遍历ls中的文件 arcpy.env.workspace = inPath + i # 脚本中最为常用的环境变量设置就是arcpy.env.workspace，该变量用于定义当前脚本的工作目录（或者称为工作空间） # print(arcpy.env.workspace) arcpy.env.scratchWorkspace = inPath + i # （为什么定义临时空间变量作用暂不明） hdfList = arcpy.ListRasters('*', 'he5') # 按名称和栅格类型返回工作空间中的栅格列表(遍历指定类型的文件)，需先设置工作空间环境 # 判断存储最终图像的文件夹是否存在，是则存储，否则创建并存储；注意if和else后一定要加冒号，格式要对齐 if os.path.exists(outDir + i[:-4]): for hdf in hdfList: # 此Hdf文件有两个波段，全部输出，相应的ExtractSubDataset_management 中也要提取对应波段 for number in range(0, 4, 2): # 有几个波段就写几个,我这里需要输出第1个和第3个波段数据 outPath = outDir + str(i[:-4]) + '\\' # 根据对 HDF 数据集创建新栅格数据集,语法 ExtractSubDataset(in_raster, out_raster, {subdataset_index}) out = arcpy.ExtractSubDataset_management(hdf, outPath + hdf[0:-4] + str(number) + ".tif", number) else: os.makedirs(outDir + i[:-4]) for hdf in hdfList: for number in range(0, 4, 2): outPath = outDir + str(i[:-4]) + '\\' out = arcpy.ExtractSubDataset_management(hdf, outPath + hdf[0:-4] + str(number) + ".tif", number) print(outPath) h5_to_tiff(inPath, outDir)这段代码参考： https://blog.csdn.net/qq_38882446/article/details/103461759 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |