0

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 ?

1 Answers1

0

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