Python | Tutorial: Intro to Cartopy

Introduction:

There are currently two main Python libraries for plotting geographic data on map: Cartopy and Basemap.  New users should use the Cartopy since support Cartopy will replace Basemap and support for Basemap is expected to wrap up in 2020.  Therefore, this tutorial will focus on Cartopy.  If you want to use Basemap (e.g., for plotting in "3D"), here are a couple of links to existing tutorials and examples to help get you started.  Furthermore, many concepts described here may be useful for understanding Basemap.

• Basemap Examples
• Basemap Tutorials 

Cartopy:

Cartopy is a geospatial plotting library built on top of Numpy and Matplotlib that makes plotting gridded data, shapefiles, and other geographic data on over 30 different map projections.  Furthermore, the Cartopy "transform" functionality makes it straightforward to convert data from one projection to another.  In this tutorial you will learn how to use Cartopy to:

• Plot GFS surface temperature data on a map
• Mask out land surfaces to plot a map of sea surface temperature and ice-cover from the GFS
• Plot a regional map of surface temperature with US states
• Use the transform function to plot GFS data on a different map projection

This tutorial accesses the NOMADS data server using the netCDF4 library, if you are unfamiliar with doing this, I recommend you see my tutorial on reading NetCDF data.

Getting Started - Grabbing GFS Data:
The following code is used to get started with this tutorial, essentially, it imports all of the required modules, pulls the GFS data from the NOMADS data server, and loads it into arrays.  Once the data is loaded, the remainder of the tutorial will focus on Cartopy. 

import numpy as np
from matplotlib import pyplot as plt
from cartopy import crs as ccrs

import cartopy.feature as cfeature

import netCDF4 as nc
import datetime as DT 

datetime=(DT.datetime.utcnow()-DT.timedelta(hours=6)).strftime('%Y%m%d')
fhour=20

nomads_path="https://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs{}/gfs_0p25_00z".format(datetime)

datafile=nc.Dataset(nomads_path)
lat=datafile.variables['lat'][:]
lon=datafile.variables['lon'][:] 

sfcT=datafile.variables['tmpsfc'] #Keep meta data for now!
sfcT_data=sfcT[fhour,:]

title=sfcT.long_name
time=datafile.variables['time']
tstart=DT.datetime.strptime(time.minimum,'%Hz%d%b%Y')
timestr=(tstart+DT.timedelta(hours=time.resolution*24.*fhour)).strftime('%c')

Now that your data is loaded, we initialize a figure, and subplot using matplotlibs "projection" keyword to link to the Cartopy library.  Once you've attached the Cartopy package to your subplot, there are a number of additional objects and classes you can add to your figure.  In the first example, I'll demonstrate by calling upon the "coastlines" attribute to plot continents on the map.

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())

Now, your figure and subplot are defined.  The Cartopy "PlateCarree()" projection is a basic cylindrical map projection that accepts coordinate information in the form of lat/lon coordinate pairs.  The lat/lon data can be either 1D or 2D, but essentially, if your coordinate information is lat/lon, your projection is PlateCarree().   You can pass 2 keywords into PlateCarree (central_longitude and globe), but in most instances, you won't need to.   Occasionally, I reset the central_longitude to 180, instead of it's default zero.

The code to plot the map is:

Z=plt.pcolormesh(lon,lat,sfcT_data,cmap='jet')

pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines()

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

And you'll get an image similar to this one:

GFS surface temperature from Cartopy

Now, lets say you wanted to mask out land areas, and add sea-ice to the map, to focus on the oceans.
Adding the sea-ice is easy, all it requires is for you to pull the icecsfc variable from the GFS data, and plot it overtop of the surface temperature map:

plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)

We use the numpy mask (ma) functionality to maskout gridcells where the ice cover is less than 10%.
To mask out land areas, we make use of the "Cartopy Feature" module, which handles shapefiles.  The Cartopy Feature module connects seemlessly to the Natural Earth Dataset, and allows for anyone to plot a host of different GIS datasets without needing to download individual shapefiles before hand.  The interface even has many pre-defined functions, including land.  This makes masking out the land areas, a breeze.

ax.add_feature(cfeature.LAND,zorder=4,color='gray')

In the above code, the "add_feature" function is used, and the pre-defined "LAND" feature from the Cartopy Feature package is supplied with a couple of basic keyword arguments.  Putting it all together:

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.)
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)
ax.add_feature(cfeature.LAND,zorder=4,color='gray')
pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines(color='k',zorder=5)

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

Running the above code will give you an image like this:

GFS SST and Sea-ice cover

Now, lets zoom in over the United States, and add states to the map.  To clip the extent, the "set_extent" function is used:

ax.set_extent([-125, -65, 25, 55],crs=ccrs.PlateCarree())

The list corresponds to lon/lat boundaries from [western-most longitude, eastern-most longitude, southern-most latitude, northern-most latitude] and the crs argument indicates which transform is being applied.  To add the states, you need to define the states from Natural Earth database. the "scale" argument

states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                         name='admin_1_states_provinces_shp',edgecolor='k')

 By changing the extent, and adding the states, you'll get a map that looks like this:

GFS SST centered over the United States

You'll notice that the "LAND" feature masks out the Great Lakes, which is unfortunate, and to my knowledge the only way to get both the land mask and the SST over lakes using GFS surface temperature data, is to use a Lakes shapefile with shapefile clipping, a technique that is beyond this simple introduction, but may be the subject of a future tutorial.

Finally, if you want to use a different map projection to plot the data, e.g., a North Polar Sterographic projection, you simply define a different projection when setting up your subplot, and use the "transform" keyword in the plotting function to specify that you are transforming your lon/lat coordinates from the PlateCarree() projection:

ax=plt.subplot(111,projection=ccrs.NorthPolarStereo())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.,transform=ccrs.PlateCarree())
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3,transform=ccrs.PlateCarree())

ax.set_extent([-180, 180, 45, 90],crs=ccrs.PlateCarree()) 

Note, that we also reset the extent to focus on the Northern Latitudes.  You should get something that looks like this:

GFS Sea-ice and SST Polar Plot

That is the extent of this introductory tutorial on how to use Cartopy to work with atmospheric data.  You can learn more from the Cartopy Website, and get a list of different map projections here Cartopy Map Projections. Future tutorials on Cartopy will focus on working with shapefiles, transforming individual points, how to work with data without lat/lon coordinates, and how to put fine details on your map.

Read More

Python | Tutorial: Introduction to Numpy

While this post has (almost) nothing to do with atmospheric science, it is an important introduction to the Numpy package, which is how you will be storing and working with data in Python.  This tutorial is heavily geared towards beginners and those who want to strengthen their background understanding of Numpy.  If you're already comfortable using Numpy to work with data, there probably isn't a whole lot new in this tutorial for you to learn, but who knows.

As stated in the introduction to Python post, Python stores vector data in lists, tuples, and dictionaries.  Numpy is an all encompassing package that transforms these types into numerical arrays and allows you to carry out efficient, vectorized, and higher order numerical functions.  To demonstrate using a simple example, if you wanted to add to vectors stored as Python lists:

list1=[5,5,5]
list2=[5,5,5]

list3=list1+list2

You end up with the result:

list3=[5, 5, 5, 5, 5, 5]

Which is not addition, but rather list concatenation.  To actually add the numbers together to get a vector of length three with "10s" for each value, you need to loop through the list.  To demonstrate how to do this with lists here is simple function that uses the zip function to access both lists at once:

def list_add(list1,list2):
    list3=[]
    for i,j in zip(list1,list2):
        list3.append(i+j)
    return list3

In the above function, a third list (creatively named "list3") is initialized with no values, and a "for loop" is used to loop through the input lists and add the list values together, storing them in list3.

list3=list_add(list1,list2)

result: [10, 10, 10]

which is the result we were looking for.  Great, we can now add vectors together using lists.  Although, addition is a pretty simple thing, and you can probably imagine the "list_add" function starts to get really complex as you grow to multidimensional arrays.  Enter Numpy.  Numpy allows you to define arrays, which behave much more intuitively for data analysis.  I fact, I don't think it's too much of a stretch to say that, almost 90-95% of the data structures you will ever work with in Python will involve using Numpy arrays.  Building off of the above example, if we repeat the exercise, only replacing "list1" and "list2" with Numpy arrays, call them "array1" and "array2":

import numpy as np

array1=np.array([5,5,5])
array2=np.array([5,5,5])

print(array1+array2)

result: [10, 10, 10]

Much simpler, in contrast to adding the two lists together, which required defining a function and using a for loop, simply adding the two arrays gives the correct solution.  Additionally, since Numpy is optimized, it is way more efficient to use Numpy array operations than it is to use loops.

Numpy array intrinsic attributes and functions

What's great about Numpy, is that it takes full advantage of Python object oriented framework when storing data.  Specifically, when you define a Numpy array, it isn't only a block of numbers, it's an entire data object with attributes and functions attached to it that make it more intuitive to work with.  To briefly demonstrate, a couple of the more useful attributes, a two-dimensional array will be defined using the np.ones function (which defines an array of ones given specific dimensions):

ones=np.ones((10,5)) #Note passing the array dimensions as a list instead of a tuple works too

This will define an array of 5 columns and 10 rows, we can confirm that using the shape attribute:

print(ones.shape)

result: (10,5)

You can also find out what kind of array it is using the dype attribute (e.g., is it an integer or floating point array).

print(ones.dype)

result: float64

You can also easily transpose the array:

ones_transpose=ones.T
print(ones.shape,ones_transpose.shape)

result:(10,5),(5,10)

Or define a new array with the same size and attributes:

twos=np.ones_like(ones)*2

This is just a small sampling of the features and functions intrinsic to any defined Numpy array, that I happen to find most useful.  There are dozens more that I won't describe here.

Numpy: My Top 5 Numpy functions for Atmospheric Science Research.

Perhaps this heading should be at the very top of this blog post, since the previous sections mostly dealt with trying to provide a cursory explanation of what makes Numpy distinct and powerful.  These are the 5 functions I use the most, not including the numpy.ma (masking) package which is extremely useful, but needs a dedicated blog post to adequately describe.

1. Numpy linspaceand arange

   Two functions for the price of one.  Basically, these functions allow you to quickly create vectors that span a range.  The difference between the two functions is best summarized as follows: 
   
   • arange takes up to three arguments: start, end, increment
   • linspace also takes up to three arguments: start, end, vector length 



lines=np.linspace(0,10,11)
arange=np.arange(0,10,2)

print(lines)
print(arange)

result (linspace): [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
result (arange): [0 2 4 6 8]


A few things to unpack from this example, note that the arange vector stops at 8, despite the second argument being 10.  This is because the non-inclusive nature of Python for the right hand side:  Basically, arange is incrementing 2 and stops before it hits 10.  Linspace on the other hand just builds a range interpolating between 0 and 10 that has 11 values.


2. Numpy Meshgrid

   I use meshgrid a lot when I need a 2D spatial array for mapping purposes, but only have two 1D vectors that define my coordinate system.  I don't know of a much better way to describe it beyond that, so I'll just throw up an example.
   
   

   lon=np.linspace(0,360,100)
   lat=np.linspace(-90,90,50)
   
   lon,lat=np.meshgrid(lon,lat)
   
   print(lon.shape)
   
   result: (50,100)
   

   
   In this example, the 1D vectors are rebroadcasted into 2D arrays taking on the shape of the original 1D vectors in each dimension.  Now, you have 2D coordinate information that is often required for certain functions (e.g., 2D interpolation, or regridding) Be warned: when using this function, if you redefine lat/lon as 2D arrays (as in this example), and you for some reason try to run meshgrid again, for instance as part of a loop or in a notebook setting, it's very easy to run out of memory. 



3. Numpy Concatenate

   Concatenate is a great way to build an array without knowing it's dimensions up front.  For example, lets say you have a collection of netCDF files with gridded data, and all of these files have matching spatial coordinates, but each file has a different number of time variables, and different time coordinates.  If you want to put all of this information into an array, you can't simply define the array size up front, you need to add to the time dimension as you loop through each file.  Concatenate allows you to simply take an array and append it to another array along a specified dimension.  To simplify the example here: 
   

   array1=np.ones((5,50,40))
   array2=np.ones((3,50,40))
   array3=np.ones((7,50,40))
   
   array_total=np.concatenate((array1,array2),axis=0)
   array_total=np.concatenate((array_total,array3),axis=0)
   
   print(np.shape(array_total))
   
   result: (15,50,40)
   

   
   See that the end result has 15 in the 0 axis dimension (5+3+7).  Note that ALL axes except the specified axis must match, or the function will crash. While concatentate (and it's relative functions hstack,vstack and dstack) are convenient, beware, it is still much better to define your entire array up front and load data into it if you know the array size.  This is because Concatentate is relatively inefficient because it requires that the entire array is reallocated into memory every time you run it, which can get really slow for large arrays.

4. Numpy Mean and Nanmean

   These are pretty self explanatory.  These functions simple take an average of an array, allowing you to specify an axis as well.  The "axis" command is really handy if you want to do a time-average while preserving the spatial dimensions, for instance.
   
   

   arr1=np.linspace(0,20,5)
   print(np.mean(arr1))
   
   result: 10.0
   
   ## set the 3rd index to a nan
   arr1[3]=np.nan
   print(np.nanmean(arr1))
   
   result: 8.75  

   
   As you noticed, once you change the 3rd index value from 15 to a nan, the function now takes the mean 0,5,10,20 (which is 8.75) instead of 0,5,10,15,20 (which is 10).  However, if you try to run plain "mean" on an array with nans in it, the returned value will also be a nan.

5. Numpy diff and gradient

   
   Theses functions are pretty similar, in that they allow you to take a gradient of a function.  diff is the simpliest of the two functions listed here, and just returns Δx along a specified dimension.  This will reduce your array length along the specified dimension by 1.  Gradient is similar, but comes with a few extra bells and whistles:
   
   • Gradient preserves the array length by using numerical approximations along the boundaries
   • Allows you to specify a sample distance for each axis (i.e., the dx for dy/dx)
   • If no "axis" argument is supplied will compute the gradient in ALL dimensions and return a tuple of length n-dimensions with each index corresponding to the gradient along that index.
   Most often, I use "gradient" especially if I need to multiply the gradient against something because I can be somewhat lazy, and I don't want to have to deal with different sized arrays/vectors. The example I will provide will be show how to use these functions to approximate the derivative of "sin(x)" which is analytically (of course) "cos(x)".  Note that this example uses matplotlib.
   
   

   import numpy as np
   from matplotlib import pyplot as plt
   
   x=np.linspace(0,2*np.pi,40)
   array=np.sin(x)
   dx=x[1]-x[0]
   
   xstag=(x[:-1]+x[1:])/2.
   
   plt.figure()
   plt.plot(x,np.gradient(array,dx),label='d/dx(sin(x)) (gradient)',lw=4.4)
   plt.plot(x,array,label='Sin(x)')
   plt.plot(xstag,np.diff(array)/dx,label="d/dx(sin(x)) (diff)",ls='--',color='r')
   plt.legend()
   plt.show()
   

   
   This will generate the following figure: 

   

   sin(x) and it's approximated derivative (cos(x))

   
   

   
   

   
   

   
   

   
   

   
   

   I played with the line colors and thicknesses to make it easier to differentiate the "diff" approximation from the "gradient" approximation.  Note that to correctly plot the "diff" function I had to create a new version of "x" called "xstag" which basically shifted the x-center points 0.5dx, and shortened the array length (n) by one simply by averaging the first n-1 array elements with the last n-1 array elements.  Also note, that I obtained dx by simply differencing the first two array elements.

Final Notes:

While this post had little to do with atmospheric science, I hope it provided a useful introduction to the Numpy object class.  For more information on Numpy and its capabilities please visit the official site.

Read More

Python | Tutorial: Introduction 1 - Introduction to Python

Introduction:

While I want to keep the focus the Python section of this blog specifically on using Python for meteorological and atmospheric science applications, and not a broad tutorial on Python in general, a few basic Python syntax and coding strategies well greatly help in understanding the tutorials presented here.  This post is very much geared towards people who are new to Python, so if you are already comfortable using Python, and/or you aren't into reading, you can probably skip this one.

The Basics (a little background on Python):

• Python is NOT a free form language, that is it has fairly strict indentation and formatting rules that have to be followed in order for your program to work.  If you're new to python (or other programming languages with these restrictions), this can be extremely frustrating at first, but you get use to it over time.  Also, almost every programming specific text editor (e.g., atom) automatically adheres your code to the python requirements. 
• Like GrADS, Python "commands" can be run within the command line, or commands can be written in Python "scripts" with the suffix ".py" 
• Python is an interpreted language, instead of compiled like, for example, Fortran.  For the purposes of the tutorials and scripts presented in this blog, the only practical application of this information is that Python scripts will run slower than compiled Fortran programs, making it less ideal for applications such as numerical modeling.
• Most data in Python is stored in lists, tuples, and dictionaries.  Many Python modules restructure these data types into "arrays" which is extremely useful for data-analysis, but it's important to know that these arrays are really lists/tuples/dictionaries at their cores. I don't want to get too far into the weeds with respect to the differences between these data structures in this post, and furthermore there are numerous tutorials available via a quick google search.  However,  I use lists and dictionaries all the time when using Python for atmospheric science applications, so it is definitely worth taking the time to learn a little bit more about these data structures.  Here are a couple of links to help get you started: Lists, Dictionaries.  Tuples are everywhere in Python for a number of different reasons, and while I encounter them relatively often, I rarely ever do so in such a way that I have to consciously be aware of the fact I'm working with a Tuple.
• Numpy, Scipy, and Matplotlib:  Most data and array operations and plotting will be performed using these three modules.  Unlike Fortran, these modules perform operations in row major, which may have no practical bearing on your usage.  However, when dealing with large datasets, it may be advantageous to store your data with row major structuring in mind. 
• Most geographic mapping is performed with either Matplotlib Basemap, or Cartopy.  While many folks like Basemap, it's slowly being phased out for Cartopy, so if you are new to Python, I strongly encourage you start using Cartopy instead of Basemap.
• Finally, Python is an extremely versatile and powerful language that is relatively easy to learn.  It has countless data analysis applications, and numerous packages developed specifically with atmospheric science in mind.  Refer to (Link to getting started page) for more information on how to get started with Python. 

Common Syntax:
This list is by no means an exhaustive list, but is just a quick reference guide on some of the basic syntax that is used in Python.

• Arrays and lists use square brackets ([]) and tuples use standard parenthesis.  Array elements are referenced using the ":" as the indexing wildcard.  For example: array[:].
• The number sign, or hashtag (#), is used to denote a comment.
• For loops and if statements do not require an "enddo" or "endif" statement at the end of the loop, simply return to the outer indent (see example below).
• Important:  Python accomplishes this by requiring "smart indent" coding practices ... which will drive you nuts if you're unfamiliar with it, but you'll learn to love it since it forces you you write more readable code.
• print statements follow this syntax: print("hello world").  print "hello world" will also work in versions 2.x.  

Example Code:

Simple syntax example:

## This simple script requires no imported modules, and demonstrates the for loop,
## Some basic list structure and conditional statements.
list=[1,2,3,4,5]

for i in list:
    print(i)

list.append(6)  ## add a 6 to the end of the list

print('--------------')
for idx in range(len(list)):
    print(idx, list[idx])
print('--------------')
list.append('Hello') ## add a string to the end of the list.

for i in list:
    if isinstance(i, str) == True:
        print(i, "is a string")
print('--------------')

Arrays and Plotting:

import numpy as np
from matplotlib import pyplot as plt

x=np.arange(0,20,1)*2.*3.141592/10
y=np.sin(x)

plt.figure(figsize=(10,5))
ax=plt.subplot(1,2,1)  ## 1 = one row, 2 = 2 columns, 1 = set 1st plot.
plt.plot(x,y) 
plt.plot(x,y,ls='none',marker='o')

## Now a 2D Plot with Wind Barbs ##

x1, y1 = np.meshgrid(x, x) 
u, v = 1.5*x1, 1.5*y1 
Z=np.sqrt(u**2.+v**.2)

ax=plt.subplot(1,2,2) ## 1 = one row, 2 = 2 columns, 2 = set 2nd plot.
ax.barbs(x1,y1,u,v,Z,length=4.5,cmap='jet') 
plt.show()

The above code should produce this image:

Figure produced by code above

def main():

As stated above, Python is an interpreted language and operations are performed from top-to-bottom.
While, this can be advantageous under many circumstances, one disadvantage is that functions, or objects must be defined above their use in the code.  Some may find it awkward to define a functions either on the fly throughout the code, or all in once place at the top.  There are two ways around this 1) You can define all of your functions in a separate Python (.py) script and import it into the script you want to call the functions from, or 2) you can define a "main" function and direct the program to call the main function:

The below code will result in an error:

import numpy as np

batman()

def batman():
    for j in range(8):
        print(np.log(-1))
    print(" BATMAN!")

The below code will run just fine:

import numpy as np

def main():
    batman()

def batman():
    for j in range(8):
        print(np.log(-1))
    print(" BATMAN!") 

if __name__ == '__main__':
    main() 

In some blog posts, I may use the above syntax, especially if I have a lot of functions to define, and in others I will use the simple top-to-bottom format, defining variables as I go.  It's just a matter of user preference and readability.

Final Notes:

While this post definitely has a lot of background text that might not be relevant to you as a scientist aiming to write "research grade" code to do cool analysis, hopefully, this short intro to Python will help you better understand how Python works at a more base level and familiarize you with some of the common vocabulary you'll encounter with the language.

Read More

Tutorial: The difference between self describing and non-self describing NETCDF files, and how to write them

In writing posts for this blog, I try my best to avoid what could be considered "mission creep," or writing coding tutorials that at best only tangentially related to GrADS.  This tutorial is verging upon what could be considered mission creep and has as much to do with teaching "good data" practices as it does with teaching GrADS skills.  Furthermore, there are no GrADS scripting lessons to be learned here, and the only code I provide is written for Python. 

All that said, in this tutorial, I will be showing you how to write NETCDF files that can be read using GrADS using both the self describing 'sdfopen' and non-self describing 'xdfopen' commands. I do this because, you may one day find yourself with some gridded data set that you want to plot using GrADS, however this data set may be in a text file format, or perhaps maybe you need to perform higher level calculations on the data that cannot be performed in GrADS before you display it.  While you could write the data out as a binary data file with a GrADS control file, maybe you are more comfortable with the NETCDF format (I know I certainly am).  Plus there are other benefits to using NETCDF instead of binary; one example is cross compatibility. NETCDF files can be read in with pretty much any data analysis software, so not only can you read your file in GrADS, but you can read it in with Matlab, IDL, R, Python, etc., as well.

What makes a NETCDF file self describing?

It's all about the metadata!  Metadata is often referred to "data about the data."  For example, metadata includes things like:

• The data and time a file was created
• What data is in the file
• Units of variables within the file
• Information about the coordinates (e.g., lat, lon vs. x,y)
• If the data file is for model output it might include information such as 

• Resolution, parameterizations, time step

Metadata is super important, and learning how to include it in any data set you are producing/sharing is a critical skill in a scientific career.  We are lucky in meteorology because there have been solid efforts to standardize metadata practices in the past couple of decades.  Now, NETCDF makes it very easy to provide good metadata, so there is no excuse not to do it.  So if nothing else, I hope to get that point across to you with this tutorial.  To really drive this point home, here is a meme.

 

To make a file self-describing (that is such that GrADS can read the file in with NO accompanying control file required using the 'sfdopen' command) the metadata including in the NETCDF file must adhere to the COARDS data standards: More information on COARDS here.  Personally, I didn't find this website all that useful when figuring out how to make self-describing NETCDF files, but it does provide some nice background.

How to write a self describing NETCDF file (SDFOPEN)

Here I will provide the key metadata required to write a self-describing NETCDF file.  The example, I provide is sub-setting the first 24 hours of GFS temperature data from the 00z December 30th 2015 run.  If you choose to try and reproduce these results, you will need to update the date.  I will us python to read in the data from the GrADS Data Server (GDS) and to write the local NETCDF file.  You could use any program of your choice to do this, so long as it reads and writes NETCDF files.  The python script is included in this tutorial.

So the most important thing to understand when writing your datafile is that your coordinate dimensions and accompanying dimension variables (e.g., lat,lon) adhere to the COARDS guidelines, otherwise GrADS will not open your file.  For example, you might get an error like:

     gadsdf: SDF file has no discernable Y coordinate

The metadata for the non-dimension variables (e.g., temperature, u,v) is more up to the person writing the file.

Assuming you know how to read in the NETCDF file and data from the GDS, in python:

   datafile='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs20151230/gfs_0p25_00z'
   cdata=netCDF4.Dataset(datafile,'r')
 
To start writing the file, we will first look at the dimensions and dimension variables.
You must include the following dimensions and variables in your NETCDF file:

• lat
• lon
• time
• lev

 In python you will make these dimensions and variables as such:

   newdata=netCDF4.Dataset(newfile,'w')
   newdata.createDimension('lat', len(lat[lmin:lmax]))
   newdata.createDimension('lon', len(lon[lomin:lomax]))
   newdata.createDimension('lev', len(lev[:]))
   newdata.createDimension('time', 8)

Where lmin,lmax,lomin,lomax are indices that mark your subset boundaries.
Then in addition to creating the dimensions, you must make the variables as well (For example latitude):

   latvar=newdata.createVariable('lat','f8',('lat'))

Now you need to populate latvar with metadata.  To determine what metadata you need, I simply printed out the 'lat' variable from within the GDS datafile as this datafile is self-describing and adheres to the COARDS standards.  If you've seen my other tutorials I often access these files using the 'sdfopen' command.  In doing that, I get the following (again this will look different if you aren't using python):

<type 'netCDF4.Variable'>
float64 lat(lat)
    grads_dim: y
    grads_mapping: linear
    grads_size: 721
    units: degrees_north
    long_name: latitude
    minimum: -90.0
    maximum: 90.0
    resolution: 0.25
unlimited dimensions:
current shape = (721,)

These are the metadata that you need to include for the latitude variable in your NETCDF file.  The dimensions will change, because you are sub setting the data, but things like 'grads_dim' and 'grads_mapping' need to be the same.  In python adding metadata to a NETCDF variable is easy:

   latvar.grads_dim='y'
   latvar.grads_mapping='linear'
   latvar.grads_size=str(len(lat[lmin:lmax]))
   latvar.units='degrees_north'
   latvar.long_name='latitude'
   latvar.minimum=str(lat[lmin])
   latvar.maximum=str(lat[lmax])
   latvar.resolution='0.25'
   latvar[:]=lat[lmin:lmax]

Again, remember with the sub setting, the lmin and lmax indices are required to match the output file.
Now if I print out MY latitude variable within the output NETCDF file it looks similar to the GDS file with only minor differences:

<type 'netCDF4.Variable'>
float64 lat(lat)
    grads_dim: y
    grads_mapping: linear
    grads_size: 200
    units: degrees_north
    long_name: latitude
    minimum: 10.0
    maximum: 60.0
    resolution: 0.25
unlimited dimensions:
current shape = (200,)

You will need to do the same thing for the other dimension variables (lon,lev,time) with a few minor differences.  I don't explicitly show the steps here to save space and redundancy.  You can find the steps written down in the accompanying .py file.

Once you have your dimension variables, you need to add the temperature variables.  In the example, I have you add both the 2 meter temperature and the temperature as a function of pressure level.  This is to show you how to add variables with differing dimensions.

As I stated earlier, there are fewer restrictions for these variables than the dimension ones.  But it is still a good practice to add metadata.  When printing out the GDS 'tmp2m' var, you get this:

<type 'netCDF4.Variable'>
float32 tmp2m(time, lat, lon)
    _FillValue: 9.999e+20
    missing_value: 9.999e+20
    long_name: ** 2 m above ground temperature [k]
unlimited dimensions:
current shape = (81, 721, 1440)

So for our variable, we will try to match this as close as possible:
  gtmpvar=newdata.createVariable('GTMP','f4',('time','lat','lon'),fill_value=9.999E20)
  gtmpvar.long_name='2 Meter Temperature'
  gtmpvar.units='K'
  gtmpvar[:]=var1[:8,lmin:lmax,lomin:lomax] ## Add the actual data from the GDS file

And that's it for the 2 meter temperature.  Now, you can do the exact same thing for the vertically varying temperature with only minor changes:
   tmpvar=newdata.createVariable('TEMP','f4',('time','lev','lat','lon'),fill_value=9.999E20)
   tmpvar.long_name='Pressure Interpolated Temperature'
   tmpvar.units='K'
   tmpvar[:]=var2[:8,:,lmin:lmax,lomin:lomax]

The big difference here, is that in the dimensions section of the createVariable function includes a dimension for 'lev'.

When you are done, close out your file and you can then open it in GrADS using the 'sdfopen' command.  Then if you display the temperature data, you should get something that looks like this:

How to read and write a NON self describing NETCDF file (xdfopen)

If you want to read a file that is not self-describing into GrADS, you need to include a control file.
Essentially, you need to treat the NETCDF file as if it were a binary data file.

Now for this tutorial, we will run through the exact same process, highlighting the differences in how we write metadata to the NETCDF file.

Starting with the dimension variables, we will no longer include the metadata that we did for the self describing file:

    newdata=netCDF4.Dataset(xdffile,'w')
    newdata.createDimension('lat', len(lat[lmin:lmax]))
    latvar=newdata.createVariable('lat','f8',('lat'))
    latvar[:]=lat[lmin:lmax]

Similarly, I'll skip the rest of the dimension variables, but as you see, all I did was define the variable and add the data to it.  For simplicity, and since it doesn't matter, I'll keep the temperature variables as they were in the self describing example.

Now the NETCDF file is created.  If however you try to open this file using the 'sdfopen' command, you will get an error.  So we need to create a control file.

The control file can be easily written using python (as in this example) or whatever program you are using.  I'm going to copy the code to write the control file below, and explain it after.

    ### Now write the control file ###
    xdfctl='C:/OpenGrADS/xdffile.ctl'
    xdflocal='xdffile.nc'
    xdf = open(xdfctl, 'w')
    xdf.write('DSET ^' + xdflocal +'\n')
    xdf.write('TITLE XDF example File 00z30dec2015 \n')
    xdf.write('UNDEF 99999.0\n')
    xdf.write('XDEF lon '+str(len(lon[lomin:lomax]))+' LINEAR ' + str(lon[lomin]) +' 0.25\n')
    xdf.write('YDEF lat '+str(len(lat[lmin:lmax]))+' LINEAR ' + str(lat[lmin]) +' 0.25\n')
    xdf.write('ZDEF lev '+str(len(lev[:]))+' LEVELS ')
    for i in lev[:]:
        xdf.write(str(int(i))+' ')
    xdf.write('\n')
    xdf.write('TDEF time 8 LINEAR 00z30dec2015 3hr\n')
    xdf.write('VARS 2\n')
    xdf.write('GTMP=>GTMP 0 99 2-meter Temperature [k]\n')
    xdf.write('TEMP=>GTMP '+str(len(lev[:]))+' 99 Temperature [K]\n')
    xdf.write('ENDVARS')
    xdf.close()

The most important take away to learn about writing a control for a NETCDF file is that you match the variable/dimension names in the control file to the variable/dimension names in the NETCDF file (case sensitive).

For example: 'GTMP=>GTMP' will only work if there is a variable in the NETCDF file called 'GTMP'.  Similarly, following the XDEF definition, you need to point to the longitude variable (lon) in the NETCDF file.  This is different than writing out control files for binary datasets which don't require anything after XDEF.

Other than that, you see some other things, the levels are defined in a special way since the vertical coordinate is not incremented linearly.  Also the lat/lon sizes are not hard coded, to allow you to make bigger or smaller spatial subsets without a lot of extra work.

Open the control file in GrADS using the 'xdfopen' command:

'xdfopen xdffile.ctl'

A few notes on using the control file:

Note that the coordinate data is not actually read
from the NETCDF file, it is inferred from the control file.  So it is important that the resolution and the starting points are right, otherwise you will get an incorrect result, even if your longitude data is correct in the NETCDF file.  For example, comparing the correct (0.25) resolution to an incorrect resolution (0.45) gives the two different maps shown on the right.  So get the resolution right!

One final note: This kind of control file only seems to work with nice rectangular datasets.  That is lat/lon coordinates that vary independently from one another.  If your lat/lon grid does not vary independently (i.e., longitude is a function of y and visa-versa), then using this method will warp tour map projection.  In these situations it is recommended that you write your data out to a binary file with the PDEF option.  I have yet to get the PDEF option to work with a NETCDF file.






So that is it for the tutorial, hopefully you have a better understanding regarding the difference between self describing and non self describing NETCDF files, and how to write and open them in GrADS.  No GrADS scripts to download here, but a full python script (version 2.7) is included that will reproduce the files and plots generated in this tutorial.  You will need the netCDF4 module installed for this script to work however.  Numpy and sys are standard I believe.  Also you may need to update the date.

Download Python Script here

Read More

Tutorial: Use basemap.gs to make a beautiful map of "SST", Sea Ice and Snow in GrADS

This tutorial will draw upon many of the skills discussed in several other tutorials on this site, e.g., how to handle multiple files at once, or how to use basemap.gs.  There isn't too much "new" information in this tutorial on how to use these different features and functions of GrADS, rather this tutorial will show you how to specifically make a very pretty map of sea-surface temperature, sea ice, and snow.  The final outcome is shown below.

For this tutorial you will need:

• color.gs 
• xcbar.gs
• colormaps.gs
• basemap.gs (plus the accompanying .asc files)

Snow Cover and Sea Ice on the Robinson Map Projection

So before we start, I want to give a full disclaimer: The plot is not actually showing any observed SST.  In this tutorial I'm going to use the 0.25 degree GFS surface temperature.  In practice the surface temperature should be roughly equal to the SST over the Oceans.

The first thing we need to do is open both the GFS data and the sea ice data.  This data will come from NOMADS and will be opened using the 'sdfopen' command.  Note, the files below are for February 2015, so keep in mind, that you will need to change the date on files if you wish to copy paste the example code below.


   'reinit'
   gfsfile='http://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs20150218/gfs_0p25_18z'
   icefile='http://nomads.ncep.noaa.gov:9090/dods/ice/ice20150218/ice.00z'
  'sdfopen 'gfsfile
  'sdfopen 'icefile


Now that the files are open we will simply set up the map.  Since we are using the Robinson Projection, we need to set the longitude to range from -180 to 180 and the latitude -90 to 90.  

   'set gxout shaded'
   'set mpdset hires'
   'set lon -180 180'
   'set lat -90 90'
   'set mproj robinson'
   'colormaps -l 272 307 0.5 -map jet' ;*Note the use of colormaps.gs


Then we display the variable:

   'd tmpsfc'
   'xcbar -fs 4'
  
After a moment, the surface temperature will be displayed following the "jet" color map.
Now, we are going to mask out the land using the land basemap from basemap.gs.

Note, in my version of basemap.gs, I did not want to include the frame around the map projection.  So I opened up basemap.gs and commented out the lines at the bottom that draw the rectangle:

   * Draw a new square frame around the plot
   * If you have 'set frame off' or 'set frame circle' before running basemap,
   * you may want to comment out the next 2 lines
   *'set line 1 1 6'
   *'draw rec 'x1' 'y1' 'x2' 'y2

We are going to use the "medium (M)" resolution option for the coastlines, and we are going to define a custom dark gray color for the land.

   'set rgb 73 80 80 80'
   'basemap L 73 0 M'         ;*L = Land, 73 = Fill Color, 0 = outline Color, M=Medium resolution

SST only!

  This will fill in the land areas with a gray map and you will have an image like this:

Now that we are done with the "SST" and the land mask, we are going to open up the sea-ice file and plot the sea ice.  It is important that you use the maskout function when plotting the sea-ice, otherwise you will overwrite all of the SST data with zeros.

  'set gxout shaded'
  'set dfile 2'                   ;*Sets Default file to file 2 (Sea ice)
  'set z 1'                        ;*Set z to 1
  'set t 1'                        ;*Set t to 1
  'set map 0 1 6'   
  'color 0 0.6 0.1 -kind dimgray->seashell->white'           ;*Note use of color.gs
  'd maskout(icecmsl,icecmsl-0.1)'                                     ;*Need to mask out ice.


After a few moments, the sea ice will plot and you will have an image like this:

Sea Ice and SST only!

Now this is all well and good, but our knowledge is limited to the ice over the oceans, we don't really get any sense of how ice and snow covers the land masses.  So the final touch that we'll add to the map is to plot the GFS snow water equivalent (SWE) using a similar color scale as ice cover.  To do that, we need to reset the file to our first file and then reset our time our vertical level.  Lastly, we plot SWE (again using maskout to avoid plotting over everything).

  'set dfile 1'
  'color 10 100 5 -kind dimgray->seashell->white'
  'set z 1'
  'set t 1'
  'set map 0 1 6'
  'd maskout(weasdsfc,weasdsfc-10)'

And after a few moments, you have a map like the one shown at the beginning of the tutorial.

You can play around with your minimum values for ice and snow, I used 0.1 for sea-ice and 10 mm for SWE.  That seemed to produce a good looking map.  Experiment as you like!  This plot is really basic, no special formatting outside of the actual plot (titles, etc).  Hopefully you enjoyed this tutorial, I know it was oddly specific, and didn't really present anything new, hopefully you learned something anyway!

Download Example Script

Read More