3

海氷データのプロットを作成しようとしています。データは EASE-North グリッドで配信されます。サンプル ファイル (HDF4) は次の場所からダウンロードできます。

ftp://n4ftl01u.ecs.nasa.gov/SAN/OTHR/NISE.004/2013.09.30/

EASE-Grid 用のカスタム プロジェクション クラスを作成しましたが、機能しているようです (海岸線はデータとうまく一致しています)。

Natural Earth 機能を追加しようとすると、空の Matplotlib フィギュアが返されます。

import gdal
import cartopy

# projection class
class EASE_North(cartopy.crs.Projection):
    
    def __init__(self):
        
        # see: http://www.spatialreference.org/ref/epsg/3408/
        proj4_params = {'proj': 'laea',
            'lat_0': 90.,
            'lon_0': 0,
            'x_0': 0,
            'y_0': 0,
            'a': 6371228,
            'b': 6371228,
            'units': 'm',
            'no_defs': ''}
        
        super(EASE_North, self).__init__(proj4_params)
        
    @property
    def boundary(self):
        coords = ((self.x_limits[0], self.y_limits[0]),(self.x_limits[1], self.y_limits[0]),
                  (self.x_limits[1], self.y_limits[1]),(self.x_limits[0], self.y_limits[1]),
                  (self.x_limits[0], self.y_limits[0]))
        
        return cartopy.crs.sgeom.Polygon(coords).exterior
        
    @property
    def threshold(self):
        return 1e5
    
    @property
    def x_limits(self):
        return (-9000000, 9000000)

    @property
    def y_limits(self):
        return (-9000000, 9000000)


# read the data
ds = gdal.Open('D:/NISE_SSMISF17_20130930.HDFEOS')

# this loads the layers for both hemispheres
data = np.array([gdal.Open(name, gdal.GA_ReadOnly).ReadAsArray()
                 for name, descr in ds.GetSubDatasets() if 'Extent' in name])

ds = None

# mask anything other then sea ice
sea_ice_concentration = np.ma.masked_where((data < 1) | (data > 100), data, 0)

# plot
lim = 3000000

fig, ax = plt.subplots(figsize=(8,8),subplot_kw={'projection': EASE_North(), 'xlim': [-lim,lim], 'ylim': [-lim,lim]})

land = cartopy.feature.NaturalEarthFeature(
            category='physical',
            name='land',
            scale='50m',
            facecolor='#dddddd',
            edgecolor='none')

#ax.add_feature(land)
ax.coastlines()

# from the metadata in the HDF
extent = [-9036842.762500, 9036842.762500, -9036842.762500, 9036842.762500]

ax.imshow(sea_ice_concentration[0,:,:], cmap=plt.cm.Blues, vmin=1,vmax=100,
          interpolation='none', origin='upper', extent=extent, transform=EASE_North())

上記のスクリプトは正常に動作し、次の結果が生成されます。

ここに画像の説明を入力

しかし、コメントを外すと、ax.add_feature(land)エラーなしで失敗し、空の図のみが返されます。明らかな何かが欠けていますか？

IPython ノートブックは次のとおりです: http://nbviewer.ipython.org/6779935

私の Cartopy ビルドは、Christoph Gohlke の Web サイト (ありがとう!) のバージョン 0.9 です。

編集：

Figure を保存しようとすると、例外がスローされます。

fig.savefig(r'D:\test.png')

C:\Python27\Lib\site-packages\shapely\speedups\_speedups.pyd in shapely.speedups._speedups.geos_linearring_from_py (shapely/speedups/_speedups.c:2270)()

ValueError: A LinearRing must have at least 3 coordinate tuples

「土地」を調べてcartopy.featureも問題はありません。すべてのポリゴンが通過し.isvalid()、すべてのリング (ext en int) が 4 つ以上のタプルで構成されています。したがって、入力形状は問題ではないようです (そして、PlateCaree() で正常に動作します)。

EASE_North に変換した後、一部のリング (南半球など) が「破損」する可能性がありますか?

編集2:

組み込みの NE 機能を削除して同じシェープファイルをロードすると (ただし、40N 未満のものはすべてクリップされます)、機能します。したがって、ある種の再投影の問題のようです。

for state in shpreader.Reader(r'D:\ne_50m_land_clipped.shp').geometries():
    ax.add_geometries([state], cartopy.crs.PlateCarree(),facecolor='#cccccc', edgecolor='#cccccc')

ここに画像の説明を入力

4

1 に答える 1

1

これはバグだと言っていたでしょう。私はadd_featurematplotlib viewLimを更新すると推測しています。その結果、画像が小さな領域にズームインします(ズームアウトしない限り、白く表示されます)。

私の頭の上から、根本的な動作は matplotlib で改善されたと思いますが、cartopy はまだ新しい viewLim 計算を利用していません。当面は、マップの範囲を次のように手動で設定することをお勧めします。

ax.set_extent(extent, transform=EASE_North())

HTH

于 2013-10-03T09:11:48.130 に答える