使用底图时更改地图投影

Question

我之前更改了底图中的地图投影，所以我知道这应该是一个简单的解决方法，但似乎没有什么对我有用。我使用了网格网格，并映射了我的 x、y 值等，但我只是得到了扭曲或疯狂的图。我认为这与我正在使用的数据被自动设置为在兰伯特保形（我不想要）上绘制的事实有关，而且尺寸的单位是公里而不是纬度和经度。我不知道从这里去哪里......

数据来源：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 值）。在我们尝试其他任何事情之前，必须先解决这个问题。
• 我不确定您为什么要通过除以 1000 来转换 x 和 y 数据（LCC 中的投影坐标？）。

无论如何，一个明智的做法如下：

下载您的数据，并使用 Pyproj 将坐标从 LCC 转换为 lon/lat：

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 进行绘图（不要使用 jet 颜色图！）

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