01 前言
最近用ENVI IDL,觉得身心俱疲,一方面是学的不深,另一方面是关于IDL的资料太少了,基本上仅有少量的陈年博客(虽然写的不错),其余基本上来源于官方文档.所以还是做好两手准备,python处理遥感影像也要抓起来,毕竟python对于人工智能的专长之于遥感影像的冲击也比较大,现在漂亮国要是再不发几颗卫星国内遥感我们这批专业生就活不下去了,论文个个都在卷计算机算法人工智能,神卷。
废话不多说,我将以先讲述其参数,至于用法我打算以GPM Imerg early数据集为例将nc文件转为tif文件来说明rasterio的用法。
(至于为什么使用rasterio而不是gdal,有一定的考虑,我认为rasterio的语法更符合python的理念,而gdal的语法等等更偏向于C语言之类的)
02 参数说明
基本形式:
rasterio.open(fp, mode='r', driver=None, width=None, height=None, count=None, crs=None, transform=None, dtype=None, nodata=None, sharing=False, **kwargs)
rasterio.open函数用于创建一个DatasetReader或DatasetWriter对象,这两个对象分别用于读取和写入栅格数据。
fp: 需要打开的文件的路径(字符串);
mode: 字符串(可选参数,默认是r,只读模式), 打开文件的模式,包括四种,分别是'r', 'w', 'r+', 'w+',分别表示只读模式, 只写, 可读可写,可读可写模式。对于后面r+和w+,区别在于如果文件路径不存在,那么w+会新建一个,而r+则会报错,共同部分即可读又可写,但是需要注意的是如果文件路径存在即文件存在那么w+会将其中内容全部清空而r+会在其后内容继续追加内容。w和w+的区别就是w只可以写入不可以读,而w+既可以读也可以写入;
driver:文件的格式,一般在r和r+模式会省略,因为函数会自动获取其文件后缀判别。若是创建文件则需要指定该参数。若是创建GeoTIFF文件,则需要指定driver=GTiff。此处的driver与gdal类似,可以参照:Raster drivers — GDAL documentation获取不同文件格式的driver参数;
width:图像的宽度,也即图像栅格矩阵的列数;
height:图像的高度,也即图像栅格矩阵的行数;
count:图像的波段数;
crs:可以是字符串/字典/crs对象。
transform: 文件的仿射变换,可以是一个Affine对象或一个包含6个元素的列表或元组。这6个元素表示(x_size, skew_y, x_upper_left, skew_x, -y_size, y_upper_left), x为lon,y为lat.
你可以自行计算出来,对于旋转系数skew_x和skew_x填写0.0即可,一般你不会使用到它的。(我后面是通过其他函数计算因为我比较懒)
由于在我们遥感影像上,坐标系的原点在左上角点,所以输入的角点信息也在左上角的经纬度坐标,由于是左上角点所以它计算是从左往右,从上往下,因此x上的分辨率正常,但是y上的分辨率需要添上负号因为从上往下纬度在减小。
dtype:文件的数据类型,字符串或者np type都是可行的。
nodata:在矩阵中无效值的像素值,可以传入int,float,nan。
sharing:bool型。是否共享句柄,多线程应当避免。当处理大量数据时,操作系统可能会耗尽可用的文件描述符,这可能导致程序崩溃或无法打开更多文件。为了避免这个问题,Rasterio维护了一个共享句柄池,可以在多个地方重用这些句柄,从而减少了打开新文件的开销。
03 用法
好了基本上说清楚了 ,至于用法我就不详细赘述,时间有限,大家自行查看即可。
以下是关于如何处理将GPM Imerg Early的NC4文件处理为GeoTIFF文件。
import netCDF4 as nc import rasterio from rasterio.transform import from_origin # preparation in_path = r'F:\ExtremePrecipitation\data\GPM IMERG Early\3B-HHR-E.MS.MRG.3IMERG.20180312-S000000-E002959.0000.V06B.HDF5.SUB.nc4' out_path = r'F:\ExtremePrecipitation\TEMP\output.tif' # get the precip_dataset, lon_dataset, lat_dataset of the nc4 file dataset = nc.Dataset(in_path) # get the dataset_writer object precipitation = dataset.variables['precipitationCal'][0, :, :] # shape=(1, 630, 510), (channels, cols, rows)) lon = dataset.variables['lon'][:] # shape=(630,), this means the cols == 630 lat = dataset.variables['lat'][:] # shape=(510,), this means the rows == 510 # get the basic info of the precipitation rows = len(lat) cols = len(lon) lon_upper_left = min(lon) lat_upper_left = max(lat) lon_res = lon[1] - lon[0] # assume the lon and lat are equally spaced lat_res = lat[1] - lat[0] # write the precipitation data into a tif file with rasterio.open(out_path, 'w', driver='GTiff', height=rows, width=cols, count=1, dtype=precipitation.dtype, crs='+proj=latlong', transform=from_origin(lon_upper_left, lat_upper_left, lon_res, lat_res)) as dst: dst.write(precipitation, 1) # write the precipitation dataset into the first band of the tif file