HDF文件转Tif 温室气体 python代码

HDF文件转tif文件 ，python代码

这个问题网上找了好多，基本用处不大，我专门找了地理学的同学给做的，对于计算机的我来说确实是个难点，话不多说

hdf中所用到的三个数据，我用的是co2浓度 和经纬度，如图

我有两个代码，以下是代码1，

import os, sys
from osgeo import gdal,ogr,osr
import numpy as np
from pyhdf.SD import SD

os.environ['PROJ_LIB'] = r'D:\python\pythonData\venv\Lib\site-packages\osgeo\data\proj'
def readhdf(filename):

    file = SD(filename)   
    lon = (file.select('lon_grid')).get()
    lat = (file.select('lat_grid')).get()
    CO2Amount = (file.select('CO2Amount_grid')).get()
    lon_res = 1
    lat_res = -1
    return CO2Amount,lon_res,lat_res,lon,lat


def writetif(data,outname,geotransform):
    nl,ns = [data.shape[0],data.shape[1]] 
    bands = 1               
    driver = gdal.GetDriverByName("GTiff")
    out_tif = driver.Create(outname, ns, nl, bands, gdal.GDT_Float32)        
    out_tif.SetGeoTransform(geotransform)     
    srs = osr.SpatialReference()
    proj_type = 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]'
    out_tif.SetProjection(proj_type)  # 给新建图层赋予投影信息       
    out_tif.GetRasterBand(1).WriteArray(data)
    del out_tif

inpath = 'E:/work/hdftotif/hdf'
outpath = 'E:/work/hdftotif/tif'


dirs = os.listdir(inpath)  # 获取指定路径下的文件
file_name = []

for x in dirs:  # 循环读取路径下的文件并筛选输出
    if os.path.splitext(x)[1] == '.hdf' or os.path.splitext(x)[1] == '.HDF':   # 筛选.HDF文件
        file_name.append(x)
for k in range(len(file_name)):
    filename = inpath+'/'+file_name[k]
    CO2Amount,lon_res,lat_res,lon,lat = readhdf(filename)      
    outname = outpath + '/' + file_name[k][0:len(file_name[k])-4] + '.tif'    
    geotransform = (np.min(lon), lon_res, 0, np.max(lat), 0, lat_res)        
    writetif(CO2Amount,outname,geotransform)
    print('file finish',k)

代码2

from osgeo import gdal, osr
# import netCDF4 as nc
import numpy as np
from pathlib import Path
from collections import namedtuple
import sys
from pyhdf import SD
import os
import sys

# os.environ['GDAL_DATA'] = '/home/superres-master/superres_env_27/share/gdal'
os.environ['PROJ_LIB'] = r'D:\python\pythonData\venv\Lib\site-packages\osgeo\data\proj'

def read_hdf(filepath):
    hdf = SD.SD(str(hdf_p))
    Lat = hdf.select('lat_grid')[:]
    Lon = hdf.select('lon_grid')[:]

    Lon_min = Lon.min()
    Lat_min = Lat.min()
    Lon_max = Lon.max()
    Lat_max = Lat.max()

    # 分辨率计算
    lat_len = len(Lat)
    lon_len = len(Lon)
    lon_res = round((Lon_max - Lon_min) / (lon_len - 1), 4)
    lat_res = round((Lat_max - Lat_min) / (lat_len - 1), 4)

    geo_info = namedtuple(
        'geo_info', 'lon_min lat_max lon_res lat_res lon_len lat_len')
    info = geo_info(lon_min=Lon_min, lat_max=Lat_max, lon_res=lon_res,
                    lat_res=lat_res, lon_len=lon_len, lat_len=lat_len)
    return hdf, info


def write_hdf_to_tif(data, tifname, geo_info):
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(tifname, geo_info.lon_len,
                           geo_info.lat_len, 1, gdal.GDT_Float32)
    # 设置影像的显示范围
    #-lat_res一定要是-的
    geotransform = (geo_info.lon_min, geo_info.lon_res, 0,
                    geo_info.lat_max, 0, -(geo_info.lat_res))
    out_ds.SetGeoTransform(geotransform)
    # 获取地理坐标系统信息，用于选取需要的地理坐标系统
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
    out_ds.SetProjection(srs.ExportToWkt())  # 给新建图层赋予投影信息
    out_ds.GetRasterBand(1).SetNoDataValue(-999)

    out_ds.GetRasterBand(1).WriteArray(data)
    out_ds.FlushCache()
    del out_ds


# def make_raster(in_ds, fn, data, nodata=None):
#     driver = gdal.GetDriverByName('GTiff')
#     out_ds = driver.Create(
#         fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, gdal.GDT_Float32)
#     srs = osr.SpatialReference()
#     srs.ImportFromEPSG(4326)  # 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
#     out_ds.SetProjection(srs.ExportToWkt())

#     out_ds.SetGeoTransform(in_ds.GetGeoTransform())
#     out_band = out_ds.GetRasterBand(1)
#     if nodata is not None:
#         out_band.SetNoDataValue(nodata)
#     out_band.WriteArray(data)


hdf_root = Path('../data/rawhdf')

tif_root = Path('../data/tif')

for hdf_p in hdf_root.glob('*'):
    hdfName = hdf_p.name
    print('{} 正在处理中'.format(hdfName))

    hdf, geo_info = read_hdf(hdf_p)
    data = hdf.select('CO2Amount_grid')[:]

    stem_name = hdf_p.stem + '.tif'
    tif_p = tif_root / stem_name
    write_hdf_to_tif(data, str(tif_p), geo_info)

两个代码最终都可以生成tif