Search code examples
projectionmatplotlib-basemapgrib

Issue with the robin projection and contourf data with Basemap


I'm using the basemap library to display spatial information from Copernicus program. The issue is i can not figure out how to project the data on the robin projection, but I do it correctly with the orthogonal projection.

So currently, I tried this :

plt.ioff()

    # adapt for location of datasources
    filePath = '../data/grib/download.grib'

    # load data
    grbs = grb.open(filePath)
    grbs.seek(0)

    data, lats, lons = (None, None, None)
    dataUnit = None
    title = None
    for g in grbs:
        data, lats, lons = g.data()

        name = g.name
        level = g.level
        pressureUnit = g.pressureUnits
        date = g.validDate

        dataUnit = g.units

        title = name + ' at ' + str(level) + ' ' + str(pressureUnit) + ' [' + str(date) + ']'
        print(title)

        break

#     mapPlot = Basemap(projection='ortho', lat_0=0, lon_0=0)
    mapPlot = Basemap(projection='robin', lat_0=0, lon_0=0, resolution='l')
    mapPlot.drawcoastlines(linewidth=0.25)

    x, y = mapPlot(lons, lats)
    mapPlot.contourf(x, y, data)
    mapPlot.colorbar(location='bottom', format='%.1f', label=dataUnit)

    plt.title(title)
    plt.show()

The orthogonal projection works correctly. But for the robin projection, I have an ... interesting pattern.

What I'm doing wrong ?


Solution

  • So i figure out how to do. I was misled but the first examples I saw.

    Here is a my code:

    
            import matplotlib
            from mpl_toolkits.basemap import Basemap, shiftgrid
    
            import matplotlib.pyplot as plt
            import numpy as np
            import pygrib as grb
    
            # Get data
            data = g['values']
            lats = g['distinctLatitudes'] # 1D vector
            lons = g['distinctLongitudes'] # 1D vector
    
            # Useful information for late
            name = g.name
            level = str(g.level) + g.pressureUnits
            date = g.validDate
            dataUnit = g.units
    
            # Parse the data
            # Shit the data to start à -180. This is important to mark the data to start at -180°
            data, lons = shiftgrid(180., data, lons, start=False)  # shiftgrid
    
            # Choose a representation (works with both)
    #         mapPlot = Basemap(projection='ortho', lat_0=0, lon_0=0)
            mapPlot = Basemap(projection='robin', lat_0=0, lon_0=0)
            mapPlot.drawcoastlines(linewidth=0.25)
    
            # Convert the coordinates into the map projection
            x, y = mapPlot(*np.meshgrid(lons, lats))
    
            # Display data
            map = mapPlot.contourf(x, y, data, levels=boundaries, cmap=plt.get_cmap('coolwarm'))
    
            # Add what ever you want to your map.
            mapPlot.nightshade(date, alpha=0.1)
    
            # Legend
            mapPlot.colorbar(map, label=dataUnit)
    
            # Title
            plt.title(name + ' at ' + str(level) + ' [' + str(date) + ']')
    
            plt.show()
    
    

    So it returns what I'm expecting. Result