python tutorial

Chapter 2 Basic Python for data manipulation

2.1 Overview:

Python is an object oriented interpreted language. It has gained popularity due to its strength and ease of use. You should read the entire tutorial at python.org. That being said, I do not expect that you will. Read as much as you can and then see if you can follow the scripts below. If you can't, read more of the tutorial.

2.2 Examples:

2.2.1 Assertions

An intensive variable is derived from two (or more) extensive variables

Mass of a compound ( M O 3 ) is an extensive variable

Volume ( V ) is an extensive variable

Concentration is an intensive variable ( M s V ).

The ideal gas law: PV = nRT ; n = PV RT

Ozone concentrations are reported in parts per billion (ppb) and the 1-h standard (.12) is in parts per million (ppm)

1 ppm =1 mol Mmol = 1 0 6 mol Mmo l air × 1 mol mo l air

1 mol = 1 ppm 1 0 6 × PV RT

2.2.2 Basic Calculations

Below is a script that converts a ppm based value to total moles. Here we will take the ozone National Ambient Air Quality Standard (NAAQS) and convert it to total moles in a model grid cell.

Code listing: Single Conversion with the Ideal Gas Law

################################################ # 1. Compute moles of ozone if the NAAQ Standard # is reached in a ground level model grid cell # at 95 degrees F ################################################ print "Starting example 1" # 1 - h Average Ozone Standard (NAAQS) o3_std_ppm = .12 # Surface air pressure is 1. atmosphere pressure = 1. # Ideal gas constant in m**3 * atm / K / mol ideal_gas_constant = 82.0574587E-6 # Volume for a single grid cell in the # 2006 TCEQ Model in m**3 area = 2000.**2 height = 33.89 volume = area * height # Typical Houston temperature temp_f = 86.9 temp_c = 5. / 9. * (temp_f - 32.) temp_k = 273. + temp_c # DegreesK # Moles of Air in grid cell air_mol = pressure * volume / ideal_gas_constant / temp_k # Mega mol of Air in grid cell at each time air_Mmol = air_mol / 1E6 # Moles of O3 at 1 - h O3 standard o3_std_mol = o3_std_ppm * air_Mmol print "\t","The number of moles of ozone in a grid cell at 86.9 F is", o3_std_mol print "\tGrid cell volume km**3:", volume / 1e3**3 print "Finished example 1"

2.2.2.1 Questions:

What happens if you change 5./9. to 5/9 ? Why?

What is the ** operator?

2.2.3 Iterative Calculations

Below is a script that converts the ozone NAAQS ppm based value to total moles for all hours of a day. During the day, the temperature varies and the script identifies the minimum number of moles. The looping is all done using the “for ... in” statements to iterate over Python lists.

Code listing: Iterative Conversion with the Ideal Gas Law

################################################ # 2. For a ground level grid cell with an hourly # varying temperature find the minimum moles of # ozone and the hour at which it occurs when the # NAAQ Standard is reached in a ground level ################################################ print "Starting example 2" # 1 - h Average Ozone Standard (NAAQS) o3_std_ppm = .12 # Surface air pressure is always 1. atmosphere pressure = 1. # Ideal gas constant in m**3 * atm / K / mol ideal_gas_constant = 82.0574587E-6 # Volume for a single grid cell in the # 2006 TCEQ Model in m**3 area = 2000.**2 height = 33.89 volume = area * height # For 24 hours add a temperature hourly_temp_f = [ 78.64, 78.26, 77.92, 77.70, 77.51, 77.34, \ 77.08, 78.44, 80.80, 82.71, 84.31, 85.61, \ 86.46, 86.90, 86.50, 85.80, 84.57, 83.59, \ 82.32, 81.04, 80.42, 80.03, 79.58, 79.17] # For each hourly Fahrenheit temperature, produce a celcius # temperature hourly_temp_c = [] for this_hour_temp_f in hourly_temp_f: hourly_temp_c.append(5. / 9. * (this_hour_temp_f - 32.)) # For each hourly celcius temperature, produce a Kelvin # temperature hourly_temp_k = [] for this_hour_temp_c in hourly_temp_c: hourly_temp_k.append(this_hour_temp_c + 273.) # Moles of Air in grid cell at each time hourly_air_mol = [] for this_hour_temp_k in hourly_temp_k: hourly_air_mol.append(pressure * volume / ideal_gas_constant \ / this_hour_temp_k) # Mega mol of Air in grid cell at each time hourly_air_Mmol = [] for this_hour_air_mol in hourly_air_mol: hourly_air_Mmol.append(this_hour_air_mol / 1E6) # Moles of O3 at 1 - h O3 standard hourly_o3_std_mol = [] for this_hour_air_Mmol in hourly_air_Mmol: hourly_o3_std_mol.append(this_hour_air_Mmol * o3_std_ppm) print "\t","The minimum moles of ozone in a grid cell is:" min_o3_std_mol=None min_o3_hour=None for this_hour_temp_f, this_hour_o3_std_mol in \ zip(hourly_temp_f, hourly_o3_std_mol): if this_hour_o3_std_mol

2.2.3.1 Questions:

What if hourly_temp_c was a tuple instead of a list (i.e. change line 32 to hourly_temp_c = ())?

Could you replace the “for ... in” statements with a “while” statement?

2.2.4 Array Calculations

Below is a script that converts the ozone NAAQS ppm based value to total moles for all hours of a day. During the day, the temperature varies and the script identifies the minimum number of moles. Here we use NumPy arrays, which allows us to do math on many elements at the same time. This removes the need for iteration (looping).

Code listing: Array Conversion with the Ideal Gas Law

################################################ # 3. Compute moles of ozone if the NAAQ Standard # is reached in a ground level model grid cell # with temperature varying for all hours of the # day using a NumPy array ################################################ print "Starting example 3" from numpy import array # For more information see www.scipy.org # 1 - h Average Ozone Standard (NAAQS) o3_std_ppm = .12 # Surface air pressure is always 1. atmosphere pressure = 1. # Exercise 1: Try changing pressure to this #pressure=array([ 0.99806756, 0.99796873, 0.99786383, 0.99765235, \ # 0.99758124, 0.99771982, 0.99781328, 0.99810308, \ # 0.99881029, 0.99928504, 0.99957174, 0.9995777 , \ # 0.99918288, 0.99895877, 0.99804389, 0.99743855, \ # 0.99738914, 0.99752951, 0.99738002, 0.99742848, \ # 0.99765474, 0.99798727, 0.99814957, 0.99845177]) # Ideal gas constant in m**3 * atm / K / mol ideal_gas_constant = 82.0574587E-6 # Volume for a single grid cell in the # 2006 TCEQ Model in m**3 area = 2000.**2 height = 33.89 volume = area * height # Hourly degrees in Fahrenheit hourly_temp_f = array([ 78.64, 78.26, 77.92, 77.70, 77.51, 77.34, \ 77.08, 78.44, 80.80, 82.71, 84.31, 85.61, \ 86.46, 86.90, 86.50, 85.80, 84.57, 83.59, \ 82.32, 81.04, 80.42, 80.03, 79.58, 79.17]) #Hourly degrees in celcius hourly_temp_c = 5. / 9. * (hourly_temp_f - 32.) #Hourly degrees in Kelvin hourly_temp_k = hourly_temp_c + 273. # Moles of Air in grid cell at each time hourly_air_mol = pressure * volume / ideal_gas_constant \ / hourly_temp_k # Mega mol of Air in grid cell at each time hourly_air_Mmol = hourly_air_mol / 1E6 # Moles of O3 at 1 - h O3 standard hourly_o3_std_mol = hourly_air_Mmol * o3_std_ppm print "\t","The minimum moles of ozone in a grid cell is:" # Use the array min method to get the minimum mol value min_o3_std_mol = hourly_o3_std_mol.min() # Create an array of booleans: true where the hourly value # equals the minimum value and false where the hourly value # is anything else is_min_o3_std_mol = hourly_o3_std_mol==min_o3_std_mol # Index the hourly temperature array with the boolean array min_o3_std_temp_f = hourly_temp_f[is_min_o3_std_mol] print "\t",min_o3_std_mol,"at", min_o3_std_temp_f,"F" print "\tGrid cell volume km**3:", volume / 1e3**3 print "Finished example 3"

2.2.4.1 Questions:

min is a python intrinsic function and a method of the numpy arrays. In python, read the help on both (i.e. help(min) and from numpy import array; x=array([0]); help(x.min)). Explain how they differ.

2.2.5 Using netCDF4 to acquire arrays

In this example, we'll use real model input data to convert the NAAQS ppm based value to total moles for all hours of a day. In this example, we'll use meteorology inputs for the CMAQ model. These inputs are in NetCDF format with IOAPI conventions. We'll learn more about model data and formats in Chapter 4. This example requires that you have access to meteorology inputs for CMAQ. These inputs can be downloaded from the http://www.cmascenter.org.

    1. Create a CMAS account (if you haven't already)
    2. login to the site
    3. click on "DOWNLOAD CENTER" on the left (yes, click)
    4. click on "model clearinghouse"
    5. choose CMAQ and submit
    6. choose CMAQ 4-6 and submit
    7. click on "Download" beside "CMAQ default data files"
    8. after the download completes untar the files (tar -xzf M3DATA...)

The untarred files will be in a folder named "data." After you have downloaded and untarred the files, you can run the script. The script will prompt you for the location of a METCRO3D file. There are two METCRO3D files, on for each day (i.e. METCRO3D_YYYYJJJ), in data/mcip3/M_36_2001.

Code listing: Model Data NAAQS Conversion with the Ideal Gas Law

################################################ # 4. Compute moles of ozone if the NAAQ Standard # is reached in all grid cells for real model # data using netCDF4 to acquire NumPy arrays ################################################ print "Starting example 4" from numpy import array, where # For more information see www.scipy.org from netCDF4 import Dataset try: # This will only work on Mac. Generally, we don't # use platform specific functions. For the tutorial, # I am making an exception. If you are doing the tutorial # on a non Mac, this will fail and raw_input will be used import EasyDialogs file_path_request = lambda: EasyDialogs.AskFileForOpen( \ message = 'Select a CMAQ METCRO3D file', \ windowTitle = 'CMAQ METCRO3D') except: file_path_request = lambda: raw_input('Type the path to a CMAQ METCRO3D file:\n') # 1 - h Average Ozone Standard (NAAQS) o3_std_ppm = .12 # Get the path to a metcro_3d file metcro3d_path = file_path_request() # Open the metcro3d file metcro3d_file = Dataset(metcro3d_path,'r') print "Data used in example:" print "\tcomes from %s" % metcro3d_path print "\tstarts on %d at %06.0f" % (metcro3d_file.SDATE, metcro3d_file.STIME) # IOAPI uses GDTYP to specify the grid type (see http://www.baronams.com/products/ioapi/GRIDS.html) if metcro3d_file.GDTYP == 2: print "\tuses an LCC projection" else: print "\tuses a non-LCC projection" # Get the pressure variable from the variables dictionary pressure_4d = metcro3d_file.variables['PRES'] # Report variable properties print "Pressure has the following dimensions:", pressure_4d.dimensions print "Pressure has the following units:", pressure_4d.units # Get the air temperature variable from the variables dictionary temperature_4d = metcro3d_file.variables['TA'] # Report variable properties print "Temperature has the following dimensions:", temperature_4d.dimensions print "Temperature has the following units:", temperature_4d.units # Get the layer top height variable from the variables dictionary layer_top_height_4d = metcro3d_file.variables['ZF'] # Report variable properties print "Layer top height has the following dimensions:", \ layer_top_height_4d.dimensions print "Layer top height has the following units:", \ layer_top_height_4d.units # Ideal gas constant in m**3 * Pa / K / mol # FYI: m**3 * N / m**2 / K / mol # FYI: m**3 * kg * m / s**2 / m**2 / K / mol # FYI: m**3 * kg / m / s**2 / K / mol IDEAL_GAS_CONSTANT = 8.314472 # The size of each grid cell in the X (Y) dimension # is defined by the property XCELL (YCELL) area = metcro3d_file.XCELL * metcro3d_file.YCELL # Volume is the product of area and height volume_4d = area * layer_top_height_4d[:] # Moles of Air in grid cell at each time hourly_air_mol = pressure_4d[:] * volume_4d / IDEAL_GAS_CONSTANT \ / temperature_4d[:] # Mega mol of Air in grid cell at each time hourly_air_Mmol = hourly_air_mol / 1E6 # Moles of O3 at 1 - h O3 standard hourly_o3_std_mol = hourly_air_Mmol * o3_std_ppm print "\t","The minimum moles of ozone in a grid cell is:" # Use the array min method to get the minimum mol value min_o3_std_mol = hourly_o3_std_mol.min() # Create an array of booleans: true where the hourly value # equals the minimum value and false where the hourly value # is anything else is_min_o3_std_mol = hourly_o3_std_mol == min_o3_std_mol # Create a tuple of index arrays. Each array corresponds # to a dimensions and the values index that dimension idx_min_o3_std_mol = is_min_o3_std_mol.nonzero() # Index the 4D temperature array with the boolean array min_o3_std_temp = temperature_4d[:][is_min_o3_std_mol] # Index the 4D pressure array with the boolean array min_o3_std_pres = pressure_4d[:][is_min_o3_std_mol] # Display results print "\t",min_o3_std_mol,"at", min_o3_std_temp, temperature_4d.units, "and", min_o3_std_pres, pressure_4d.units print "\t ", temperature_4d.dimensions print "\tThe minimum occurs in cell:", idx_min_o3_std_mol print "\tGrid cell volume km**3:", volume_4d[idx_min_o3_std_mol].squeeze() / 1e3**3 print "Finished example 4"