2017-11-16 9 views
1

cartopy、matplotlib、およびimshowを使用して、等間隔(緯度/経度)データの正方形グリッドをプロットしようとしています。データはデータラインを横切っており、正しく動作するようにマップを取得することに問題がありました。imshowとcartopyのデータ交差点を含む画像で

ここに私の問題の例です:

import numpy as np 
import cartopy.crs as ccrs 
import matplotlib.pyplot as plt 


lat = np.arange(6000)*0.02 + (-59.99) 
lon = np.arange(6000)*0.02 + (85.01) 

dat = np.reshape(np.arange(6000*6000),[6000,6000]) 

tran = ccrs.PlateCarree() 
proj = tran 

plt.figure(figsize=(8,8)) 

ax = plt.axes(projection=proj) 
print([lon[0],lon[-1],lat[0],lat[-1]]) 
ax.imshow(dat, extent=[lon[0],lon[-1],lat[0],lat[-1]],transform=tran,interpolation='nearest') 
ax.coastlines(resolution='50m', color='black', linewidth=2) 
ax.gridlines(crs=proj,draw_labels=True) 
plt.show() 

tran = ccrs.PlateCarree(central_longitude=180) 
proj = tran 

plt.figure(figsize=(8,8)) 

ax = plt.axes(projection=proj) 
print([lon[0]-180,lon[-1]-180,lat[0],lat[-1]]) 
ax.imshow(dat, extent=[lon[0]-180,lon[-1]-180,lat[0],lat[-1]],transform=tran,interpolation='nearest') 
ax.coastlines(resolution='50m', color='black', linewidth=2) 
ax.gridlines(crs=tran,draw_labels=True) 
plt.show() 

は、最初のプロットは180Eで降りチョッピング、この画像が得られます。 chopped at 180

二フィックスマップの問題を、しかし、グリッドのダニは今間違っています: 180 now zero

私は再投射を試みましたが、私は(tran!= proj)と思っていますが、見かけ上ハングアップしすぎていました。

私は基本的には下の画像が必要ですが、適切なラベルが付いています。私は、より多くの地理的に離れたデータを重ねて表示するつもりだから、今はちょうどハックのように見えるのではなく、正確にそれをやりたいのです。

+0

はcartopyを持っていない、しかし、あなたはちょうどあなたが 'ax.set_xticklabels'で欲しいものにラベルを変更することはできませんか? –

+0

私の答えに関するコメント/質問はありますか? – swatchai

答えて

0

Cartopyを使用すると、地図を横切るデータ線を描画することは常に困難です。ここにあなたが望む地図をプロットするコードがあります。

import numpy as np 
import cartopy.crs as ccrs 
import matplotlib.pyplot as plt 

# demo raster data 
n1 = 300 
m1 = 0.4 
lat = np.arange(n1)*m1 + (-59.99) 
lon = np.arange(n1)*m1 + (85.01) 
dat = np.reshape(np.arange(n1*n1), [n1,n1]) 

cm_lon=180 # for central meridian 

tran = ccrs.PlateCarree(central_longitude = cm_lon) 
proj = tran 

plt.figure(figsize=(8,8)) 
ax = plt.axes(projection=proj) 

ext = [lon[0]-cm_lon, lon[-1]-cm_lon, lat[0], lat[-1]] 
#print(ext) 

ax.imshow(dat, extent=ext, \ 
       transform=tran, interpolation='nearest') 

ax.coastlines(resolution='110m', color='black', linewidth=0.5, zorder=10) 

# this draws grid lines only, must go beyond E/W extents 
ax.gridlines(draw_labels=False, xlocs=[80,100,120,140,160,180,-180,-160,-140]) 

# this draw lables only, exclude those outside E/W extents 
ax.gridlines(draw_labels=True, xlocs=[100,120,140,160,180,-160]) 

plt.show() 

結果のマップ:

output