更新: いくつかのデータをマッピングしようとしています。グリッド内の基準点から測定された一連の後方方位角 (baz) があります。バズに沿った大円が交差するグリッド上のすべての点を見つけたいと思います。これを行うには、グリッド内の各ポイントを繰り返し処理し、そのポイントと参照ポイントの間の予想される後方方位角を計算し、測定された各 baz と比較します。両者の差が小さい (2 度未満) 場合は、その点に重みを付けます。次に、すべてを地図に載せました。私が使用するコードは以下のとおりですが、結果は少し奇妙に見えます。どこが間違っているのか誰か知っていますか?
from matplotlib.colorbar import ColorbarBase
import matplotlib.cm as cm
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import mpl_toolkits.basemap.pyproj as pyproj
llcrnrlon = -30.0
llcrnrlat = 45.0
urcrnrlon = 0.0
urcrnrlat = 65.0
lon_0 = (urcrnrlon + llcrnrlon) / 2.
lat_0 = (urcrnrlat + llcrnrlat) / 2.
lat = 51.58661577 # reference point
lon = -9.18822525
# Generate random back-azimuths.
baz = zeros((20))
for i in xrange(len(baz)):
baz[i] = random.randint(200,230)
####################################################################
## Set up the map background.
m = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,
resolution='i',projection='lcc',lon_0=lon_0,lat_0=lat_0)
m.drawcoastlines()
m.fillcontinents()
# draw parallels
m.drawparallels(np.arange(10,70,10),labels=[1,0,0,0])
# draw meridians
m.drawmeridians(np.arange(-80, 25, 10),labels=[0,0,0,1])
# Plot station locations.
x, y = m(lon, lat) # array ref points
m.plot(x,y,'ro', ms=5)
####################################################################
## Set up the grids etc.
glons = np.linspace(llcrnrlon, urcrnrlon, 100)
glats = np.linspace(llcrnrlat, urcrnrlat, 100)
# Convert to map coords.
xlons, ylats = m(glons, glats)
# create grid for pcolormesh.
grid_lon, grid_lat = np.meshgrid(xlons, ylats)
# create weights for pcolormesh.
weights = np.zeros(np.shape(grid_lon))
# create grid of lat-lon coords for baz calculation.
gln, glt = np.meshgrid(glons, glats)
####################################################################
## calculate baz from grid_lon, grid_lat to lon, lat. If less
## than error weight grid point.
# method for BAZ calculation via pyproj.
def get_baz(lon1, lat1, lon2, lat2):
g = pyproj.Geod(ellps='WGS84')
az, baz, dist = g.inv(lon1, lat1, lon2, lat2)
return baz
# BAZ calcultion for each point in grid.
ll=0
for mBAZ in baz:
for i in xrange(len(gln)):
for k in xrange(len(gln[i])):
nbaz = get_baz(lon, lat, gln[i][k], glt[i][k])
nbaz += 180
if abs(nbaz - mBAZ) < 2:
weights[i][k] = 1
ll+=1
# plot grid.
m.pcolormesh(grid_lon, grid_lat, weights, cmap=plt.cm.YlOrBr)
plt.colorbar()
plt.show()
以下の元の質問は、現在では古くなっています。
いくつかのデータをマッピングしようとしています。各方向の値 (周波数) の範囲を示すデータセットがあります。特定の方位角に沿った各グリッド ポイントが特定の周波数の電力によって重み付けされるように、それらをグリッド上にプロットしたいと考えています。次のように、ベースマップを使用してマップを作成し、その上にグリッドをプロットしました。
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
from shoot import *
llcrnrlon = -20.0
llcrnrlat = 45.0
urcrnrlon = 10.0
urcrnrlat = 65.0
lon_0 = (urcrnrlon + llcrnrlon) / 2.
lat_0 = (urcrnrlat + llcrnrlat) / 2.
m = Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,
resolution='i',projection='lcc',lon_0=lon_0,lat_0=lat_0)
## Set up the grid.
glons = np.linspace(-20,10,50)
glats = np.linspace(45, 65, 50)
xlons, ylats = m(glons, glats)
grid_lon, grid_lat = np.meshgrid(xlons, ylats)
pwr = np.zeros((50,50))
m.drawcoastlines()
m.fillcontinents()
# draw parallels
m.drawparallels(np.arange(10,70,10),labels=[1,0,0,0])
# draw meridians
m.drawmeridians(np.arange(-80, 25, 10),labels=[0,0,0,1])
lats = [54.8639587, 51.5641564]
lons = [-8.1778180, -9.2754284]
x, y = m(lons, lats) # array ref points
# Plot station locations.
m.plot(x,y,'ro', ms=5)
m.pcolormesh(grid_lon, grid_lat, pwr)
次に、この素敵なサイトで見つけたいくつかの機能を使用して、必要な大円を撃ちます
glon1 = lons[0]
glat1 = lats[0]
azimuth = 280.
maxdist = 200.
great(m, glon1, glat1, azimuth, color='orange', lw=2.0)
plt.show()
ただし、線をプロットするだけでは不十分です。大円が交差するグリッド ポイントを見つけて、それらに値を割り当てることができるようにしたいと考えています。誰もこれについて行く方法を知っていますか??
