从 GRIB2 文件移动数据

Question

我已经从 NCEP and I am having trouble being able to convert the coordinates to plot them using matplotlib, using the custom convertXY function from this post Plot GDAL raster using matplotlib Basemap 成功打开了一个 grib2 文件。

我得到了我所期望的，但只有一半的世界，我可以通过从我的 xmin 和 xmax 中减去 180.0 来解决它，但我想我失去了坐标转换问题是我没有移动数据，可能使用 mpl_toolkits 中的 shiftgrid，但我也无法使该功能正常工作，有什么建议吗？

这是一张没有减法的地图图片：

这是我从 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()

Answer 1

这是我带来的适用于所有投影的内容：

from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid
import osr, gdal
import matplotlib.pyplot as plt
import numpy as np

# Read the data and metadata
# Pluviocidad.
#ds = gdal.Open( 'C:\Users\Paula\Downloads\enspost.t00z.prcp_24hbc (1).grib2', gdal.GA_ReadOnly )
# Sea Ice
ds = gdal.Open( 'D:\Downloads\seaice.t00z.grb.grib2', gdal.GA_ReadOnly )
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()

xres = gt[1]
yres = gt[5]

xsize = ds.RasterXSize
ysize = ds.RasterYSize

# get the edge coordinates and add half the resolution 
# to go to center coordinates
xmin = gt[0] + xres * 0.5
xmax = gt[0] + (xres * xsize) - xres * 0.5
ymin = gt[3] + (yres * ysize) + yres * 0.5
ymax = gt[3] - yres * 0.5

ds = None

xx = np.arange( xmin, xmax + xres, xres )
yy = np.arange( ymax + yres, ymin, yres )

data, xx = shiftgrid( 180.0, data, xx, start = False )

# Mercator
m = Basemap(projection='merc',llcrnrlat=-85,urcrnrlat=85,\
            llcrnrlon=-180,urcrnrlon=180,lat_ts=0,resolution='c')

x, y = m(*np.meshgrid(xx,yy))

# plot the data (first layer) data[0,:,:].T
im1 = m.pcolormesh( x, y, data, shading = "flat", cmap=plt.cm.jet )

# annotate
m.drawcountries()
m.drawcoastlines(linewidth=.5)

plt.show()

从 GRIB2 文件移动数据

Shifting data from a GRIB2 file

python

gis

numpy

gdal

matplotlib-basemap