从netCDF文件中提取数据包含在shape文件的边界内

AAAAA:

我有以下shape文件netCDF文件

我想从所包含shape文件的边界内的netCDF文件中提取数据。

你有我如何能做到这一点的任何建议?

shape文件对应于SREX区域11 北欧(NEU)和netCDF文件是CMIP6气候模式数据输出(UA变量)的一个例子。我需要的输出必须是在创建NetCDF格式。


更新

到目前为止,我试图创建使用NCL和CDO一个的NetCDF面具,这个面具向原的NetCDF数据集。这里下面的步骤(和NCL脚本):

#################
## 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 
#################

输出几乎是正确的。见下图。但是shape文件的南部边界(SREX 11,NEU)不正确捕获。所以我想,有什么问题在生成的NetCDF面具NCL脚本。

在这里输入图像描述

巴特:

再使用一些旧脚本/代码,我很快就想出了这个为一个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()

本文收集自互联网,转载请注明来源。

如有侵权,请联系 [email protected] 删除。

编辑于
0

我来说两句

0 条评论
登录 后参与评论

相关文章

来自分类Dev

有条件地将元数据包含在Android清单文件中

来自分类Dev

Python从包含多个地标的KML文件中提取数据

来自分类Dev

从.wav文件中提取数据

来自分类Dev

使用路径从对象提取数据-将所有数据包含在数组中

来自分类Dev

NCO:使用NCO ncks从NetCDF文件中提取变量

来自分类Java

写JSON文件从网页中提取的数据

来自分类Python

如何从python中的文件中提取数据

来自分类Dev

从tesseract hocr xhtml文件中提取数据

来自分类Dev

从python文件中提取并绘制数据

来自分类Dev

标头和边界包含在上载的文件Django Rest框架中

来自分类Linux

从文件中提取值

来自分类Java

从包含文件名的字符串中提取模式

来自分类Dev

如何从gulp包含的回调中提取文件名?

来自分类Dev

如何从文件中提取数据以用作C中的变量

来自分类Dev

使用python从xml文件中提取数据时出错

来自分类Dev

如何在python中提取文件数据

来自分类Dev

无法从工厂AngularJS中的JSON文件中提取数据

来自分类Dev

从文本文件中提取模式之间的数据

来自分类Dev

Golang:从1个CSV文件中提取数据到anthoer

来自分类Dev

从Access数据库的“附件”字段中提取文件

来自分类Dev

使用python从JSON文件中提取部分数据

来自分类Dev

用PHP从.osm XML文件中提取数据

来自分类Dev

从文本文件中提取特定数据

来自分类Dev

从下载的文件中提取源元数据

来自分类Dev

U-SQL-从复杂的嵌套json文件中提取数据

来自分类Dev

如何检测3D对象何时包含在由4个坐标定义的2D边界内?

来自分类Dev

从CSV文件中提取包含来自另一个CSV文件的任何值的行

来自分类Linux

从文件中提取特定列?

来自分类Dev

从单个文件中提取值

TOP 榜单

热门标签

归档