Python3.GDAL写入GeoTiff的一个小坑

从hdf文件提取的影像写入GeoTiff时，总是出错变成全是噪声的雪花图像该怎么办？我辛辛苦苦折腾了一个下午终于搞定了，但一点都开心不起来。

GDAL

GDAL是一个厉害的开源栅格空间数据转换库，ArcGIS、Google Earth等软件都用的它。像保存GeoTiff这种工作用它来做最合适不过了。

GDAL成为OSGeo的子项目后，作为一个Python的第三方库，GDAL被包含在osgeo这个名字空间里，所以安装时用的是

conda install gdal

使用的时候却是

from osgeo import gdal

复现问题

我现在有一个从hdf文件中用netCDF4读取的名为channel12的矩阵，表示10~54°N，70~140°E，分辨率为0.1°的一景图像，在Spyder的变量窗口显示这样

然后我用下面这段代码将它写入channel12.tif

from osgeo import gdal, osr

driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create('channel12.tif', 701, 441, 1, gdal.GDT_Float32)
dataset.SetMetadataItem('AREA_OR_POINT', 'Point')
dataset.SetGeoTransform([69.95, 0.1, 0, 54.05, 0, -0.1])
# 获取地理坐标系统信息，用于选取需要的地理坐标系统
sr = osr.SpatialReference()
sr.SetWellKnownGeogCS('WGS84')
# 给新建图层赋予投影信息
dataset.SetProjection(sr.ExportToWkt()) 
dataset.GetRasterBand(1).WriteArray(channel12)
dataset.FlushCache()  # 将数据写入硬盘
dataset = None

关于AREA_OR_POINT我最后讲，先讲大问题：这段代码可以正常运行并生成channel12.tif，但是用ArcMap打开它是这样的

检查过channel12矩阵中数值范围是191.1~319.5，用matplotlib画出来的图像也是正常的，所以问题出在写入GeoTiff这一步中

于是一个念头闪过：会不会是我装的GDAL有bug？因为conda环境的默认包管理策略比较保守，第三方库的版本经常比实际最新版本落后一点。这个愚蠢的念头浪费了我一个下午——我单独卸载重装GDAL过、建了一个虚拟python环境按PYPI的GDAL页面推荐的方式安装过、用whl文件安装过……这些尝试都没有解决问题！

我好难啊！

网上别人都没有这个问题，我用numpy创建的比这一景大的ndarray写入tif都没有问题，偏偏就是这一景图像会出现这可恶的雪花噪声问题！

我绝望地排查一切细节时，检查了一下channel12.dtype，结果是这样

这是个啥呀？！变量窗口明明显示的是float32啊！修改一下dtype

channel12 = channel12.astype(np.float32)

然后重新运行一次写入channel12.tif的代码，boom，成功了！！

AREA_OR_POINT

关于AREA_OR_POINT，是用来说明栅格究竟表示的是面积AREA还是点POINT的一个Metadata，默认是AREA。无论是哪个，在ArcMap中打开的效果是一样的。

但是相关的不得不提的是dataset.SetGeoTransform中第一和第三个参数分别表示最左上角（即第一个）栅格的左上角的经、纬度，第二和第六个参数分别表示从左到右的经度间隔和从上到下的纬度间隔。

本文中的dataset.SetGeoTransform表示：图像最左上角的经纬度为69.95°E，54.05°N。AREA_OR_POINT=POINT表示最左上角栅格值是属于栅格中心（70°E，55°N），若AREA_OR_POINT=AREA表示最左上角栅格值是属于栅格左上角（69.95°E，54.05°N）。

反思

造成这次浪费这么多时间的主要原因是细节知识不扎实：我只知道netCDF4读取的数据是numpy.ndarray，但是没有意识到numpy.ndarray中的数值类型是非常丰富的，并不像MATLAB那样简单只有几种常见的数值类型。

事实上，研究一下>f4确实也是float32，只不过是大端模式存储的，而在我的系统默认是小端存储，即<f4等价于np.float32如图

继续努力吧，这辈子都学不完哒