Search code examples
arraysprojectioncartopy

Generating multiple maps from 2d array from netcdf


I need to generate multiple temperature plots for every month of every year, spanning from 2002 to 2018. I have managed to develop one plot for all of the data in 2002 (about 6 months). I had import my netcdf files into an array and slice the 3d array to a 2d array of lat lon lst. Using pcolormesh I plotted one singular lon, lat, data_2d plot but can't figure out how to define the months so I can plot them separately. The data is not available online but I am looking for a general function or command that can iterate through the months and plot them separately onto a map.

  import netCDF4 
    import numpy as np                             
    import xarray as xr
    import pandas as pd
    import os
    from netCDF4 import Dataset
    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    print ('All packages imported')
    # make an empty np array
    data_grid = []
    # list of files 
    days = pd.date_range (start='4/7/2002', end='31/7/2002')
    print ('Initiate for loop')

    # for loop iterating through %date range to find path files
    for day in days:
        # defining what paths should look like
            # %i = integer, %02d = integer with two characters, for year month and day
        path = "/data/atsr/MMDB/AQUA_MODIS_L3E/2.00/%i/%02d/%02d/ESACCI-LST-L3C-LST-MODISA-LONDON_0.01deg_1DAILY_DAY-%i%02d%02d000000-fv2.00.nc" % (day.year, day.month, day.day, day.year, day.month, day.day)
        print(path)
        # os fetches contents 
        if os.path.isfile(path):    
            # open file and define lst daily mean
            f = Dataset(path)
            lst = f['lst'][:]
            # populate numpy array, lst slice ignoring first index and accounting for NaN
            data_grid.append(f['lst'][0,:,:].filled(np.nan))
            f.close()
        else: print ("%s not found" % path)

    # stack array into three dimensions (lat, lon, time)
    data_3d = np.dstack(data_grid)  

    # calculate mean of each grid point pixel and ignore NaN values 
        # make it 2d
    data_2d = np.nanmean(data_3d, axis = 2)

# define lat lon from last file 
    f = netCDF4.Dataset ('/data/atsr/MMDB/AQUA_MODIS_L3E/2.00/2018/12/31/ESACCI-LST-L3C-LST-MODISA-LONDON_0.01deg_1DAILY_DAY-20181231000000-fv2.00.nc')
    print (f)

    # define lat and lons including all indicies 
    lat = f['lat'][:]
    lon = f['lon'][:]

    # plot data
    # set up a map and size
    fig = plt.figure(figsize=(12,12))
    ax = fig.add_subplot(projection=ccrs.PlateCarree())
    # define the coordinate system that the grid lons and grid lats are on
    mapped_array = ax.pcolormesh(lon, lat, data_2d)

    # add title 
    plt.title ('Aqua Modis London', fontsize=12) 
    # set axis labels
    plt.xlabel('Latitude', fontsize=10)
    plt.ylabel('Longitude',fontsize=10) 
            
    # add features
    ax.coastlines()
    ax.add_feature(cfeature.BORDERS, edgecolor='black')
    # add colorbar
    divider = make_axes_locatable(ax)
    ax_cb = divider.new_horizontal(size="5%", pad=0.1, axes_class=plt.Axes)
    fig.add_axes(ax_cb)
    plt.colorbar(mapped_array, cax=ax_cb)

    #add lat lon grids and ticks
    gl = ax.gridlines(draw_labels=True, alpha=0.1)
    gl.top_labels = False
    gl.right_labels = False    
    # show and save   
    plt.show() 
    plt.close()

Solution

  • I would recommend using the xarray package for collecting all of your data into a single array. First put all of the netCDFs in the same directory, then merge with:

    import xarray as xr
    f = xr.open_mfdataset('/data/atsr/MMDB/AQUA_MODIS_L3E/2.00/ESACCI-LST-L3C-LST-MODISA-LONDON_0.01deg_1DAILY_DAY*.nc')
    

    From there it looks like you want the monthly mean:

    f_monthly = f.resample(time = 'MS', skipna = True).mean()
    

    Then it is pretty straightforward to loop through each of the months/years and plot.

    f_monthly.sel(time = '2018-12').plot()