5

いくつかのポリゴン(カナダの州) があり、 を使用して読み取り、GeoPandasこれらを使用してマスクを作成し、2 次元の緯度経度グリッド上のグリッド データに適用したいと考えています ( を使用してnetcdfファイルから読み取りますiris)。最終的な目標は、特定の州のデータのみを残し、残りのデータをマスクすることです。したがって、マスクは、州内のグリッド ボックスの場合は 1 になり、州外のグリッド ボックスの場合は 0 または NaN になります。

ポリゴンは、こちらのシェープファイルから取得できます: https://www.dropbox.com/s/o5elu01fetwnobx/CAN_adm1.shp?dl=0

私が使用している netcdf ファイルは、 https ://www.dropbox.com/s/kxb2v2rq17m7lp7/t2m.20090815.nc?dl=0 からダウンロードできます。

ここには2つのアプローチがあると思いますが、両方に苦労しています：

1) 多角形を使用して緯度経度グリッドにマスクを作成し、これを Python 以外の多くのデータファイルに適用できるようにします (推奨)。

2) 多角形を使用して、読み込まれたデータをマスクし、関心のある州内のデータのみを抽出して、対話的に操作します。

これまでの私のコード:

import iris
import geopandas as gpd

#read the shapefile and extract the polygon for a single province
#(province names stored as variable 'NAME_1')
Canada=gpd.read_file('CAN_adm1.shp')
BritishColumbia=Canada[Canada['NAME_1'] == 'British Columbia']

#get the latitude-longitude grid from netcdf file
cubelist=iris.load('t2m.20090815.nc')
cube=cubelist[0]
lats=cube.coord('latitude').points
lons=cube.coord('longitude').points

#create 2d grid from lats and lons (may not be necessary?)
[lon2d,lat2d]=np.meshgrid(lons,lats)

#HELP!

助けやアドバイスをありがとう。

更新: 以下の @DPeterK の優れたソリューションに従って、元のデータをマスクして、次のようにすることができます。

カナダの州のシェープファイルを使用してマスクされた netcdf ファイルからの温度データ

4

2 に答える 2

8

順調にスタートできたようです！シェープファイルから読み込まれたジオメトリは、さまざまな地理空間比較方法を公開します。この場合、そのcontains方法が必要です。これを使用して、立方体の水平グリッド内の各ポイントがブリティッシュ コロンビア州のジオメトリ内に含まれているかどうかをテストできます。(これは高速な操作ではないことに注意してください!) この比較を使用して 2D マスク配列を作成し、キューブのデータに適用したり、他の方法で使用したりできます。

上記を実行するための Python 関数を作成しました。この関数は、立方体とジオメトリを受け取り、立方体の (指定された) 水平座標のマスクを生成し、そのマスクを立方体のデータに適用します。機能は以下です。

def geom_to_masked_cube(cube, geometry, x_coord, y_coord,
                        mask_excludes=False):
    """
    Convert a shapefile geometry into a mask for a cube's data.

    Args:

    * cube:
        The cube to mask.
    * geometry:
        A geometry from a shapefile to define a mask.
    * x_coord: (str or coord)
        A reference to a coord describing the cube's x-axis.
    * y_coord: (str or coord)
        A reference to a coord describing the cube's y-axis.

    Kwargs:

    * mask_excludes: (bool, default False)
        If False, the mask will exclude the area of the geometry from the
        cube's data. If True, the mask will include *only* the area of the
        geometry in the cube's data.

    .. note::
        This function does *not* preserve lazy cube data.

    """
    # Get horizontal coords for masking purposes.
    lats = cube.coord(y_coord).points
    lons = cube.coord(x_coord).points
    lon2d, lat2d = np.meshgrid(lons,lats)

    # Reshape to 1D for easier iteration.
    lon2 = lon2d.reshape(-1)
    lat2 = lat2d.reshape(-1)

    mask = []
    # Iterate through all horizontal points in cube, and
    # check for containment within the specified geometry.
    for lat, lon in zip(lat2, lon2):
        this_point = gpd.geoseries.Point(lon, lat)
        res = geometry.contains(this_point)
        mask.append(res.values[0])

    mask = np.array(mask).reshape(lon2d.shape)
    if mask_excludes:
        # Invert the mask if we want to include the geometry's area.
        mask = ~mask
    # Make sure the mask is the same shape as the cube.
    dim_map = (cube.coord_dims(y_coord)[0],
               cube.coord_dims(x_coord)[0])
    cube_mask = iris.util.broadcast_to_shape(mask, cube.shape, dim_map)

    # Apply the mask to the cube's data.
    data = cube.data
    masked_data = np.ma.masked_array(data, cube_mask)
    cube.data = masked_data
    return cube

2D マスクだけが必要な場合は、上記の関数がキューブに適用される前にそれを返すことができます。

元のコードでこの関数を使用するには、コードの最後に次を追加します。

geometry = BritishColumbia.geometry
masked_cube = geom_to_masked_cube(cube, geometry,
                                  'longitude', 'latitude',
                                  mask_excludes=True)

これで何もマスクされない場合は、キューブとジオメトリが異なる範囲で定義されている可能性があります。つまり、キューブの経度座標が 0° ～ 360° の範囲であり、ジオメトリの経度値が -180° ～ 180° の範囲である場合、封じ込めテストは を返しませんTrue。これは、キューブの範囲を次のように変更することで修正できます。

cube = cube.intersection(longitude=(-180, 180))
于 2017-12-14T12:17:11.850 に答える