Maps

We tell stories with maps about global warming, biodiversity; we can design more livable cities, track the spread of epidemics. That makes a difference. 
Jack Dangermond

A map is not the territory it represents, but, if correct, it has a similar structure to the territory, which accounts for its usefulness.
Alfred Korzybski

I wanted to plot sea surface temperature (SST). Actually, I wanted more than that. I wanted an easy way to plot global SST. I wanted a lazy way to make a nice plot. It turns out it's not hard, but more difficult than I naively expected.

Get Some Data

NOAA's National Centers for Environmental Information (NCEI) hosts public access archives for environmental data on Earth. Among the many data sets are high resolution satellite measurements of SST. The data, daily OISST, I'm interested in has a spatial grid resolution of 0.25° and temporal resolution of 1 day. It uses Advanced Very High Resolution Radiometer (AVHRR) infrared satellite SST data. The first problem that you might run into using NOAA data is that the NOAA website has dead links. The link from the first page to the OISST data is dead. The second problem is that the data comes in two different formats. Prior to April 1, 2017 OISST data is in NetCDF3, After that date, the data is in NetCDF4.

NetCDF (Network Common Data Form) is a machine-independent data format for array-oriented data. It's a complicated format, but there are libraries for reading NetCDF files and extracting the data arrays. The file format contains extensive meta-data describing the data. 

 For R, there is the netcf4 library. This code will read a NetCDF4 file whose name is stored in nc_files. It extracts SST, latitude, and longitude data from the file.

  require("ncdf4")
  nc=nc_open(nc_filename)
  temp=ncvar_get(nc, 'analysed_sst')
  lon=ncvar_get(nc, 'lon')
  lat=ncvar_get(nc, 'lat')
  nc_close(nc)   

Python has the netCDF4 library.

    import netCDF4
    file = netCDF4.Dataset(filename)
    lat = file.variables['lat'][:]
    lon = file.variables['lon'][:]
    data = file.variables['analysed_sst'][0,:,:] - 273.15
    file.close()

According to the netCDF4 docs, I should be able to read a url with the netCDF4.Dataset constructor rather than a file name, but when I used the url of a file that I know to exist, I get a file not found error. It's easy enough to work around that problem using urllib.

Plotting SST in R

I downloaded some SST data for Jan. 24 2018. My first attempt was to use the R maps package. Once you have the SST data, it pretty easy to create a contour plot using maps. Here's the whole thing.

plot_world = function(nc_filename, title='') {
  # filename - the path to an NetCDF4 file
  # title - an optional title passed to plot
  
  require("ncdf4")
  require('maps')
  
  # get the SST data
  nc=nc_open(nc_filename)
  temp=ncvar_get(nc, 'analysed_sst')  # SST 
  lon=ncvar_get(nc, 'lon')
  lat=ncvar_get(nc, 'lat')
  nc_close(nc)
  
  # plot the contours
  maps::map(database="world", fill=TRUE, col="light blue")
  maps::map.axes()
  filled.contour(x=lon, y=lat, z=temp-273.15,
                 col=rev(heat.colors(24)),
                 key.title=title('Temp C'),
                 plot.title=title(main=title, xlab='Longitude', ylab='Latttude'))
}

It's easy, but the results aren't completely satisfying.

Maybe I could do better with ggplot2. Here's an attempt.

ggplot_world = function(nc_filename, title='') {
  require("ncdf4")

  # get the SST data
  nc=nc_open(nc_filename)
  temp=ncvar_get(nc, 'analysed_sst')  # SST 
  lon=ncvar_get(nc, 'lon')
  lat=ncvar_get(nc, 'lat')
  nc_close(nc)
  
  # create a df to hold the data. Each temp has a latitude and logitude
  sst = data.frame(lon=rep(lon,each=dim(temp)[2]), 
                   lat=rep(lat, dim(temp)[1]), temp=as.vector(t(temp-273.15)))
  
  require(ggplot2)
  mapWorld = borders("world", colour="gray50", fill="gray50")
  mp = ggplot() + mapWorld
  pm<- mp + geom_tile(data=sst, aes(x=lon, y=lat, fill=temp)) + 
    scale_fill_distiller(palette="Spectral") +
    labs(title = title, x = 'Latitude', y='Longitude', fill='Temp.')
  pm
}

Here's what it produces.

Plotting SST in Python

I could probably tweak the ggplot2 code to make the colors look better but I decided to try Python.

Python has a map library called Basemap for plotting maps with matplotlib. It's pretty easy to use and has a large number of features like different map projections.

Here's some code, partly adopted from here.

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import netCDF4
import urllib
import tempfile
import sys
import argparse

def main():
    url = GetArgs()   # get the data url from the command line

    # use a temp file to hold teh downloaded data
    tf = tempfile.NamedTemporaryFile()
    filename, headers = urllib.request.urlretrieve(url, tf.name)

    # read the NetCDF4 file
    file = netCDF4.Dataset(filename)
    lat = file.variables['lat'][:]
    lon = file.variables['lon'][:]
    data = file.variables['analysed_sst'][0,:,:] - 273.15
    file.close()

    # init Basemap projection (Miller Cylindrical)
    # set max and min lattitude and longitude and resolution.
    m = Basemap(projection='mill', llcrnrlon=lon.min(), urcrnrlon=lon.max(), llcrnrlat=lat.min(), urcrnrlat=lat.max(), resolution='l')

    # draw the basic map
    m.drawcoastlines(linewidth=0.1)
    m.fillcontinents(color='None', lake_color='White')
    m.drawmapboundary(fill_color='None')

    m.drawparallels(np.arange(-90, 90, 30), labels=[1, 0, 0, 0], linewidth=0.1, color='k', fontsize=10)
    m.drawmeridians(np.arange(-180, 180, 60), labels=[0, 0, 0, 1], linewidth=0.1, color='k', fontsize=10)

    # plot the contours
    Lon, Lat = np.meshgrid(lon, lat)
    x, y = m(Lon, Lat)
    cs = m.pcolormesh(x, y, data, shading='flat', cmap=plt.cm.jet)
    cbar = m.colorbar(cs)
    cbar.outline.set_linewidth(1)
    cbar.ax.get_yaxis().labelpad = 15
    cbar.ax.set_ylabel('Temperature Deg. C', rotation=270)
    plt.title('SST 2018-01-24', fontsize=9)
    plt.show()

Here is the result using Jan. 24 data.

The Python code can be downloaded here.