Python使用小教程03——tiff/grib/nc/hdf/txt数据的读写

工作需要，要从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

还有吗还有吗？没了吧。。。
嗯，那就酱吧