PODAAC

Pyresample is a Python package for resampling (reprojection) of earth observing satellite data. Pyresample handles both resampling of gridded data (e.g. geostationary satellites) and swath data (polar orbiting satellites). Pyresample can use multiple processor cores for resampling. Pyresample supports masked arrays.

Here is the plot and code using MODIS Aqua Level 2 dataset for the day September 1, 2016. Only half day of data is used here due to the image size limitation.

Code: Select all

import matplotlib.pyplot as plt
import pyresample as pr
from pyresample import image, geometry
from netCDF4 import Dataset
from pylab import *
import glob
import os

# Global Parameters
dirname = 'YOUR DATAFILE DIRECTORY'
dlon = 360.0/(24.0*60.0)*3.0
mydpi = 100

k = 0
for filename in sorted(glob.glob(dirname+'/*.nc')):
  fname = os.path.basename(filename)

  #Read init data
  ncin = Dataset(filename, 'r')
  lons = ncin.variables['lon'][:]
  lats = ncin.variables['lat'][:]
  sst = ncin.variables['sea_surface_temperature'][:]
  ncin.close()

  if k == 0:
    alat0 = np.average(lats)
    alon0 = -150.0

  alat = 23.5
  alon = alon0 - dlon*k
  alon2 = alon-180.0

  if alon < -180.0:
    alon = alon + 360.0
  if alon2 < -180.0:
    alon2 = alon2 + 360.0

  area_def = geometry.AreaDefinition('area1', 'Test1', 'area1',
        {'a': '6370998.0', 'units': 'm',
        'lat_0': ''+str(alat),
        'lon_0': ''+str(alon), 'proj': 'ortho'},
        640, 320,
        [-6370998.0, -6370998.0, 6370998.0, 6370998.0])
        #[-5326849.0625,-5326849.0625,5326849.0625,5326849.0625])
        #[-1370912.72, -909968.64, 1029087.28, 1490031.36])

  area_def_2 = geometry.AreaDefinition('area2', 'Test2', 'area2',
        {'a': '6370998.0', 'units': 'm',
        'lat_0': ''+str(-alat),
        'lon_0': ''+str(alon2), 'proj': 'ortho'},
        640, 320,
        [-6370998.0, -6370998.0, 6370998.0, 6370998.0])

  plt.figure(frameon=False)

  swath_def2 = pr.geometry.SwathDefinition(lons, lats)

  if k == 0:
    sst_cat  = sst
    sst1  = sst
    swath_def = swath_def2
    swath_def1 = swath_def2
  else:
    swath_def = swath_def1.concatenate(swath_def2)
    swath_def1 = swath_def
    sst_cat = np.concatenate((sst1, sst), axis=1)
    sst1 = sst_cat

  if k >= 0:
    left = 0.05
    bottom = 0.15
    width = 0.9
    height = 0.8
    fig = plt.figure(figsize=(380/mydpi, 290/mydpi), dpi=mydpi)
    plt.subplots_adjust(bottom=0.16, left=.02, right=.99, top=.92, hspace=.0)
    fig.subplots_adjust(wspace=0, top=0.92)
    plt.axis('off')

    bmap = pr.plot.area_def2basemap(area_def)
    bmng = bmap.bluemarble()
    result = pr.kd_tree.resample_nearest(swath_def, sst_cat, area_def, radius_of_influence=20000, fill_value=None)
    col = bmap.imshow(result, origin='upper', vmin=250, vmax=310)

    plt.subplot(122)
    plt.axis('off')

    bmap_2 = pr.plot.area_def2basemap(area_def_2)
    bmng_2 = bmap_2.bluemarble()
    result = pr.kd_tree.resample_nearest(swath_def, sst_cat, area_def_2, radius_of_influence=20000, fill_value=None)
    col = bmap_2.imshow(result, origin='upper', vmin=250, vmax=310)

    cbaxes = fig.add_axes([0.15,0.15,0.73,0.03])
    cbar = plt.colorbar(cax=cbaxes, orientation='horizontal')
    cbar.ax.set_xlabel('Sea Surface Temperature (K)',size=10)
    cbar.ax.tick_params(labelsize=8)

    plt.savefig("modis_a_l2_small_{0}".format(str(k).rjust(4, "0")), num_meridians=10, num_parallels=10, bbox_inches='tight')

    plt.clf()

  k = k + 1