使用底圖時更改地圖投影

Question

我之前更改過底圖中的地圖投影，所以我知道這應該是一個簡單的修復方法，但似乎沒有任何適用於我的方法。我已經使用了meshgrid，並且映射了我的x，y值等等，而我只是被扭曲或者瘋狂的繪製。我認爲這與我正在使用的數據自動設置爲繪製在Lambert Conformal（我不想要）上的事實有關，也是以km爲單位而非緯度和經度。我不知道在哪裏何去何從......

數據來源：http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/RAP/CONUS_13km/RR_CONUS_13km_20150415_0600.grib2/GC.html

這是我工作的代碼。我有一堆評論的東西，我一直在嘗試沒有運氣。

所有的

import numpy as np 
import math as m 
import urllib2 
import time 
import datetime as dt 
import matplotlib.pyplot as plt 
from mpl_toolkits.axes_grid1 import make_axes_locatable 
from mpl_toolkits.basemap import Basemap, shiftgrid 
from matplotlib.colors import LinearSegmentedColormap 
from pydap.client import open_url 
from pydap.proxy import ArrayProxy 
import scipy 

data_url = 'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/RAP/CONUS_13km/RR_CONUS_13km_' + '20150415' + '_' + '01' + '00.grib2/GC' 

print('Getting Data from URL:\n\n "{0}"\n'.format(data_url)) 

# Create Array of all data from URL 
dataset = open_url(data_url) 

# Map Projection Info 

proj_attributes = dataset['LambertConformal_Projection'].attributes 
rsphere = proj_attributes['earth_radius'] 

lat_0 = proj_attributes['latitude_of_projection_origin'] 
lon_0 = proj_attributes['longitude_of_central_meridian'] 
lat_1 = proj_attributes['standard_parallel'] 

llcrnrlat = 16.28100013732909 # (1,1) 
llcrnrlon = 360-126.18 # (1,1) 

urcrnrlat = 55.552133975329625 # (614,428) 
urcrnrlon = 360-59.15590040502627 # (614,248) 

x = np.array(dataset['x'][:]) 
y = np.array(dataset['y'][:]) 

def xy_converter(var): 
    """ 
    Downloads entered variable (x or y) coordinates 
    and converts from m to km. Inputs for var 
    should be 'x' or 'y'. 
    """ 
    values = dataset[var][:] 
    data_array = values * 1000 
    newarray = data_array + abs(data_array.min()) 
    return newarray 

# Download x & y coord. and convert m to km 
x = xy_converter('x') 
y = xy_converter('y') 

# Temp Contour 
temp_2m = dataset['Temperature_height_above_ground'].array[1,:,:,:]-273. 
temp_2m = temp_2m * (9./5.) + 32. 
temp_2m = temp_2m.squeeze() 

#plot 
fig = plt.figure(figsize=(11,11)) 
ax = fig.add_subplot(1,1,1) 

map = Basemap(projection='lcc', lat_0 = lat_0, lon_0 = lon_0, 
            llcrnrlon = llcrnrlon, llcrnrlat = llcrnrlat, 
            urcrnrlat = urcrnrlat, urcrnrlon = urcrnrlon, 
            area_thresh = 1000., rsphere = rsphere, resolution='i') 
map.drawcoastlines(linewidth=0.3) 
map.drawcountries(linewidth=0.3) 
map.drawcounties(linewidth=0.1) 
map.drawstates(linewidth=0.3) 
map.drawmapboundary(linewidth=0.5) 

#lons,lats = basemap_parameters.map(x,y) 
#lon,lat = basemap_parameters.map(lons,lats,inverse=True) 

#ny = range(len(y)); nx = range(len(x)) 
#ny = temp_2m.shape[0]; nx = temp_2m.shape[1] 
#lons, lats = map(x, y) # get lat/lons of ny by nx evenly space grid. 
#xx, yy = map(lons, lats) 

levels = np.linspace(-42,122,320) 
ticks = [-60,-50,-40,-30,-20,-10,0,10,20,30,40,50,60,70,80,90,100,110,120] 

plot = plt.contourf(x,y,temp_2m,levels,cmap='jet',extend='both') 

# Set Colorbar Text Color 
color_bar = map.colorbar(plot) 
color_bar.set_ticks(ticks) 

# CONUS 
plt.xlim(x[30],x[440]) 
plt.ylim(y[30],y[290])  
plt.savefig('/home/public_html/conus_temp.png', dpi=100, bbox_inches='tight', pad_inches = .05) 

Answer 1

• 首先，你不應該定義底圖實例作爲map - 它的陰影Python的內置map()功能，只會造成混亂。
• 其次，您下載的x和y值是不同的長度（分別爲451和337值）。在我們嘗試其他任何事情之前，必須先解決這個問題。
• 我不知道你爲什麼除以轉換您的X和Y數據（？在LCC投影座標）由1000

在任何情況下，明智的做法如下：

從LCC下載您的數據，並將其轉換座標經度/緯度使用Pyproj：

import pyproj 
lc = pyproj.Proj("+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs") 
lons, lats = lc(x, y, inverse=True) 
# lons & lats are now unprojected WGS84 

轉換你的座標，你所需地圖投影座標：

現在

# now you can get ll and ur lons and lats 
# set up your basemap with whatever projection you'd like, e.g. 
m = Basemap(
    projection = 'merc', 
    llcrnrlon = lllon, llcrnrlat = lllat, urcrnrlon = urlon, urcrnrlat = urlat, 
    resolution='h') 
# project lons & lats into map coordinates 
projected_lons, projected_lats = m(x, y) 

我們可以做其他的東西：

我不知道，如果你的溫度數據已插入到網格。如果不是：

from matplotlib.mlab import griddata 
# set up a square grid with the same extents as our measured data 
numcols, numrows = 1000, 1000 
xi = np.linspace(projected_lons.min(), projected_lons.max(), numcols) 
yi = np.linspace(projected_lats.min(), projected_lats.max(), numrows) 
# get lon and lat coords of our grid points 
xi, yi = np.meshgrid(xi, yi) 
# interpolate 
x, y, z = projected_lons, projected_lats, temp2m 
zi = griddata(x, y, z, xi, yi) 

設置你的其他地圖細節往常一樣，用contourf圖（不使用噴氣顏色表！）

conf = m.contourf(xi, yi, zi, cmap='coolwarm', ax=ax)