Notebook

Given the altitude of a CGRF grid cell, calculate the pressure at sea level.

In [1]:

%matplotlib inline
from __future__ import division
import matplotlib.pyplot as plt
import netCDF4 as nc
import numpy as np
from salishsea_tools import (viz_tools,
                             
    bathy_tools,
    nc_tools, stormtools,
)
import requests
from xml.etree import cElementTree as ElementTree
import cStringIO

Load CGRF data¶

In [2]:

f=nc.Dataset('/ocean/dlatorne/MEOPAR/CGRF/NEMO-atmos/slp_y2006m02d04.nc')
press=f.variables['atmpres']
f=nc.Dataset('/ocean/dlatorne/MEOPAR/CGRF/NEMO-atmos/t2_y2006m02d04.nc')
T=f.variables['tair']
lat=f.variables['nav_lat']
lon=f.variables['nav_lon']

Correction provided by Frederic Dupont¶

Assuming constant air density, pressure correction is

$ p_{corr}=\rho_a g z_1 $

where $z_1$ altitude of the grid cell, $\rho_a = 1.22$ kg/m$^3$ and $g=9.80665$ m/s$^2$. Then, the the sea level pressure is

$p_{slp} = p + p_{corr}$

In [3]:

#load correction
fname='/ocean/nsoontie/MEOPAR/pcorrection.nc'
fP = nc.Dataset(fname)
corr = fP.variables['PCOR']

pfred = press[:] + corr[:]

Alternative Correction¶

Given the altitude $z_1$, pressure $p_1$ and temperature $T_1$, of an air parcel we can estimate the sea level pressure as:

$p_{slp}=p_1\left(\gamma \frac{z_1}{T_1}+1\right)^{\frac{g}{\gamma R}} $

where $g$ is the acceleration due to gravity, $R$ is the specific gas constant, and $\gamma$ is the temperature lapse rate of the atmosphere.

The assumptions are:

The temperature of the atmosphere decreases with height at a constant rate $\frac{dT}{dz}=−\gamma$.
The atmosphere is an ideal gas, $p=\rho RT$
The atmosphere is in hydrostatic equilibrium $\frac{d p}{dz}=−\rho g$.

We have taken $R = 287$ J/kg/K, $g = 9.81$ m/s$^2$, and $\gamma = 0.0098$ deg/m.

In [4]:

#correction function
R = 287 #ideal gas constant
g = 9.81 #gravity
gam = 0.0098 #lapse rate(deg/m)

def slp_correct(Z,P,T):
    ps = P*(gam*(Z/T) +1)**(g/gam/R)
    return ps

In [5]:

#load altitude
fname='/data/nsoontie/MEOPAR/tools/I_ForcingFiles/Atmos/altitude_CGRF.nc'
fP = nc.Dataset(fname)
Z = fP.variables['alt']

#correct pressure
palt = slp_correct(Z[:],press[:],T[:])

Compare with observations near YVR¶

In [6]:

#prepare for getting EC Pressure at YVR
def get_pressure_obs(station, start_day, end_day, month, year):
    station_ids = {
        'Sandheads': 6831,
        'YVR': 889,
        'PointAtkinson': 844,
        'Victoria': 10944,
        'CampbellRiver': 145,
        'PatriciaBay': 11007
    }
    wind_spd, pressure = [], []
    url = 'http://climate.weather.gc.ca/climateData/bulkdata_e.html'
    query = {
        'timeframe': 1,
        'stationID': station_ids[station],
        'format': 'xml',
        'Year': year,
        'Month': month,
        'Day': 1,
    }
    response = requests.get(url, params=query)
    tree = ElementTree.parse(cStringIO.StringIO(response.content))
    root = tree.getroot()
    raw_data = root.findall('stationdata')
    for record in raw_data:
        day = int(record.get('day'))
        hour = int(record.get('hour'))
        year = int(record.get('year'))
        month = int(record.get('month'))
        selectors = (
            (day == start_day - 1 and hour >= 16)
            or
            (day >= start_day and day < end_day)
            or
            (day == end_day and hour < 16)
        )
        if selectors:
            try:
                wind_spd.append(float(record.find('windspd').text))
            except TypeError:
                wind_spd.append(0)
            if station == 'YVR':
                try:
                    pressure.append(float(record.find('stnpress').text) )
                except ValueError:
                    pressure.append(0)
    wind_spd= np.array(wind_spd) * 1000 / 3600
    pressure = np.array(pressure)
    
    return wind_spd, pressure

In [7]:

#Load Data
[wind,press_o] = get_pressure_obs('YVR', 4, 4,2,2006) 

# CGRF Point close to YVR
point =[463, 526]


#plot
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(press_o*1000,label='observations')
ax.plot(press[:,point[0],point[1]],label='uncorrected CGRF')
ax.plot(pfred[:,point[0],point[1]],label='Freds correction')
ax.plot(palt[:,point[0],point[1]],label='alternative correction')
ax.legend(loc=0)
ax.set_xlabel('time on Feb 4 in hours'); ax.set_ylabel('Pressure [Pa]')
ax.set_title('Pressure at Feb 4, 2006 near YVR')

Out[7]:

<matplotlib.text.Text at 0x346df10>

Longer Times Series¶

In [8]:

d=15;
P=[]; PFRED=[]; PALT=[]
# CGRF Point close to YVR
point =[463, 526]
for t in np.arange(2,d+1):
    da= "%02d" % (t,)
    string='y2006m02d' +str(da) +'.nc'
    f=nc.Dataset('/ocean/nsoontie/MEOPAR/CGRF/NEMO-atmos/slp_' +string)
    press=f.variables['atmpres']; P.append(press[:,point[0],point[1]])
    pfred = press[:] +corr[:]; PFRED.append(pfred[:,point[0],point[1]])
    f=nc.Dataset('/ocean/nsoontie/MEOPAR/CGRF/NEMO-atmos/slp_corr_'+string)
    palt=f.variables['atmpres']
    PALT.append(palt[:,point[0],point[1]])
 
l = d-1    
PALT=np.array(PALT).reshape(l*24,)
P=np.array(P).reshape(l*24,)
PFRED=np.array(PFRED).reshape(l*24,)
    
[wind,press_o] = get_pressure_obs('YVR', 2, d,2,2006) 
time = np.array(stormtools.convert_date_hours(np.arange(l*24), '02-Feb-2006'))

#plot
fig,ax=plt.subplots(1,1,figsize=(10,5))
ax.plot(time,press_o*1000,label='observations')
ax.plot(time,P,label='uncorrected CGRF')
ax.plot(time,PFRED,label='Freds correction')
ax.plot(time,PALT,label='alternative correction')
ax.legend(loc=0)
ax.set_xlabel('time'); ax.set_ylabel('Pressure [Pa]')
ax.set_title('Pressure near YVR')

Out[8]:

<matplotlib.text.Text at 0x3683a90>

In [ ]: