This page demonstrates Python tips and tricks that I use in my everyday programming as an atmospheric science graduate student.
-Brian Blaylock

Tuesday, August 30, 2016

Dealing with python datetime arrays

Sometimes it's a pain that you can't perform datetime operations on each element of a numpy array. Here are a couple tipes

1. List of a range of datetime objects

The best way to make a list of datetime objects or numpy array is to use this tip:


which creates a list of datetime objects from today going back a day until it hits the number of days range limit. You can change the - to a + to increment one day in the future, or change the days argument to hours, etc.

2. Convert list of datetime strings to python datetime objects

When you read in date and time data from a file, it will often come in as a string. Often, you'll read these in as a list. To convert that list of strings to a list of python datetime objects do this:

3. Convert list of epoch times to python datetime objects


Alternatively, if the dates come in as a list of epoch times, use this...


4. Convert list of datetime strings to python datetime objects using a numpy vectorized function

This method is slightly appears to be slightly faster than method discussed in 2.


In the past I would create an empty array and append to it with a for loop like this...
but this method takes longer to compute, especially for long arrays.

Monday, August 29, 2016

Decoding Rawinsonde Sounding Observations

Have you ever seen upper air sounding data that looks like this...I'm trying to write a Python script to decode this kind of stuff. I'm not sure what to make of this or how the best way to decode it is.
There are some hints here: http://www.ofcm.gov/fmh3/pdf/13-app-e.pdf
and here: http://www.meteor.iastate.edu/classes/mt311/extras/Codul-TEMP.pdf
and here: http://www.rwic.und.edu/metresources/raob.html

USUS41 KWBC 191400

IIAA USBOI02 69144 99441 71160 15645 09333

99910 12464 00000 00144 ///// ///// 92796 ///// ///// 85536

18876 01008 70162 10476 34521 50588 07170 33030 40756 20769

33043

88999

77999

31313 40708 81429=

NNNN

UMUS41 KWBC 191400

IIBB USBOI02 69148 99441 71160 15645 09333

00910 12464 11870 18073 22813 18477 33808 18278 44727 10874

55694 10476 66590 00468 77501 06969 88387 22969

21212 11910 00000 22840 01010 33795 05010 44712 36019 55712

36019 66641 31518 77449 34038 88387 33041

31313 40708 81429

41414 /////

21212 11910 00000 22840 01010 33795 05010 44712 36019 55712

36019 66641 31518 77449 34038 88387 33041=

NNNN

Tuesday, August 2, 2016

Find the nearest latitude and longitude grid box for a point

🐍 An updated notebook is available for this demonstration on GitHub.

Also look at my demonstration using the KDTree method (scipy.spatial.KDTree.query) to find nearest neighbor for many points.



I deal a lot with weather model output (specifically HRRR model data) and often need to verify a data point with observational data. Weather stations we get observations from are not located at the exact grid points as the model's grid points. The hard part is finding the index value of a latitude and longitude grid point nearest to the observation location. This is made easy with numpy's maximum function (not max).

First, you need to have three arrays, one for latitude, one for longitude, and one for the environment variable such as temperature. Load these values into a numpy array. I use pygrib to read model output in grib2 format.

lat = latitude  # a 2D numpy array of your latitudes
lon = longitude # a 2D numpy array of your longitudes
temp = temperature # a 2D numpy array of your temperatures, or other variable

Next you need to know the latitude and longitude for the observation point. We are looking for the nearest grid point in the lat and lon arrays for that grid point.

stn_lat = station's latitude # eg. 40
stn_lon = station's longitude # eg. -110

Now find the absolute value of the difference between the  station's lat/lon with every point in the grid. This tells us how close a point is to the particular latitude and longitude.

abslat = np.abs(lat-stn_lat)
abslon= np.abs(lon-stn_lon)

If you plot abslat and abslon, you'll see what we have done...we have created two new arrays with that tell us how close each latitude and longitude is to the grid point we are interested in.
abslat: the blue line represents the nearest grid latitude line our station is located.
abslon: the blue line represents the nearest grid longitude our station is located.
Now we need to combine these two results.
We will use numpy.maximum, which takes two arrays and finds the local maximum.

c = np.maximum(abslon,abslat)

Which is an array that looks like this...
c: an array expressing the local maximums between the abslat and abslon arrays
The dark blue point is the nearest grid location to our station. You can find the index location on the grid of this by using the min function.

latlon_idx = np.argmin(c)

Now, this latitude/longitude index value is the index for a flattened array, so when you look for that same index value in, say, your temperature array, you should flatten the array to pluck out the value at that point.

grid_temp = temp.flat[latlon_idx]

If you don't like flattened arrays, you can also get the row/column index like this
x, y = np.where(c == np.min(c))
grid_temp = temp[x[0], y[0]

And there you have it!

Python Matplotlib Superscript

In scientific notation it is more proper to express fractional units with superscipts like ms^-1 instead of m/s.  To do this in labels in python's Matplotlib there is some special formating you need to do...

plt.ylabel(r'Wind Speed (ms$\mathregular{^{-1}}$)')

this will print "Wind Speed (ms^-1)" with the -1 as a superscript on the ylabel. The \mathregular expression makes the font not a fancy italic font. See here: http://stackoverflow.com/questions/21226868/superscript-in-python-plots