01 警告说明
首先贴出相关代码:
out_file_name = 'Rs_{:4.0f}{:02.0f}.tiff'.format(year, month)
out_path = os.path.join(out_dir, out_file_name)
mem_driver = gdal.GetDriverByName('MEM')
mem_ds = mem_driver.Create('', len(lon), len(lat), 1, gdal.GDT_Float32)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
mem_ds.SetProjection(srs.ExportToWkt())
mem_ds.SetGeoTransform(geo_transform)
mem_ds.GetRasterBand(1).WriteArray(cur_Rs)
mem_ds.GetRasterBand(1).SetNoDataValue(-9999) # 设置无效值
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(out_path, len(lon), len(lat), 1, gdal.GDT_Float32)
ds.SetGeoTransform(geo_transform)
ds.SetProjection(srs.ExportToWkt())
gdal.Warp(ds, mem_ds, cropToCutline=True, cutlineDSName=mask_path, xRes=out_res, yRes=out_res,
resampleAlg=gdal.GRA_Cubic, srcNodata=-9999, dstNodata=-9999)
警告代码在gdal.Warp
函数中产生。
警告: All options related to creation ignored in update mode
,很明显,正如它所说,所有和创建相关的options
将会在更新模式下忽略。
这里的options
根据gdal
源代码中:
02 解决
应该可以理解为gdal.Warp
函数中的相关参数,换言之,我们指定的参数例如前面代码中的:
gdal.Warp(ds, mem_ds, cropToCutline=True, cutlineDSName=mask_path, xRes=out_res, yRes=out_res,
resampleAlg=gdal.GRA_Cubic, srcNodata=-9999, dstNodata=-9999)
也就是说参数应该发生了某种冲突,导致有些参数无法起作用而忽略了,下面是该警告导致的输出影像显示:
可以发现,至少cropToCutline
、dstNodata
两个参数是没有正常作用的,掩膜似乎是作用的。周围黑色背景均为0,说明dstNodata
没有起作用,而且似乎不应该是0,因为我确定我的栅格数据中不存在0值,无效值也是赋值的是-65535。另外影像周围很多多余黑色背景,说明它没有将shp矩形范围外的栅格矩阵丢弃,因此cropToCutline
也是失效的。
调试发现,确实自己该的。gdal.Warp
的参数和前面的代码存在矛盾导致的警告,进而导致错误影像产生。
在代码中:
ds = driver.Create(out_path, len(lon), len(lat), 1, gdal.GDT_Float32)
ds.SetGeoTransform(geo_transform)
ds.SetProjection(srs.ExportToWkt())
一开始我是不存在传入仿射参数和投影信息的,但是命令行警告:
ERROR 1: Unable to compute a transformation between pixel/line and georeferenced coordinates for E:\FeaturesTargets\uniform\Rs\Rs_198307.tiff. There is no affine transformation and no GCPs. Specify transformation option DST_METHOD=NO_GEOTRANSFORM to bypass this check.
于是我顺手给加上了,但是投影信息是正常的,而传入的geo_transform
却是和后面的掩膜相矛盾,掩膜后的栅格矩阵必不可能是这个仿射参数。于是我想起来了,原来只需要传入输出路径即可,如下:
out_file_name= 'Rs_{:4.0f}{:02.0f}.tiff'.format(year, month)
out_path = os.path.join(out_dir, out_file_name)
mem_driver = gdal.GetDriverByName('MEM')
mem_ds = mem_driver.Create('', len(lon), len(lat), 1, gdal.GDT_Float32)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)
mem_ds.SetProjection(srs.ExportToWkt())
mem_ds.SetGeoTransform(geo_transform)
mem_ds.GetRasterBand(1).WriteArray(cur_Rs)
mem_ds.GetRasterBand(1).SetNoDataValue(-9999) # 设置无效值
gdal.Warp(out_path, mem_ds, cropToCutline=True, cutlineDSName=mask_path, xRes=out_res, yRes=out_res,
resampleAlg=gdal.GRA_Cubic, srcNodata=-9999, dstNodata=-9999)
无需建立两个ds对象出来。
但是在调试过程中,我一开始认为是在MEM
中创建影像导致的问题,于是将代码修改为:
mem_driver = gdal.GetDriverByName('GTiff')
mem_ds = mem_driver.Create(os.path.join(out_dir, 'temp.tiff'), len(lon), len(lat), 1, gdal.GDT_Float32)
结果发现输出的temp.tiff
文件在GIS中显示异常:
对的,就是完全空的状态,我在代码中使用matplotlib
显示如下:
说明应该是gdal在存入的时候导致的错误,最后找了很久没有发现原因,本来不打算研究了,因为虽然中间数据虽然存在显示问题,但是由于最后的输出还是正常的如下(我已经修改了此前的错误关于geo_transform
):
当然,结局是好的,最后发现的原因应该是python的垃圾回收机制这方面的原因,换句话就是代码执行到这里时
mem_ds.GetRasterBand(1).SetNoDataValue(-9999)
我们想当然理解为数据集应该传入到了硬盘中,但实际上它还在内存中。这是很正常的,我以为它会在代码结束后自动将内存中的数据写入到硬盘,但是实际上并没有。
往常我是一般会写:
mem_ds = None
或者
mem_ds.FlushCache()
前者通过python的垃圾回收机制将内存中的东西写入到硬盘后再进行回收,后者则是gdal
自带方法有相关写入内存的函数。
但是由于有时我发现写不写似乎并没有导致一些异常发生,且就算有异常一般是发生在命令行/jupyter/调试过程中,这个过程我认为是因为写入数据实际上还在内存,当我结束或者关闭时它会自动写入,事实上它确实如我所想这么实现了。因此我误以为所有数据集写入其实无需进行写入硬盘操作,因为当程序结束等会自动进行,但是事实证明我错了,或许在更为复杂的情况会导致没有来得及写入硬盘导致输出结果错误的情况发生。
所以,养成习惯,显式调用它们吧,这些错误有时非常致命当你不清楚它时。