如何在自定义地图上不重采样绘制海岸线

Question

我想展示一张卫星图像（最好是用Python来做，不过其他方法也可以）。这张图像是一个浮点数参数P，大小是NxM，每个像素都有对应的经纬度信息（同样是NxM的大小）。所以我想要做到以下几点：

(1) 创建一个参数P的图像，并且要有对应的颜色比例尺。这个图像不应该被重新采样，所以它的大小应该保持在NxM。

(2) 在这个图像上显示海岸线。

目前，我可以用PIL来完成第(1)点。我也可以使用basemap库来显示图像和海岸线，但我不知道怎么在不改变投影的情况下，保持图像的原始大小NxM来实现。

补充说明：参数P本身并不包含海岸线的信息。只有像素的经纬度（lat, lon）可以用来叠加海岸线。海岸线的坐标可以从gshhs获取，比如说。gshhs实际上是在basemap库中使用的。

Answer 1

对于那些仍在寻找答案的朋友，我现在正在使用的方法（目的很相似——我在尝试测试卫星数据的地理定位）需要一个陆地掩模。

网上有很多陆地掩模的数据集，每个数据集都有不同的规则和特点。我在用Python处理netCDF4格式的数据，而我的陆地掩模是一个网格化的.nc数据集，其中海洋的值为1，陆地的值为0。

在处理我的卫星数据时，我会把每个纬度和经度的值乘以陆地掩模中每度的元素数量。在我的情况下，每度有120个元素，所以

lon_inds = (lons*120).astype(int)
lat_inds = (lats*120).astype(int)

更一般的写法是把120替换为

len(lons)/360
len(lats)/180

。如果使用numpy数组（这正是Python netCDF4模块的情况），这两种操作几乎可以瞬间完成。

接下来，我创建了一个自己的掩模：它必须和数据数组有相同的维度（对于不太了解卫星数据的人来说，数据、纬度和经度数组的维度都是一样的）：

my_mask = np.zeros(data.shape, dtype=int)

现在我们需要做的就是在掩模中替换海岸线的值。这是通过遍历lat_inds和lon_inds数组来完成的，查找陆地掩模中

landmask[lon_inds[i,j],lat_inds[i,j]]

的值，并将

mask[i,j]

的值改为1，如果任何邻居的值

landmask[lon_inds[i,j]-1,lat_inds[i,j]]
landmask[lon_inds[i,j]+1,lat_inds[i,j]]
landmask[lon_inds[i,j],lat_inds[i,j]-1]
landmask[lon_inds[i,j],lat_inds[i,j]+1]

不等于0（当然，通过添加对角线邻居的单元格可以生成更平滑的海岸线，但这通常不是必要的，因为你应该使用的陆地掩模数据集的空间分辨率应该比你的卫星数据更高）。

Answer 2

如果你只是想增强陆地和水域之间的边界，使用高通滤波器可能是个不错的选择。

比如，先从这张图开始，叫做Lena：

然后应用一个高通滤波器：

_{(来源: kiranisingh at sites.google.com)}

接着把高通滤波的结果叠加到原图上：

(更多细节和例子可以在这里找到)。

你可以在scipy中找到滤波器，具体可以查看这里。