使用卷积窗口隔离沿海网格（python）(Isolating coastal grids using convolving window (python))

使用卷积窗口隔离沿海网格（python）(Isolating coastal grids using convolving window (python))

我试图从数据集中分离出沿海网格。 我能够通过单独检查陆地网格的每个邻居来做到这一点，如果8个邻居中的一个是海洋网格，那么我正在保存网格。 这是我写的函数：

def coastalGrids(ds): grid=~out[0,:,:].mask #out is a 3D masked array (time,lat,lon) #At each grid where there is land, check all its 8 neighbors, #and if any of them are ocean, save the grid to the dataset coastline=np.zeros(grid.shape, dtype=bool) for i in range(1,len(lat)-1): for j in range(1,len(lon)-1): if grid[i,j]==True: if (grid[i+1,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i+1,j]!=True): coastline[i,j]= grid[i,j] elif (grid[i+1,j-1]!=True): coastline[i,j]= grid[i,j] elif (grid[i,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i,j-1]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j-1]!=True): coastline[i,j]= grid[i,j] return coastline

我想知道是否：

使用scipy的卷积窗口函数，有这样一个不那么丑陋的方式吗？ 扩展这样的功能，有没有办法将陆地网格隔离半径10英里的海岸线？

谢谢！

I'm trying to isolate coastal grids from a dataset. I am able to do it by individually checking each neighbor of a land grid, and if one of the 8 neighbors is an ocean grid, then I'm saving the grid. Here's the function I've written:

def coastalGrids(ds): grid=~out[0,:,:].mask #out is a 3D masked array (time,lat,lon) #At each grid where there is land, check all its 8 neighbors, #and if any of them are ocean, save the grid to the dataset coastline=np.zeros(grid.shape, dtype=bool) for i in range(1,len(lat)-1): for j in range(1,len(lon)-1): if grid[i,j]==True: if (grid[i+1,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i+1,j]!=True): coastline[i,j]= grid[i,j] elif (grid[i+1,j-1]!=True): coastline[i,j]= grid[i,j] elif (grid[i,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i,j-1]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j+1]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j]!=True): coastline[i,j]= grid[i,j] elif (grid[i-1,j-1]!=True): coastline[i,j]= grid[i,j] return coastline

I wish to know if:

There is a less ugly way of doing this, say using scipy's convolving windows functions? Extending such function, is there a way to isolate land grids by a radius of say 10 grids from the coastline?

Thanks!

最满意答案

使用图像形态学算子

您当前正在做的事情相当于原始布尔数组与其二进制侵蚀之间的差异。

这与形态梯度密切相关，但它的定义略有不同。

无论如何，假设我们有一个非常简单的岛屿：

import numpy as np import matplotlib.pyplot as plt y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 fig, ax = plt.subplots() ax.pcolormesh(land, cmap='gray', edgecolor='gray', antialiased=True) plt.show()

我们可以通过侵蚀该岛来计算海岸线：

import scipy.ndimage erosion = scipy.ndimage.binary_erosion(land)

然后看看哪里有什么区别：

coast = land != erosion

---

正方形与对角线连接

默认情况下，这使用“正方形”连接。 换句话说，它不会将对角线视为触摸。 默认结构（aka“footprint”）看起来像：

[[0, 1, 0], [1, 1, 1], [0, 1, 0]]

在你的代码中，你假设对角线是触摸的。 在这种情况下，您需要“完整”连接，以及如下所示的结构：

[[1, 1, 1], [1, 1, 1], [1, 1, 1]]

为此，我们将指定类似于：

erosion = scipy.ndimage.binary_erosion(land, structure=np.ones((3,3)))

作为一个完整的例子：

import numpy as np import matplotlib.pyplot as plt import scipy.ndimage y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 erosion = scipy.ndimage.binary_erosion(land, structure=np.ones((3,3))) coast = land != erosion fig, ax = plt.subplots() ax.pcolormesh(coast, cmap='gray', edgecolor='gray', antialiased=True) plt.show()

使用形态学渐变

您也可以考虑使用形态学梯度算子。 这是给定输入和连接占用空间的二进制扩张和二进制侵蚀之间的差异。

在你的情况下，它还将包括边界陆地的海洋像素以及与海洋交界的陆地像素。 实际上，它会给你更厚的边界。

举个例子：

import numpy as np import matplotlib.pyplot as plt import scipy.ndimage y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 coast = scipy.ndimage.morphological_gradient(land, footprint=[[0, 1, 0], [1, 1, 1], [0, 1, 0]]) fig, ax = plt.subplots() ax.pcolormesh(coast, cmap='gray', edgecolor='gray', antialiased=True) plt.show()

Using Image Morphology Operators

What you're currently doing is equivalent to the difference between the original boolean array and its binary erosion.

This is closely related to the morphological gradient, but it's a slightly different definition.

At any rate, let's say we have a very simple island:

import numpy as np import matplotlib.pyplot as plt y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 fig, ax = plt.subplots() ax.pcolormesh(land, cmap='gray', edgecolor='gray', antialiased=True) plt.show()

We can calculate the coastline as you've defined it by eroding the island:

import scipy.ndimage erosion = scipy.ndimage.binary_erosion(land)

And then seeing where there are any differences:

coast = land != erosion

---

Square vs. Diagonal Connectivity

By default, this uses "square" connectivity. In other words, it doesn't count diagonals as touching. The default structure (a.k.a. "footprint") looks like:

[[0, 1, 0], [1, 1, 1], [0, 1, 0]]

In your code, you're assuming diagonals are touching. In that case, you want "full" connectivity, and a structure that looks like:

[[1, 1, 1], [1, 1, 1], [1, 1, 1]]

To do that, we'd specify something similar to:

erosion = scipy.ndimage.binary_erosion(land, structure=np.ones((3,3)))

As a complete example:

import numpy as np import matplotlib.pyplot as plt import scipy.ndimage y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 erosion = scipy.ndimage.binary_erosion(land, structure=np.ones((3,3))) coast = land != erosion fig, ax = plt.subplots() ax.pcolormesh(coast, cmap='gray', edgecolor='gray', antialiased=True) plt.show()

Using a Morphological Gradient

You might also consider using the morphological gradient operator. It's the difference between the binary dilation and the binary erosion for a given input and connectivity footprint.

In your case, it would also include the sea pixels that border land, as well as the land pixels that border sea. Effectively, it will give you thicker borders.

As an example:

import numpy as np import matplotlib.pyplot as plt import scipy.ndimage y, x = np.mgrid[-10:10:20j, -10:10:20j] land = np.hypot(x, y) < 7 coast = scipy.ndimage.morphological_gradient(land, footprint=[[0, 1, 0], [1, 1, 1], [0, 1, 0]]) fig, ax = plt.subplots() ax.pcolormesh(coast, cmap='gray', edgecolor='gray', antialiased=True) plt.show()