从包含在 shapefile 边界内的 netcdf 文件中提取数据

Posted 2023-03-23

【中文标题】从包含在 shapefile 边界内的 netcdf 文件中提取数据

【英文标题】：extract data from netcdf file contained within a shapefile's boundaries

【发布时间】：2020-05-30 15:38:41

【问题描述】：

我有以下shapefile 和netcdf file。

我想从 netcdf 文件中提取包含在 shapefile 边界内的数据。

您对我如何实现这一点有什么建议吗？

shapefile 对应于 SREX 区域 11 North Europe (NEU)，netcdf 文件是 CMIP6 气候模型数据输出的示例（ua 变量）。我想要的输出必须是 netcdf 格式。

---

更新

到目前为止，我尝试使用 NCL 和 CDO 创建一个 netcdf 掩码，并将此掩码应用于原始 netcdf 数据集。下面是步骤（和NCL scripts）：

#################
## remove plev dimension from netcdf file
cdo --reduce_dim -copy nc_file.nc nc_file2.nc

## convert longitude to -180, 180
cdo sellonlatbox,-180,180,-90,90 nc_file2.nc nc_file3.nc

## create mask 
ncl create_nc_mask.ncl

## apply mask
cdo div nc_file3.nc shape1_mask.nc nc_file4.nc 
#################

输出几乎是正确的。见下图。但是 shapefile 的南部边界（SREX 11，NEU）没有被正确捕获。所以我想生成 netcdf 掩码的 NCL 脚本有问题。

【问题讨论】：

您是专门寻找 nco/cdo/ncl 解决方案，还是也欢迎其他语言（例如 Python）？

感谢@Bart。欢迎任何语言:)

多么慷慨的赏金????

使用 GDAL 的解决方案很简单。将您的 NetCDF 文件视为“栅格”（用 GIS 的说法）。您必须使用与 Shapefile 相同的 CRS 转换为 GeoTIFF。也是微不足道的。见：gis.stackexchange.com/questions/118236/…

【参考方案1】：

重新使用一些旧的脚本/代码，我很快想出了一个 Python 解决方案。它基本上只是遍历所有网格点，并检查每个网格点是在形状文件中的多边形内部还是外部。结果是变量mask（带有True/False 的数组），可用于屏蔽您的NetCDF 变量。

注意：这使用Numba（所有@jit 行）来加速代码，尽管在这种情况下这并不是真正必要的。如果您没有 Numba，您可以将它们注释掉。

import matplotlib.pyplot as pl
import netCDF4 as nc4
import numpy as np
import fiona
from numba import jit

@jit(nopython=True, nogil=True)
def distance(x1, y1, x2, y2):
    """
    Calculate distance from (x1,y1) to (x2,y2)
    """
    return ((x1-x2)**2 + (y1-y2)**2)**0.5

@jit(nopython=True, nogil=True)
def point_is_on_line(x, y, x1, y1, x2, y2):
    """
    Check whether point (x,y) is on line (x1,y1) to (x2,y2)
    """

    d1 = distance(x,  y,  x1, y1)
    d2 = distance(x,  y,  x2, y2)
    d3 = distance(x1, y1, x2, y2)

    eps = 1e-12
    return np.abs((d1+d2)-d3) < eps

@jit(nopython=True, nogil=True)
def is_left(xp, yp, x0, y0, x1, y1):
    """
    Check whether point (xp,yp) is left of line segment ((x0,y0) to (x1,y1))
    returns:  >0 if left of line, 0 if on line, <0 if right of line
    """

    return (x1-x0) * (yp-y0) - (xp-x0) * (y1-y0)

@jit(nopython=True, nogil=True)
def is_inside(xp, yp, x_set, y_set, size):
    """
    Given location (xp,yp) and set of line segments (x_set, y_set), determine
    whether (xp,yp) is inside polygon.
    """

    # First simple check on bounds
    if (xp < x_set.min() or xp > x_set.max() or yp < y_set.min() or yp > y_set.max()):
        return False

    wn = 0
    for i in range(size-1):

        # Second check: see if point exactly on line segment:
        if point_is_on_line(xp, yp, x_set[i], y_set[i], x_set[i+1], y_set[i+1]):
            return False

        # Calculate winding number
        if (y_set[i] <= yp):
            if (y_set[i+1] > yp):
                if (is_left(xp, yp, x_set[i], y_set[i], x_set[i+1], y_set[i+1]) > 0):
                    wn += 1
        else:
            if (y_set[i+1] <= yp):
                if (is_left(xp, yp, x_set[i], y_set[i], x_set[i+1], y_set[i+1]) < 0):
                    wn -= 1

    if wn == 0:
        return False
    else:
        return True

@jit(nopython=True, nogil=True)
def calc_mask(mask, lon, lat, shp_lon, shp_lat):
    """
    Calculate mask of grid points which are inside `shp_lon, shp_lat`
    """

    for j in range(lat.size):    
        for i in range(lon.size):
            if is_inside(lon[i], lat[j], shp_lon, shp_lat, shp_lon.size):
                mask[j,i] = True


if __name__ == '__main__':

    # Selection of time and level:
    time = 0
    plev = 0

    # Read NetCDF variables, shifting the longitudes
    # from 0-360 to -180,180, like the shape file:
    nc = nc4.Dataset('nc_file.nc')
    nc_lon = nc.variables['lon'][:]-180.
    nc_lat = nc.variables['lat'][:]
    nc_ua  = nc.variables['ua'][time,plev,:,:]

    # Read shapefile and first feature
    fc = fiona.open("shape1.shp")
    feature = next(iter(fc))

    # Extract array of lat/lon coordinates:
    coords = feature['geometry']['coordinates'][0]
    shp_lon = np.array(coords)[:,0]
    shp_lat = np.array(coords)[:,1]

    # Calculate mask
    mask = np.zeros_like(nc_ua, dtype=bool)
    calc_mask(mask, nc_lon, nc_lat, shp_lon, shp_lat)

    # Mask the data array
    nc_ua_masked = np.ma.masked_where(~mask, nc_ua)

    # Plot!
    pl.figure(figsize=(8,4))
    pl.subplot(121)
    pl.pcolormesh(nc_lon, nc_lat, nc_ua, vmin=-40, vmax=105)
    pl.xlim(-20, 50)
    pl.ylim(40, 80)

    pl.subplot(122)
    pl.pcolormesh(nc_lon, nc_lat, nc_ua_masked, vmin=-40, vmax=105)
    pl.xlim(-20, 50)
    pl.ylim(40, 80)

    pl.tight_layout()

编辑

要将掩码写入 NetCDF，可以使用以下内容：

nc_out = nc4.Dataset('mask.nc', 'w')
nc_out.createDimension('lon', nc_lon.size)
nc_out.createDimension('lat', nc_lat.size)

nc_mask_out = nc_out.createVariable('mask', 'i2', ('lat','lon'))
nc_lon_out = nc_out.createVariable('lon', 'f8', ('lon'))
nc_lat_out = nc_out.createVariable('lat', 'f8', ('lat'))

nc_mask_out[:,:] = mask[:,:]  # Or ~mask to reverse it
nc_lon_out[:] = nc_lon[:]     # With +180 if needed
nc_lat_out[:] = nc_lat[:]

nc_out.close()

【讨论】：

谢谢。有用。如何将 nc_ua_masked、nc_lat 和 nc_lon 变量保存为唯一的 netcdf 文件？

也...nc_ua_masked 数组不是具有每日 ua 值的完整数组。我没有解决我的问题。

即使将掩码保存为带有 lat、long 和 TRUE/FALSE 变量的 netcdf 文件也会很好

我添加了一个如何将掩码写入 NetCDF 的示例。它使用short 将掩码写为0-1，NetCDF4 似乎不支持布尔值（？）。作为替代方案，也可以在所有级别和时间步长上（在 NetCDF 文件的副本上）应用掩码。

【参考方案2】：

到目前为止，我想出了这个（我知道它不是完整的解决方案）：

1) 要打开 shapefile 和 nc 文件，您需要安装两个包：

pip3 install pyshp
pip3 install netCDF4

2) 那么这就是你在 python 中导入它们的方式：

import shapefile
from netCDF4 import Dataset

3) 从 shapefile 读取数据：

with shapefile.Reader("shape1.dbf") as dbf:
    print(f'dbf\n')
    print(f'bounding box: dbf.bbox')
    shapes = dbf.shapes()
    print(f'points: shapes[0].points')
    print(f'parts: shapes[0].parts')
    print(f'fields: dbf.fields')
    records = dbf.records()
    dct = records[0].as_dict()
    print(f'record: dct')

这会给你输出：

shapefile Reader
    1 shapes (type 'POLYGON')
    1 records (4 fields)

bounding box: [-10.0, 48.0, 40.0, 75.0]
points: [(-10.0, 48.0), (-10.0, 75.0), (40.0, 75.0), (40.0, 61.3), (-10.0, 48.0)]
parts: [0]
fields: [('DeletionFlag', 'C', 1, 0), ['NAME', 'C', 40, 0], ['LAB', 'C', 40, 0], ['USAGE', 'C', 40, 0]]
record: 'NAME': 'North Europe [NEU:11]', 'LAB': 'NEU', 'USAGE': 'land'

4) 读取nc文件：

nc_fdata = Dataset('nc_file.nc', 'r')

5) 使用这个辅助函数来查看里面的内容：

def ncdump(nc_fid, verb=True):
    def print_ncattr(key):
        try:
            print("\t\ttype:", repr(nc_fid.variables[key].dtype))
            for ncattr in nc_fid.variables[key].ncattrs():
                print('\t\t%s:' % ncattr, repr(nc_fid.variables[key].getncattr(ncattr)))
        except KeyError:
            print("\t\tWARNING: %s does not contain variable attributes" % key)

    # NetCDF global attributes
    nc_attrs = nc_fid.ncattrs()
    if verb:
        print("NetCDF Global Attributes:")
        for nc_attr in nc_attrs:
            print('\t%s:' % nc_attr, repr(nc_fid.getncattr(nc_attr)))
    nc_dims = [dim for dim in nc_fid.dimensions]  # list of nc dimensions
    # Dimension shape information.
    if verb:
        print("NetCDF dimension information:")
        for dim in nc_dims:
            print("\tName:", dim)
            print("\t\tsize:", len(nc_fid.dimensions[dim]))
            print_ncattr(dim)
    # Variable information.
    nc_vars = [var for var in nc_fid.variables]  # list of nc variables
    if verb:
        print("NetCDF variable information:")
        for var in nc_vars:
            if var not in nc_dims:
                print('\tName:', var)
                print("\t\tdimensions:", nc_fid.variables[var].dimensions)
                print("\t\tsize:", nc_fid.variables[var].size)
                print_ncattr(var)
    return nc_attrs, nc_dims, nc_vars

nc_attrs, nc_dims, nc_vars = ncdump(nc_fdata)

我猜你需要名为“ua”的变量，因为它同时包含经度和纬度地址。

因此，为了构造遮罩，您必须从“ua”中提取所有内容，其中经度和纬度位于 shapefile 的边界框值之间。

【讨论】：

这如何回答这个问题？阅读 NetCDF 和 shape 文件是问题的一部分（几乎每个在线教程都涉及）。

在第一行中，我写道这不是一个完整的解决方案。用户想要从文件中提取数据，我帮助他如何在 python 中做到这一点，并提示他如何继续。

我的意见：你把最困难的部分留了下来（只有 "solve" 琐碎的部分），在这种情况下我不会称之为答案。跨度>

这并不是为了解决问题，我只是想给用户一个先机。另外，评论太多了，所以我不得不写成“答案”。