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: Tips on handeling more than 1 file at a time in GrADS

While it is possible to use one control file to work with multiple data or binary files by using a few options in the .ctl format (specifically for files of identical spatial dimensions but for different time periods), I find it much easier to move around numerous .ctl files in GrADS.  This tutorial is going to show you how to open more than 1 file at a time as well as address the main hiccups encountered when handling more than one file.

To start out we will need two files to work with.  For simplicity we will use two files with the same spatial dimensions and variables.  Before we start, since we are dealing with more than one file in this tutorial, we will use the 'reinit' command to close all of our current files as to not get ourselves tripped up.  Then, we'll start by opening these two files.

    'reinit'
    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_00z'
    'sdfopen http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_12z'

Once our two files are open, we will use the 'q files' (not to be confused with q file) command, to list our open files.
 
Output:

    File 1 : GFS half degree (0.5x0.5) fcst starting from 00Z04june2013, downloaded June 04 04:32 UTC 
    Descriptor: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_h_00z
    Binary: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_00

    File 2 : GFS half degree (0.5x0.5) fcst starting from 12Z04june2013, downloaded June 04 16:33 UTC
    Descriptor: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_h_12z
    Binary: http://nomads.ncep.noaa.gov:9090/dods/gfs_hd/gfs_hd20130604/gfs_hd_12

So basically, you just get a little information regarding your open files, the most relevant of which is your file number.  Since GrADS only works with one file at a time, you need to specify which file you are working with.  This is best done using the 'set dfile' command.  This changes your default file around. So:

    'set dfile 1'  - Points to file 1, in this case 00z GFS
    'set dfile 2'  - Points to file 2, in this case 12z GFS


Once your dfile is set, you treat GrADS as you normally would with just one file.  So, we will continue this example by plotting CAPE from the GFS output, we will start with file 1.

    'set dfile 1'
    'set gxout shaded'
    'set lat 22 50'
    'set lon 230 295'
    'set mpdset hires'
    'set t 5'
    'd capesfc'

Now, lets move on to file 2, and plot CAPE at the same time as a model comparison.  However, since the time of the model initialization has switched, when we work with file 2, we will point to the actual time, not the timestep as we did in file one.  So our command will be 'set time', not 'set t'  First, we will need to get our time though, and we will do that with the 'q dims' command.

By doing that, we get our time is equal to: 12Z13MAY2013.  Now we simply switch files:

    'set dfile 2'
    'set time 12Z04JUN2013'
    'set gxout contour'
    'd capesfc'

CAPE at 12Z June 4th from 00z GFS (shaded) and 12z GFS contoured

Now, we have our comparison map, since these files were identical spatially, we did not need to change the lat/lon dimensions around.  However, I will mention that in general it is not a bad practice to reset all of your space and time dimensions every time you switch from one file to the next, just to ensure that you are plotting everything correctly.

Now, this brings me into the special case:  If both of your files have the same spatial/time dimensions (for example if you were comparing two runs of the WRF, or different GFS ensemble members), you can short cut without having to go through the 'set dfile' commands.

Similar to the GrADS array syntax you can simply put a '.filenumber' at the end of your plotted variable:

   'd var.1'  - Plots var from file 1
   'd var.2'  - Plots var from file 2


Again, only use this shortcut if your model dimensions, are identical, e.g., grid spacing, timestep, etc.  The main problem you will run into handling multiple files in GrADS is an "Entire Grid Undefined" error.  This is simply because you didn't reset the dimensions to fit with your file.  If you get this error, simply do a 'q dims' command, see where x,y,z,t fall outside of your file dimensions.

So that is all for the basics of file handling in GrADS.  Hope you found it useful!

Download Example Script 

Read More

Tutorial: GrADS Hovmöller diagram

Hovmöller diagrams are very useful for determining the speed of various atmospheric disturbances, as it contours a variable in both time and space.  Typically, the spatial variable is on the x-axis and the time variable on the y-axis.

It's very easy to plot a simple Hovmöller diagram in GrADS.  To do this, you only need to set your varying dimensions to longitude and time.  The axes will arrange longitude on x, and time on y.  In this tutorial, the example will be a Hovmöller diagram of 500mb meridonal wind from the GFS model.  The block of code below will set up and plot this data.

    'sdfopen http://nomads.ncep.noaa.gov:80/dods/gfs_2p5/gfs_2p520130428/gfs2p5_12z'

    'set parea 1.5 10.0 0.5 8.0'

    'set lat 35'
    'set lon -135 -60'
    'set lev 500'
    'set gxout shaded'
    'set display color white'
    'set t 1 33'
    'clear'
    'set yflip on'
    'd vgrdprs'
    'cbar'

A few things to note about this block of code:

1. Both latitude and Height are fixed.  This allows for longitude and time to vary
2. The command 'set yflip' flips the y-axis around so that time starts at the top, rather than the bottom of the plot (normal Hovmöller display).
3. The first number in the 'set parea' command is 1.5 giving the left side of the plot a large margin in order to contain the labeled dates on the y-axis.  A smaller margin will cut these off.

The above block of code will create this image:

Hovmöller Diagram of meridional wind at 500mb

Notice, that you can easily track the speed and frequency of the atmosphere disturbances as they propagate across the domain by taking the slope of the wind packets.

Download this script

Read More

Tutorial: Vertical Cross section with topography in pressure coordinates

Vertical cross sections are very easy to plot in GrADS, requiring only that you fix either the latitude or longitude value, and vary the height levels.  For example, an west to east cross section is easily set up using the following commands.

     'set lat 40.0'
     'set lon -125 -68'
     'set lev 1000 100'

A lot of cross sections however like to add in the topography, giving the cross section greater structure.  The method to do this in GrADS is not terribly hard or difficult, it is just a little tricky.

The first step is to fix your page area using the 'set parea' command.  This is because the in order to overlay topography onto your cross section, we need to change around the domain, but keep domains plotting in the same spot.  Setting the page area assures this.  Once your page are is set, just plot your desired cross section variable as usual.

Now it's time to overlay the topography.  In pressure coordinates, the topography can be represented by the surface pressure (surface pressure decreases with increased topographical height).  Since the surface pressure is a 2D and not a 3D variable, you have to fix the height in GrADS to the lowest model level (usually done by 'set lev 1000', or 'set z 1').  Now that you have done this, you need to fix the y-axis of the plot to match the y-axis of your cross section.  This is done with the 'set vrange' command.  The best way to do this is to set it exactly as you would have set the levs when plotting your cross section variable.

For example, if you typed in 'set lev 1000 100' for your cross section then you would type in 'set vrange 1000 100' for your topography.

This way, your topography is set to plot on the same scale as your cross section.  Now, we are ready to plot.  To do this, we make use of the 'set gxout linefill' option.  Then we simply fill the gap between the surface pressure and the baseline pressure of 1000.  The block of code to do this looks is shown below.

    'set lev 1000'
    'set vrange 1000 100'
    'set gxout linefill'
    'set lfcols 1 1'
    'd pressfc/100;const(pressfc/100,1000)'

The 'set lfcols' command is the syntax for setting the color of the topography.  The end result for an east west cross-section through the United States will look something like this:

Example of Cross section with topography

An additional way to plot topography is by using the 'set gxout bar' command, instead of the 'set gxout linefill' command.  The syntax is similar, but with minor differences.  You set the vrange the same, but instead of having to fill between two values, you now only need to plot the variable.  Visually, this gives your topography a more discrete, blocky look to it.

Syntax for plotting topography with 'set gxout bar':

     'set gxout bar'
     'set baropts filled'
     'set ccolor 1'
     'set lev 1000'
     'set vrange 1000 100'
     'd pressfc/100'

Resultant Image:

The real difficulty comes with rying to plot a vertical cross section on a log scale.  Since surface pressure is only a 2D variable, it is unaffected when setting the y-axis to be on a log scale.  As a result, neither of the above methods are capable of correctly plotting the topography on a log scale.  So in order to do this, we need to essentially trick GrADS into plotting the topography using a different method.

This is a little tricky to explain, but the basic idea is you simply plot the geopotential height at contour intervals so close together that they overlap and appear as one entity.  Since we are focused on topography, we subtract the surface geopotential height (elevation) from the geopotential height.  Also, since we are only interested in plotting the topography, any value >0 is masked out.  Hopefully, this short paragraph gives you some indication of what this next block of code is doing.

    'set gxout contour'
    'set clab off'
    size=750
    int=10
    i=0
   'hgt=hgtprs-hgtsfc'
  
    while(i*size>=-5000)
       j=i-1
      'set ccolor 1'
      'set cthick 10'
      'set cstyle 1'
      'set cmax 'i*size
      'set cmin 'j*size
      'set cint 'int
      'd hgt'
      i=i-1
   endwhile

Now that you can look at the code, I can explain it a little more thoroughly.  Basically, we need a while loop because GrADS cannot set enough contour levels at one time, so we have to continuously reset the contour levels until we cover the full range of values.  Within the loop we reset the line thickness and the line style to ensure continuity through the plotting.

The result of this code is a correctly plotted vertical cross section on a log scale with topography.

Vertical cross section (with terrain) plotted on a log scale

Hopefully this tutorial has provided you with a greater understanding of how to plot vertical cross sections in GrADS (at least in pressure coordinates).


Download this script

Read More

Tutorial: Shapefile Basics

The capability of GrADS to read shapefiles brings a whole new dimension to the plotting software.  It provides you with a whole new range of possibilities when it comes to plotting data.  With the use of shapefiles, you are no longer limited to generic maps with data plotted.  Shapefiles provide you the possibility of plotting everything from county boundaries, to roads and rivers, to congressional districts, weather watches and warnings, and much more.  This tutorial will give you a very basic introduction into using shapefiles in GrADS.

It is recommended (not required): that you have both color.gs and cbar.gs in your script library to plot the RADAR color scale and associated colorbar.

For this tutorial, we will simply plot county boundaries on a map of the US.  Before we get started, we will first need to download the proper shapefiles.  The county shapefiles are available from NOAA here.  This site should provide you with a table showing recent versions of the shapefile, download and unzip it.  The list of files included within the zipped folder shows four files, these are the files from the April 2013 version:

• c_02ap13.prj
• c_02ap13.shp
• c_02ap13.shx
• c_02ap13.dbf 

Most shapefiles I have encountered have some combination of these files, and because each shape tends to come with multiple files, I have found it useful to create a separate folder to store my shapefiles in, aptly named: Shapefiles.

Now that you have the county shapefiles in your proper folder we can begin with the tutorial.  The first, and possibly most important thing to know about shapefiles is that you need to plot something from your GrADS data before you can draw your shapefiles.  If you don't you will get an error that says: "No Scaling Environment."   

So, in order to plot our counties, we need to plot a variable first, which means we will need a file.  For this example we will use the RAP model using the 'sdfopen' command.

     'sdfopen http://nomads.ncep.noaa.gov:9090/dods/rap/rap20130507/rap_00z'

Now that our data file is open, the rest of the script is very short and most of it is just setting up the domain and display options.  If you have read through any previous tutorials, then next block of code will look familiar.


    'clear'
    'set lat 22 51'
    'set lon -128 -65'
    'set gxout shaded'
    'set mpdset hires'
    'set map 0'     ;*Sets map outline color to Black

With the domain set up, we just need to plot the variable.  For this example, we will plot the model predicted RADAR reflectivity.  I included in this block of code with a good approximation of the NWS color scale for reflectivity using color.gs.

 *If you have color.gs*
     'color 5 75 5 -kind white->dodgerblue->blue->lime->green->darkgreen->yellow->goldenrod->orange->firebrick->red->darkred->fuchsia->indigo'
**********************
     'd refcclm'
     'cbar'
  
Now that the RADAR is plotted, its time to draw the shapefile.  This is perhaps the easiest part of the tutorial.  You simply use the 'draw shp' command and point to the file path.  In this example:

    'draw shp Shapefiles/c_03ap12.shp'

This command is all you need to draw the counties on the map.  If you wish to change the properties of the shapefile you can do so using a few additional commands.  The outline color can be set using the 'set line command.  The fill properties can be set using the 'set shpopts' command.  By default, shapefiles are not filled.

    'set line 0'       ;*Where the outline color is set to 0
     'set shpopts 15'    ;*Fills shapefile, where 15 is the number of the desired fill color. A value of -1 turns off fill.

Depending on your exact output options, your final map should look similar to this.

RAP Reflectivity map with counties

Using shapefiles in GrADS can get much more complex.  However, for the purposes of this tutorial, I wanted to provide a real basic introduction to shapefiles in GrADS.  I may, at a later time, write up more advanced tutorials for shapefiles.  In the mean time, if you are interested in learning more about shapefiles in GrADS, I encourage you to check out this resource to get a more comprehensive look.

Again, you can download shapefiles for pretty much everything you could imagine, here is just a small list of shapefile resources to get you started.

• http://www.nws.noaa.gov/geodata/ 
• http://www.nws.noaa.gov/gis/shapepage.htm 
• http://www.statsilk.com/maps/download-free-shapefile-maps 

 
Download this Script

Read More