工作需要,要从ec、cimiss、h8官网上下载并处理各种数据,涉及tiff、grib、nc格式数据的读写,这里做一个简单的总结
01 读数据
打开tiff主要是用的gdal(netCDF4也可以的)
def read_tiff(path):
"""
读取 TIFF 文件
"""
dataset = gdal.Open(path)
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
im_bands = dataset.RasterCount # 波段数
im_proj = dataset.GetProjection() # 获取投影信息
im_geotrans = dataset.GetGeoTransform() # 获取仿射矩阵信息
im_data = dataset.ReadAsArray(0, 0, im_width, im_height) # 获取数据
logging.info('read tiff success')
return im_data, im_width, im_height, im_bands, im_geotrans, im_proj
打开nc数据就用netCDF4,这里是处理的ec数据
def read_nc(path):
data = nc.Dataset(path)
lon = np.array(data['longitude'])
lat = np.array(data['latitude'])
time_num = len(np.array(data['time']))
st2 = np.array(data['skt'])
打开grib2数据比较麻烦,我的思路是wgrib2转nc转tiff,其实也有直接打开grib2的包的,pygrib,但我在win下调用不起来,就放弃了。。。这里是处理的cldas数据
def grib2_to_nc(grib_file, nc_file):
cmd = 'wgrib2.exe ' + grib_file + ' -netcdf ' + nc_file
os.system(cmd)
print("done!")
打开HDF文件用netCDF4(nc)就好,这里处理的FY4A数据,返回的是通道亮温值
def get_bt(fdi_file, channel, invalid_value):
fdi = nc.Dataset(fdi_file, 'r', format='NETCDF4_CLASSIC')
channel_name = 'NOMChannel' + str(channel)
cal_name = 'CALChannel' + str(channel)
bt_lut = np.array(fdi.variables[cal_name])
bt_lut = np.append(bt_lut, np.nan)
dn = np.array(fdi.variables[channel_name])
dn[(dn == invalid_value) | (dn >= np.size(bt_lut))] = bt_lut.shape[0] - 1
bt = bt_lut[dn]
fdi.close()
return bt
打开txt数据,把它转成二维数组,格点数据
def txt_to_array(path):
f = open(path)
line = f.readline()
line_list = line.split("|")
line_array = [float(i) for i in line_list]
data = np.array(line_array).reshape(n_lat, n_lon)
return data
02 写数据
把已经读好的二维数组写成tiff
单波段:如果你的需求只是WGS84的话,下面这一段就可以了,但是如果你还需要别的投影信息,那就要再研究一下仿射矩阵的参数和投影的设置了
def write_to_tiff(data, out_file):
ds = gdal.GetDriverByName('Gtiff').Create(out_file, n_lon, n_lat, 1, gdal.GDT_Float32)
geotransform = (lon_min, lon_res, 0, lat_max, 0, -lat_res)
ds.SetGeoTransform(geotransform)
srs = osr.SpatialReference() #获取地理坐标系统信息,用于选取需要的地理坐标系统
print(type(srs))
print(srs)
srs.ImportFromEPSG(4326) # 定义输出的坐标系为"WGS 84",AUTHORITY["EPSG","4326"]
ds.SetProjection(srs.ExportToWkt()) # 给新建图层赋予投影信息
ds.GetRasterBand(1).WriteArray(data) # 将数据写入内存,此时没有写入硬盘
del ds
多波段:还可以通过SetNoDataValue设置no_data
def write_tiff(im_data, im_width, im_height, im_bands, im_geotrans, im_proj, datatype = gdal.GDT_Byte,
out_path=None, return_mode='TIFF'):
if out_path:
dataset = gdal.GetDriverByName('GTiff').Create(out_path, im_width, im_height, im_bands, datatype)
else:
dataset = gdal.GetDriverByName('MEM').Create('', im_width, im_height, im_bands, datatype)
# 写入数据
if dataset is not None:
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(im_proj)
# 写入矩阵
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
还有吗还有吗?没了吧。。。
嗯,那就酱吧