我已经成功地从NCEP打开了一个grib2文件,并且我很难使用matplotlib
转换坐标来绘制它们,使用本文Plot GDAL raster using matplotlib Basemap中的自定义convertXY
函数。你知道吗
我得到了我所期望的,但是只有一半的人,我可以通过从我的xmin
和xmax
中减去180.0来解决它,但是我失去了坐标转换,我猜问题是我没有移动数据,可能使用shiftgrid
从mpl_toolkits
,但是我也不能让函数工作,有什么建议吗?你知道吗
这是一张没有减法的地图图像:
下面是我从xmin
和xmax
变量中减去180.0得到的结果:
您可以从以下网址下载我正在使用的grib2文件: https://drive.google.com/open?id=1RsuiznRMbJNpNsrQeXEunvVsWZJ0tL2d
from mpl_toolkits.basemap import Basemap
import osr, gdal
import matplotlib.pyplot as plt
import numpy as np
def convertXY(xy_source, inproj, outproj):
# function to convert coordinates
shape = xy_source[0,:,:].shape
size = xy_source[0,:,:].size
# the ct object takes and returns pairs of x,y, not 2d grids
# so the the grid needs to be reshaped (flattened) and back.
ct = osr.CoordinateTransformation(inproj, outproj)
xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T))
xx = xy_target[:,0].reshape(shape)
yy = xy_target[:,1].reshape(shape)
return xx, yy
# Read the data and metadata
ds = gdal.Open(r'D:\Downloads\flxf2018101912.01.2018101912.grb2')
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
xres = gt[1]
yres = gt[5]
# get the edge coordinates and add half the resolution
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmin -= 180.0
xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5
xmax -= 180.0
ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5
ymax = gt[3] - yres * 0.5
ds = None
# create a grid of xy coordinates in the original projection
xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
# Create the figure and basemap object
fig = plt.figure(figsize=(12, 6))
m = Basemap(projection='robin', lon_0=0, resolution='c')
# Create the projection objects for the convertion
# original (Albers)
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)
# Get the target projection from the basemap object
outproj = osr.SpatialReference()
outproj.ImportFromProj4(m.proj4string)
# Convert from source projection to basemap projection
xx, yy = convertXY(xy_source, inproj, outproj)
# plot the data (first layer)
im1 = m.pcolormesh(xx, yy, data[0,:,:].T, cmap=plt.cm.jet)
# annotate
m.drawcountries()
m.drawcoastlines(linewidth=.5)
plt.show()
以下是我提供的适用于所有投影的内容:
相关问题 更多 >
编程相关推荐