Level 2 Visualization Using Python Pyresample Library (2)

Level 2 Visualization Using Python Pyresample Library (2)

Postby yiboj » Fri Nov 18, 2016 9:10 am

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.

modis_a_l2_small.gif
MODIS Aqua L2 Animation
modis_a_l2_small.gif (1.98 MiB) Viewed 4709 times


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
yiboj
 
Posts: 130
Joined: Mon Mar 30, 2015 11:22 am

Return to Visualization

cron