Simple map from NetCDF file or OPeNDAP dataset using NetCDF4

Extract a specific time step from a netcdf file or opendap dataset using netCDF4 and plot without spatial coordinates.

Import requirements, taking care to put the "%matplotlib inline" before "import matplotlib" to avoid problems on jupyterhub servers



In [1]:

    
%matplotlib inline   
import matplotlib.pyplot as plt
import numpy as np
import netCDF4
import datetime as dt

Create a dictionary to hold short titles in keys, OPeNDAP Data URLS in values



In [2]:

    
models = {'GFS CONUS 80km':'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/CONUS_80km/Best',
          'NAM CONUS 20km':'http://thredds.ucar.edu/thredds/dodsC/grib/NCEP/NAM/CONUS_20km/noaaport/Best'}

Pick a model



In [3]:

    
name = 'NAM CONUS 20km'
url = models[name]



In [4]:

    
nc = netCDF4.Dataset(url)
ncv = nc.variables
ncv.keys()









    Out[4]:





[u'x',
 u'y',
 u'time',
 u'reftime',
 u'time1',
 u'reftime1',
 u'layer_between_two_heights_above_ground_layer',
 u'layer_between_two_pressure_difference_from_ground_layer',
 u'height_above_ground',
 u'height_above_ground1',
 u'layer_between_two_pressure_difference_from_ground_layer1',
 u'layer_between_two_pressure_difference_from_ground_layer2',
 u'isobaric',
 u'LambertConformal_Projection',
 u'time1_bounds',
 u'layer_between_two_heights_above_ground_layer_bounds',
 u'layer_between_two_pressure_difference_from_ground_layer_bounds',
 u'layer_between_two_pressure_difference_from_ground_layer1_bounds',
 u'layer_between_two_pressure_difference_from_ground_layer2_bounds',
 u'Precipitable_water_entire_atmosphere',
 u'Total_cloud_cover_entire_atmosphere',
 u'Pressure_surface',
 u'Convective_inhibition_surface',
 u'Convective_Available_Potential_Energy_surface',
 u'Temperature_height_above_ground',
 u'Maximum_temperature_height_above_ground',
 u'Convective_inhibition_layer_between_two_pressure_difference_from_ground_layer',
 u'Convective_Available_Potential_Energy_layer_between_two_pressure_difference_from_ground_layer',
 u'Minimum_temperature_height_above_ground',
 u'Dew_point_temperature_height_above_ground',
 u'Wind_direction_from_which_blowing_height_above_ground',
 u'Wind_speed_height_above_ground',
 u'u-component_of_wind_height_above_ground',
 u'v-component_of_wind_height_above_ground',
 u'Relative_humidity_height_above_ground',
 u'Mean_Sea_Level_Pressure_NAM_Model_Reduction_msl',
 u'Absolute_vorticity_isobaric',
 u'Pressure_reduced_to_MSL_msl',
 u'Thunderstorm_probability_surface_3_Hour_Accumulation',
 u'Convective_precipitation_surface_3_Hour_Accumulation',
 u'Total_precipitation_surface_3_Hour_Accumulation',
 u'Probability_of_freezing_precipitation_surface_3_Hour_Accumulation',
 u'Probability_of_precipitation_surface_3_Hour_Accumulation',
 u'Percent_of_frozen_precipitation_surface_3_Hour_Accumulation',
 u'Snow_depth_surface_3_Hour_Accumulation',
 u'Temperature_layer_between_two_pressure_difference_from_ground_layer',
 u'Storm_relative_helicity_layer_between_two_heights_above_ground_layer',
 u'Parcel_lifted_index_to_500_hPa_layer_between_two_pressure_difference_from_ground_layer',
 u'u-component_of_wind_layer_between_two_pressure_difference_from_ground_layer',
 u'v-component_of_wind_layer_between_two_pressure_difference_from_ground_layer',
 u'Relative_humidity_layer_between_two_pressure_difference_from_ground_layer',
 u'Best_4_layer_lifted_index_layer_between_two_pressure_difference_from_ground_layer']



In [5]:

    
# take a look at the "metadata" for the variable "u"
var_name = 'Temperature_height_above_ground'
var = ncv[var_name]
print var









    



<type 'netCDF4._netCDF4.Variable'>
float32 Temperature_height_above_ground(time, height_above_ground, y, x)
    long_name: Temperature @ Specified height level above ground
    units: K
    description: Temperature
    missing_value: nan
    grid_mapping: LambertConformal_Projection
    coordinates: reftime time height_above_ground y x 
    Grib_Variable_Id: VAR_7-0-2-11_L105
    Grib1_Center: 7
    Grib1_Subcenter: 0
    Grib1_TableVersion: 2
    Grib1_Parameter: 11
    Grib1_Parameter_Name: TMP
    Grib1_Level_Type: 105
    Grib1_Level_Desc: Specified height level above ground
unlimited dimensions: 
current shape = (137, 1, 257, 369)
filling off



In [6]:

    
var.shape









    Out[6]:





(137, 1, 257, 369)



In [7]:

    
# Desired time for snapshot
# ....right now (or some number of hours from now) ...
start = dt.datetime.utcnow() + dt.timedelta(hours=+2)
# ... or specific time (UTC)
#start = dt.datetime(2014,9,9,0,0,0) + dt.timedelta(hours=-1)

Identify the time coordinate from the "coordinates" attribute



In [8]:

    
var.coordinates









    Out[8]:





u'reftime time height_above_ground y x '

We want time, not reftimes. Extract the index closest to the specified time



In [9]:

    
time_var = ncv['time']
itime = netCDF4.date2index(start,time_var,select='nearest')

Access data at specified time step and level



In [10]:

    
lev = 0
t = var[itime,lev,:,:]
t.shape









    Out[10]:





(257, 369)

Make a nice time stamp string



In [11]:

    
tm = netCDF4.num2date(time_var[itime],time_var.units)
daystr = tm.strftime('%Y-%b-%d %H:%M UTC')

Make a pcolormesh plot of the data



In [12]:

    
fig=plt.figure(figsize=(18,10))
plt.pcolormesh(t)
cbar = plt.colorbar()
plt.title('%s:%s, %s, Level %d' % (name, daystr, var_name, lev));