Level 2 Contour Plot Issue Using Python

Programming scripts on how to read, analyze and access PO.DAAC data

Level 2 Contour Plot Issue Using Python

Postby yiboj » Wed May 16, 2018 8:08 pm

There is an issue/bug with python matplotlib package with the filled contour plot as shown in Figure 1. The horizontal lines runs across the map when the longitude data values are mixed across the boundary (0 or 360 degree).
Here is the code to create Figure 1, which shows the wrong contour plot by using RSS SMAP Level 2 data on April 1, 2015.
Code: Select all
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from pylab import *

# Read in SMAP sea surface salinity L2 granule file
filename = 'RSS_smap_SSS_L2C_r00872_v02.0_70km.nc'
ncin = Dataset(filename, 'r')
lons = ncin.variables['cellon'][:]
lats = ncin.variables['cellat'][:]
sss = ncin.variables['sss_smap'][:]
ncin.close()

lats = lats[:,:,1]
lons = lons[:,:,1]
sss = sss[:,:,1]
lats = lats.squeeze()
lons = lons.squeeze()
sss = sss.squeeze()

# Setup plot
plt.figure(figsize=(2, 1.12), dpi=300)

# Contour levels
clevs = np.linspace(31, 39, 101)

cmap=plt.get_cmap("jet")
cmap.set_under("black")
cmap.set_over("DarkViolet")

m = Basemap(projection='cyl', llcrnrlon=0.0, llcrnrlat=-90,
        urcrnrlon=360.0, urcrnrlat=90)
m.bluemarble()

# Filled contour
cs=m.contourf(lons,lats, sss, clevs, cmap=cmap, extend='both')

# Create extended colorbar
cb = m.colorbar(cs, 'right', size='2%', pad='0.5%')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=4)
cb.set_label('SSS (PSU)', fontsize=4,fontweight="bold")
cb.set_ticks(range(30,40,1))

plt.title("RSS_smap_SSS_L2C_r00872_v02.0_70km.nc", fontsize=4,fontweight="bold")

plt.subplots_adjust(left=0.05, right=0.94, top=0.98, bottom=0.01)

plt.show()

l2_bad.png
Figure 1, wrong contour plot
l2_bad.png (241.46 KiB) Viewed 2007 times


The easy and fast way to fix the contour plot is to separate the data before and after the boundary (360.0 degree). User may also fix the contour plot routine contourf in the matplotlib package. Here is the code to create Figure 2, which shows the correct contour plot.
Code: Select all
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from pylab import *

# Read in SMAP sea surface salinity L2 granule file
filename = 'RSS_smap_SSS_L2C_r00872_v02.0_70km.nc'
ncin = Dataset(filename, 'r')
lons = ncin.variables['cellon'][:]
lats = ncin.variables['cellat'][:]
sss = ncin.variables['sss_smap'][:]
ncin.close()

lats = lats[:,:,1]
lons = lons[:,:,1]
sss = sss[:,:,1]
lats = lats.squeeze()
lons = lons.squeeze()
sss = sss.squeeze()

sss1 = sss.copy()
sss2 = sss.copy()

idx1 = np.where(lons>=180.0)
sss1[idx1] = np.nan

idx2 = np.where(lons<180.0)
sss2[idx2] = np.nan

# Setup plot
plt.figure(figsize=(2, 1.12), dpi=300)

# Contour levels
clevs = np.linspace(31, 39, 101)

cmap=plt.get_cmap("jet")
cmap.set_under("black")
cmap.set_over("DarkViolet")
m = Basemap(projection='cyl', llcrnrlon=0.0, llcrnrlat=-90,
        urcrnrlon=360.0, urcrnrlat=90)

m.bluemarble()


# Filled contour
cs=m.contourf(lons,lats, sss1, clevs, cmap=cmap, extend='both')
cs=m.contourf(lons,lats, sss2, clevs, cmap=cmap, extend='both')

# Create extended colorbar
cb = m.colorbar(cs, 'right', size='2%', pad='0.5%')
cb.ax.set_yticklabels(cb.ax.get_yticklabels(), fontsize=4)
cb.set_label('SSS (PSU)', fontsize=4,fontweight="bold")
cb.set_ticks(range(30,40,1))

plt.title("RSS_smap_SSS_L2C_r00872_v02.0_70km.nc", fontsize=4,fontweight="bold")

plt.subplots_adjust(left=0.05, right=0.94, top=0.98, bottom=0.01)

l2_good.png
Figure 2, correct contour plot
l2_good.png (286.57 KiB) Viewed 2007 times
yiboj
 
Posts: 77
Joined: Mon Mar 30, 2015 11:22 am

Re: Level 2 Contour Plot Issue Using Python

Postby richli » Fri May 25, 2018 6:02 pm

Because basemap is deprecated, I lightly modified the script to use cartopy as well. Additionally, a context manager is used to automatically close the file and I simplified a few other minor details. The plot isn't exactly the same but the script at least illustrates one way to use cartopy.

Code: Select all
#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from netCDF4 import Dataset

# Read in SMAP sea surface salinity L2 granule file
filename = 'RSS_smap_SSS_L2C_r00872_v02.0_70km.nc'
with Dataset(filename, 'r') as ncin:
    lons = ncin.variables['cellon'][:]
    lats = ncin.variables['cellat'][:]
    sss = ncin.variables['sss_smap'][:]

lats = lats[:, :, 1]
lons = lons[:, :, 1]
sss = sss[:, :, 1]

sss1 = sss.copy()
sss2 = sss.copy()

sss1[lons >= 180.0] = np.nan
sss2[lons < 180.0] = np.nan

# Setup plot
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()},
                       figsize=(6, 4))
ax.set_global()
ax.stock_img()

# Filled contours
clevs = np.linspace(31, 39, 101)
for sal in (sss1, sss2):
    cs = ax.contourf(lons, lats, sal, clevs, cmap='magma', extend='both',
                     transform=ccrs.PlateCarree())

# Create extended colorbar
cb = fig.colorbar(cs, ax=ax, pad=0.02, shrink=0.7)
cb.set_label('SSS (PSU)')
cb.set_ticks(range(30, 40))

ax.set_title("RSS_smap_SSS_L2C_r00872_v02.0_70km.nc")

fig.savefig("smap_r00872.png", bbox_inches='tight')


(Sorry, I can't figure out how to attach the image to this post.)
richli
 
Posts: 1
Joined: Fri May 25, 2018 4:36 pm


Return to Data Recipes