这个问题网上找了好多，基本用处不大，我专门找了地理学的同学给做的，对于计算机的我来说确实是个难点，话不多说
hdf中所用到的三个数据，我用的是co2浓度 和经纬度，如图
我有两个代码，以下是代码1，
import os, sysfrom osgeo import gdal,ogr,osrimport numpy as npfrom pyhdf.SD import SDos.environ['PROJ_LIB'] = r'D:pythonpythonDatavenvLibsite-packagesosgeodataproj'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,latdef 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_tifinpath = '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 ncimport numpy as npfrom pathlib import Pathfrom collections import namedtupleimport sysfrom pyhdf import SDimport osimport sys# os.environ['GDAL_DATA'] = '/home/superres-master/superres_env_27/share/gdal'os.environ['PROJ_LIB'] = r'D:pythonpythonDatavenvLibsite-packagesosgeodataproj'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, infodef 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