1

I have successfully opened a grib2 file from 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.

I got what I expect, but only for half of the world, I can solve it by subtracting 180.0 from my xmin and xmax, but then I lose the coordinate conversion, I guess the problem is that I am not shifting the data, possibly using shiftgrid from mpl_toolkits, but I can not get the function to work either, any suggestions?

Here is an image of the map without the subtraction:

enter image description here

Here is what I get when I subtract 180.0 from the xmin and xmax variables:

Bad transformation of coordinates

You can download the grib2 file I am using from: 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()
0

Here is what I came with that works with all projections:

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()

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.