上次写《如何从中央气象台全国雷达拼图中提取dbz》的文章之后，就有很多朋友对于全国雷达拼图的投影信息给我提供了一些建议和帮助，后来同为大气所毕业的何爽爽师妹给我提供了一份她之前自己试出来的投影参数，于是我测了一下，吻合度很高。在她的参数的基础上，我尝试着找出了可以用于反算地理坐标的边界参数，最终实现了使用GDAL反解png地理坐标的目标，在此对何师妹表示感谢，下面我将我这道题的解题过程大致过程记录如下。
基于所猜测的投影参数绘制样例地图
为了测试所猜测的投影参数是否与原始图片的投影参数吻合，我们就需要使用猜测的投影参数绘制中国的国界线图，然后再用该图与原始图片进行比对，若二者在国界线上的吻合度足够高，我们就认为该投影参数是可用的。
基于上述的思路，我们首先需要根据猜测的投影参数绘制图片，根据何师妹提供的参数，原始数据为兰伯特正形投影，中央经度是110°E，两条标准纬线是30°N和62°N（对应的中央纬度是39°N），我根据该参数基于cartopy编写了以下代码：
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cnmaps import get_map, draw_map
east = 126.9
west = 83.6
north = 52.1
south = 13.9
fig = plt.figure(dpi=200)
ax = fig.add_subplot(111, projection=ccrs.LambertConformal(central_longitude=110, standard_parallels=[30,62], central_latitude=39))
draw_map(get_map('中国2'), color='r', linewidth=0.5)
ax.set_extent([west, east, south, north])
plt.savefig('sample.png')
效果图：
这里需要说明一下，上述代码中出现了东南西北的边界参数：
east = 126.9 west = 83.6 north = 52.1 south = 13.9
这是一个已经经过反复推敲和试验后得到的边界范围参数，该参数可以较好地与原始图片的尺寸完成匹配，但是注意，这个边界参数并非真实的四角坐标参数，不可以直接用于地理坐标系反解。
将样例地图和原始地图进行比对
我们将原始雷达图的底图和我们绘制的样例图同时输入Photoshop中，分开两个图层，将样例图放在上层并设为半透明状态，调整其大小使其尽量与原始图片的边界吻合。最终可以看到，二者的匹配度相当的高。
因此我们可以认为该投影参数是可用的，那么怎么进行坐标系反解呢？当然首先想到的就是GDAL，我们可以通过gdal_translate命令将PNG图片转为带有坐标信息的GeoTIFF格式，但是至少需要传入投影参数和边界坐标，投影参数目前已知，而边界坐标目前未知，需要我们找到在兰伯特正形投影坐标系下的四角坐标值（而不是地理坐标系的坐标）。好在我们可以在cartopy中直接获取该值，可以在绘图的过程中执行：
xmin, xmax = ax.get_xlim() ymin, ymax = ax.get_ylim()
注意上述代码需要在ax.set_extent执行之后再执行。这样我们可以获得的四角坐标是：左上：(xmin, ymax)，左下:（xmin, ymin），右上（xmax, ymax），右下（xmax, ymin）。而对于gdal_translate来说，它需要传入的-a_ullr参数，其实只有左上和右下的坐标参数，经过反复试验和比对，最终我确定的数值为 -3116700.0385205294, 1693154.412068708, 2017567.3399639379, -2737364.450775645。
GDAL大法反解PNG的地理坐标系并生成GeoTIFF
根据上面的步骤，我们已经能够获取用于反解地理坐标的足够的信息，可以下载一个全国雷达拼图的图片（例如名字叫做SEVP_AOC_RDCP_SLDAS_EBREF_ACHN_L88_PI_20210802214800001.PNG ），然后执行以下命令
gdal_translate -of GTiff -a_srs "+proj=lcc +lat_1=30 +lat_2=62 +lat_0=39.0 +lon_0=110 +a=6378137 +b=6378137 +units=m +no_defs" -a_ullr -3116700.0385205294 1693154.412068708 2017567.3399639379 -2737364.450775645 SEVP_AOC_RDCP_SLDAS_EBREF_ACHN_L88_PI_20210802214800001.PNG SEVP_AOC_RDCP_SLDAS_EBREF_ACHN_L88_PI_20210802214800001.tif
可以获得一个带有地理坐标参数的GeoTIFF文件，该文件可以直接用QGIS或ArcGIS打开，下图展示了用QGIS打开该图层，设为半透明后叠加在Google地图底图上的效果，由于带有投影和地理参数，QGIS在展示时自动将其转为了Web墨卡托投影（EPSG: 3857）效果：
可以看到其匹配效果相当的好。至此，我们基本已经解决了反解全国雷达图中最难的一关，其他投影转换、生成nc之类的工作都是很简单的事情了，后续我想让这个反解功能的实用性发挥得更好，因此计划用业余时间开发一套可以一键反解全国雷达图的Python扩展包（附带命令行工具）。感兴趣的可以持续关注一下。
近期评论