Home
cosmo_query is a set of routines to automatically download COSMO data from the CSCS servers for specific regions, dates, resolutions (COSMO-1, COSMO-2 or COSMO-7) and simulation modes (analysis or forecast) as well as to extract data profiles at fixed altitude levels, fixed latitude, longitude and to interpolate COSMO data along radar RHI and PPI scans.
These routines require only few dependencies that you should most likely already have installed on your system (numpy, scipy, pyproj, netCDF4, h5py and paramiko)
You can simply clone the library by using
git clone [email protected]:wolfensb/cosmo_query.git
Alternatively you can download the compressed folder under the Files tab of the GitHub project.
Then cd to the downloaded folder, make sure that you are in the same directory as setup.py and install the library with
sudo python3 setup.py install
or
sudo python setup.py install
depending if you are using python2.7 or python3.
All dependencies should be installed automatically. In case they are still missing you will need to install them manually
sudo pip install netCDF4 numpy scipy pyproj
General imports
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineSpecific imports
from cosmo_query import config, SSH, Query
from cosmo_query import save_netcdf, load_netcdf, save_hdf5, load_hdf5
from cosmo_query import extract, coords_profileFirst step is to establish a connection to CSCS servers, this connection can then be used to retrieved COSMO data
We initiate a connection to ela.cscs.ch, with username and password specified in the config.py file (password is not needed if ssh key is defined)
connection = SSH(config.ELA_ADRESS,config.USERNAME,config.PASSWORD)Connecting to server...
We also need to open a channel to kesch.cscs.ch since field_extra (the tool that is used to extract and convert from Grib files) is not available on ela.cscs.ch
connection.open_channel(config.KESCH_ADRESS,config.USERNAME,config.PASSWORD)Now that the connection is setup we can create a Query instance
query = Query(connection)And use this query to retrieve some data
variables = ['T','P'] # we want temperature and pressure, see next cell to know which variables you can get
date = '2016-05-31 12:30' # for the 31th May 2016 at 12h30
model_res = 'fine' # at high resolution (COSMO-2 in this case)
mode = 'analysis' # In analysis mode
coord_bounds = ([6.6,45.8],[8.4,46.6]) # Over an area covering roughly the Valais
data = query.retrieve_data(variables, date, model_res = 'fine',
mode = 'analysis', coord_bounds = coord_bounds)Output will be written to '/users/wolfensb/PMSL.grb'.
Extracting table of content...
Filtering...
Output will be written to '/users/wolfensb/PMSL_cosmo_2.nc'.
Extracting table of content...
Converting to NETCDF...
(163, 233, 136, 187)
Output will be written to '/users/wolfensb/filtered0.grb'.
Extracting table of content...
Filtering...
Output will be written to '/users/wolfensb/crop0.grb'.
Extracting table of content...
Cropping...
Output will be written to '/users/wolfensb/convert0.nc'.
Extracting table of content...
Converting to NETCDF...
Output will be written to '/users/wolfensb/filtered1.grb'.
Extracting table of content...
Filtering...
Output will be written to '/users/wolfensb/crop1.grb'.
Extracting table of content...
Cropping...
Output will be written to '/users/wolfensb/convert1.nc'.
Extracting table of content...
Converting to NETCDF...
To know which variables are in the file, you can either check the file logs in your "path_to_library"/cosmo_info/file_logs/ folder or you can use the function get_available_vars
print(query.get_available_vars(date = '2016-05-31 12:30', model_res = 'fine', mode = 'analysis'))We can save this data to a netcdf
save_netcdf(data,'myfile.nc')And load it again from this file
data = load_netcdf('myfile.nc')We can then use this data to extract a profile, here are some examples
Here we get the temperature at 2000 and 4000 m
prof1 = extract(data,['T'],'level',[2000,4000])/usr/local/lib/python2.7/dist-packages/COSMO_query/extract.py:139: RuntimeWarning: invalid value encountered in divide
var[:][indexes_smaller[i],j,k])/dist*(cv-zlevels[indexes_smaller[i],j,k])
The output is a dictionary, with one key 'T' (which corresponds to the variables you wanted), the variable prof1['T'] is a special numpy arrays with one additional class variables coordinates (an ordered dictionary containing all coordinates, ordered based on the dimension of the numpy array, in this case heights,lon, lat).
prof1['T'].coordinates.keys()['heights', 'lat', 'lon']
Here we get the profile of pressure and temperature at 46.2 degrees latitude
prof2 = extract(data,['P','T'],'lat',46.2)The output is a dictionary, with two keys 'T' and 'P' (which corresponds to the variables you wanted). In this case the coordinates are:
prof2['T'].coordinates.keys()['heights', 'lon/lat']
Where heights is an
Here we get the profile of pressure at 46.2 degrees longitude
prof3 = extract(data,['T'],'lon',7)Structure of output is very similar to longitudinal profile
We can also get a transect (vert. profile) through an arbitrary set of lon/lat coordinates. If you want a rectilinear profile you can use the coords_profile function
# Transect going from Martigny to Visp with a step of 2000 m
coords = coords_profile([7.07,46.1],[7.88,46.29],step=2000)
# or transect going from Martigny to Visp with 100 equidistant points
coords2 = coords_profile([7.07,46.1],[7.88,46.29],npts = 100)Output is a
If wou want a non-rectilinear profile, you will have to construct your
prof3 = extract(data,['T'],'lonlat',coords)Coordinates are the same as for lat and lon profiles. Besides coordinates the outputs arrays have an additional variable called attributes which contains additional useful info. In this case there is only one attribute distance which gives the distance along the profile.
It is also possible to simulate PPI profiles from a radar, to do this you need to provide a dictionary with the necessary parameters (mandatory ones are in bold)
- beamwidth : 3dB beamwidth of the radar (in degrees)
- elevation : elevation angle of the PPI (in degrees)
- rrange : 1D array containing the distance of all radar gates
- rpos : position of the radar in the form [lon,lat,altitude]
- azimuth : 1D array containing all azimuth angles in the scan. Default is from 0 to 360 with a step corresponding to the 3dB beamwidth.
- npts_quad : tuple of the form [npts_hor, npts_ver], the interpolation scheme uses a quadrature scheme to interpolate over the antenna diagram. This defines the number of integration points in horizontal and vertical directions. Use [1,1] if you want no integration. Default is [3,3]
- refraction_method : To compute the propagation of the radar beam, two schemes can be used: if refraction_scheme = 1, the standard 4/3 Earth radius scheme will be used, if refraction_scheme = 2, a more accurate scheme based on a refractivity profile estimated from the model will be used; this requires QV (water vapour), T (temperature) and P (pressure) to be present in the file. Default is 1
In the example below we simulate a PPI scan at 10 degree elevation in Martigny
options_PPI = {}
options_PPI['beamwidth'] = 1.5 # 1.5 deg 3dB beamwidth, as MXPol
options_PPI['elevation'] = 10
options_PPI['rrange'] = np.arange(200,20000,75) # from 0.2 to 20 km with a res of 75 m
options_PPI['npts_quad'] = [3,3] # 3 quadrature points in azimuthal, 3 in elevational directions
options_PPI['rpos'] = [7.0923,46.1134,500] # Radar position in lon/lat/alt
options_PPI['refraction_method'] = 1 # 1 = standard 4/3, 2 = ODE refraction
ppi_T = extract(data,['T'],'PPI',options_PPI) Could not read 'azimuth' field in input dict.
Using default one: np.arange(0,360,beamwidth_3dB)
Could not read 'refraction_method' field in input dict.
Using default one: 1 (4/3 Earth radius method)
Processing azimuth ang. 0.0
Processing azimuth ang. 1.5
Processing azimuth ang. 3.0
Processing azimuth ang. 4.5
Processing azimuth ang. 6.0
Processing azimuth ang. 7.5
Processing azimuth ang. 9.0
Processing azimuth ang. 10.5
Processing azimuth ang. 12.0
Processing azimuth ang. 13.5
Processing azimuth ang. 15.0
Processing azimuth ang. 16.5
Processing azimuth ang. 18.0
Processing azimuth ang. 19.5
Processing azimuth ang. 21.0
Processing azimuth ang. 22.5
Processing azimuth ang. 24.0
Processing azimuth ang. 25.5
Processing azimuth ang. 27.0
Processing azimuth ang. 28.5
Processing azimuth ang. 30.0
Processing azimuth ang. 31.5
Processing azimuth ang. 33.0
Processing azimuth ang. 34.5
Processing azimuth ang. 36.0
Processing azimuth ang. 37.5
Processing azimuth ang. 39.0
Processing azimuth ang. 40.5
Processing azimuth ang. 42.0
Processing azimuth ang. 43.5
Processing azimuth ang. 45.0
Processing azimuth ang. 46.5
Processing azimuth ang. 48.0
Processing azimuth ang. 49.5
Processing azimuth ang. 51.0
Processing azimuth ang. 52.5
Processing azimuth ang. 54.0
Processing azimuth ang. 55.5
Processing azimuth ang. 57.0
Processing azimuth ang. 58.5
Processing azimuth ang. 60.0
Processing azimuth ang. 61.5
Processing azimuth ang. 63.0
Processing azimuth ang. 64.5
Processing azimuth ang. 66.0
Processing azimuth ang. 67.5
Processing azimuth ang. 69.0
Processing azimuth ang. 70.5
Processing azimuth ang. 72.0
Processing azimuth ang. 73.5
Processing azimuth ang. 75.0
Processing azimuth ang. 76.5
Processing azimuth ang. 78.0
Processing azimuth ang. 79.5
Processing azimuth ang. 81.0
Processing azimuth ang. 82.5
Processing azimuth ang. 84.0
Processing azimuth ang. 85.5
Processing azimuth ang. 87.0
Processing azimuth ang. 88.5
Processing azimuth ang. 90.0
Processing azimuth ang. 91.5
Processing azimuth ang. 93.0
Processing azimuth ang. 94.5
Processing azimuth ang. 96.0
Processing azimuth ang. 97.5
Processing azimuth ang. 99.0
Processing azimuth ang. 100.5
Processing azimuth ang. 102.0
Processing azimuth ang. 103.5
Processing azimuth ang. 105.0
Processing azimuth ang. 106.5
Processing azimuth ang. 108.0
Processing azimuth ang. 109.5
Processing azimuth ang. 111.0
Processing azimuth ang. 112.5
Processing azimuth ang. 114.0
Processing azimuth ang. 115.5
Processing azimuth ang. 117.0
Processing azimuth ang. 118.5
Processing azimuth ang. 120.0
Processing azimuth ang. 121.5
Processing azimuth ang. 123.0
Processing azimuth ang. 124.5
Processing azimuth ang. 126.0
Processing azimuth ang. 127.5
Processing azimuth ang. 129.0
Processing azimuth ang. 130.5
Processing azimuth ang. 132.0
Processing azimuth ang. 133.5
Processing azimuth ang. 135.0
Processing azimuth ang. 136.5
Processing azimuth ang. 138.0
Processing azimuth ang. 139.5
Processing azimuth ang. 141.0
Processing azimuth ang. 142.5
Processing azimuth ang. 144.0
Processing azimuth ang. 145.5
Processing azimuth ang. 147.0
Processing azimuth ang. 148.5
Processing azimuth ang. 150.0
Processing azimuth ang. 151.5
Processing azimuth ang. 153.0
Processing azimuth ang. 154.5
Processing azimuth ang. 156.0
Processing azimuth ang. 157.5
Processing azimuth ang. 159.0
Processing azimuth ang. 160.5
Processing azimuth ang. 162.0
Processing azimuth ang. 163.5
Processing azimuth ang. 165.0
Processing azimuth ang. 166.5
Processing azimuth ang. 168.0
Processing azimuth ang. 169.5
Processing azimuth ang. 171.0
Processing azimuth ang. 172.5
Processing azimuth ang. 174.0
Processing azimuth ang. 175.5
Processing azimuth ang. 177.0
Processing azimuth ang. 178.5
Processing azimuth ang. 180.0
Processing azimuth ang. 181.5
Processing azimuth ang. 183.0
Processing azimuth ang. 184.5
Processing azimuth ang. 186.0
Processing azimuth ang. 187.5
Processing azimuth ang. 189.0
Processing azimuth ang. 190.5
Processing azimuth ang. 192.0
Processing azimuth ang. 193.5
Processing azimuth ang. 195.0
Processing azimuth ang. 196.5
Processing azimuth ang. 198.0
Processing azimuth ang. 199.5
Processing azimuth ang. 201.0
Processing azimuth ang. 202.5
Processing azimuth ang. 204.0
Processing azimuth ang. 205.5
Processing azimuth ang. 207.0
Processing azimuth ang. 208.5
Processing azimuth ang. 210.0
Processing azimuth ang. 211.5
Processing azimuth ang. 213.0
Processing azimuth ang. 214.5
Processing azimuth ang. 216.0
Processing azimuth ang. 217.5
Processing azimuth ang. 219.0
Processing azimuth ang. 220.5
Processing azimuth ang. 222.0
Processing azimuth ang. 223.5
Processing azimuth ang. 225.0
Processing azimuth ang. 226.5
Processing azimuth ang. 228.0
Processing azimuth ang. 229.5
Processing azimuth ang. 231.0
Processing azimuth ang. 232.5
Processing azimuth ang. 234.0
Processing azimuth ang. 235.5
Processing azimuth ang. 237.0
Processing azimuth ang. 238.5
Processing azimuth ang. 240.0
Processing azimuth ang. 241.5
Processing azimuth ang. 243.0
Processing azimuth ang. 244.5
Processing azimuth ang. 246.0
Processing azimuth ang. 247.5
Processing azimuth ang. 249.0
Processing azimuth ang. 250.5
Processing azimuth ang. 252.0
Processing azimuth ang. 253.5
Processing azimuth ang. 255.0
Processing azimuth ang. 256.5
Processing azimuth ang. 258.0
Processing azimuth ang. 259.5
Processing azimuth ang. 261.0
Processing azimuth ang. 262.5
Processing azimuth ang. 264.0
Processing azimuth ang. 265.5
Processing azimuth ang. 267.0
Processing azimuth ang. 268.5
Processing azimuth ang. 270.0
Processing azimuth ang. 271.5
Processing azimuth ang. 273.0
Processing azimuth ang. 274.5
Processing azimuth ang. 276.0
Processing azimuth ang. 277.5
Processing azimuth ang. 279.0
Processing azimuth ang. 280.5
Processing azimuth ang. 282.0
Processing azimuth ang. 283.5
Processing azimuth ang. 285.0
Processing azimuth ang. 286.5
Processing azimuth ang. 288.0
Processing azimuth ang. 289.5
Processing azimuth ang. 291.0
Processing azimuth ang. 292.5
Processing azimuth ang. 294.0
Processing azimuth ang. 295.5
Processing azimuth ang. 297.0
Processing azimuth ang. 298.5
Processing azimuth ang. 300.0
Processing azimuth ang. 301.5
Processing azimuth ang. 303.0
Processing azimuth ang. 304.5
Processing azimuth ang. 306.0
Processing azimuth ang. 307.5
Processing azimuth ang. 309.0
Processing azimuth ang. 310.5
Processing azimuth ang. 312.0
Processing azimuth ang. 313.5
Processing azimuth ang. 315.0
Processing azimuth ang. 316.5
Processing azimuth ang. 318.0
Processing azimuth ang. 319.5
Processing azimuth ang. 321.0
Processing azimuth ang. 322.5
Processing azimuth ang. 324.0
Processing azimuth ang. 325.5
Processing azimuth ang. 327.0
Processing azimuth ang. 328.5
Processing azimuth ang. 330.0
Processing azimuth ang. 331.5
Processing azimuth ang. 333.0
Processing azimuth ang. 334.5
Processing azimuth ang. 336.0
Processing azimuth ang. 337.5
Processing azimuth ang. 339.0
Processing azimuth ang. 340.5
Processing azimuth ang. 342.0
Processing azimuth ang. 343.5
Processing azimuth ang. 345.0
Processing azimuth ang. 346.5
Processing azimuth ang. 348.0
Processing azimuth ang. 349.5
Processing azimuth ang. 351.0
Processing azimuth ang. 352.5
Processing azimuth ang. 354.0
Processing azimuth ang. 355.5
Processing azimuth ang. 357.0
Processing azimuth ang. 358.5
Processing azimuth ang. 360.0
Output is in polar coordinates with coordinates = rrange, azimuth
print(ppi_T['T'].coordinates.keys())['rrange', 'azimuth']
Its attributes are the height of the radar gates (altitude), their lon and lat, the elevation angle as well as the frac_power of every radar gate which indicates which fraction of the incident power reached the radar gate. 0 means that all the signal was lost (i.e. complete beam-blockage), 1 means that no signal was lost.
print(ppi_T['T'].attributes.keys())['heights', 'lon', 'lat', 'elevation', 'frac_power']
plt.contourf(ppi_T['T'].coordinates['azimuth'],ppi_T['T'].coordinates['rrange'],ppi_T['T'],
levels = np.arange(240,295,1),extend='both')
plt.colorbar(label = 'Temp. [K]')
plt.xlabel('Azimuth [degree]')
plt.ylabel('Range [m]')<matplotlib.text.Text at 0x7f85765628d0>
Missing areas (white pixels) correspond to areas with total beam blockage
plt.contourf(ppi_T['T'].coordinates['azimuth'],ppi_T['T'].coordinates['rrange'],
ppi_T['T'].attributes['frac_power'],levels = np.arange(0,1,0.1),extend='both')
plt.colorbar(label = 'Fraction of power')
plt.xlabel('Elevation [degree]')
plt.ylabel('Range [m]')<matplotlib.text.Text at 0x7f8576e35a50>
Note that the COSMO DEM is quite imprecise, so the topography of COSMO might be higher than the real topography...for example real altitude of this point is 458.2 but closest COSMO point is at 626.28 m altitude...
It is also possible to simulate RHI profiles from a radar. The outputs are exactly the same as for PPI except that the elevation argument is now an optional list (with a default of 0 to 90 with a step corresponding to the 3dB beamwidth) and that the azimuth argument is a mandatory scalar.
Below we simulate a RHI scan from Martigny with an azimuth angle of 47 degrees corresponding to the direction of the valley.
options_RHI = {}
options_RHI['beamwidth'] = 1.5 # 1.5 deg 3dB beamwidth, as MXPol
options_RHI['azimuth'] = 47
options_RHI['rrange'] = np.arange(200,10000,75) # from 0.2 to 10 km with a res of 75 m
options_RHI['npts_quad'] = [3,3] # 3 quadrature points in azimuthal, 3 in elevational directions
options_RHI['rpos'] = [7.0923,46.1134,500] # Radar position in lon/lat/alt
options_RHI['refraction_method'] = 1 # 1 = standard 4/3, 2 = ODE refraction
rhi_T = extract(data,['T'],'RHI',options_RHI) Could not read 'elevation' field in input dict.
Using default one: np.arange(0,90,beamwidth_3dB)
Processing elevations ang. 0.0
Processing elevations ang. 1.5
Processing elevations ang. 3.0
Processing elevations ang. 4.5
Processing elevations ang. 6.0
Processing elevations ang. 7.5
Processing elevations ang. 9.0
Processing elevations ang. 10.5
Processing elevations ang. 12.0
Processing elevations ang. 13.5
Processing elevations ang. 15.0
Processing elevations ang. 16.5
Processing elevations ang. 18.0
Processing elevations ang. 19.5
Processing elevations ang. 21.0
Processing elevations ang. 22.5
Processing elevations ang. 24.0
Processing elevations ang. 25.5
Processing elevations ang. 27.0
Processing elevations ang. 28.5
Processing elevations ang. 30.0
Processing elevations ang. 31.5
Processing elevations ang. 33.0
Processing elevations ang. 34.5
Processing elevations ang. 36.0
Processing elevations ang. 37.5
Processing elevations ang. 39.0
Processing elevations ang. 40.5
Processing elevations ang. 42.0
Processing elevations ang. 43.5
Processing elevations ang. 45.0
Processing elevations ang. 46.5
Processing elevations ang. 48.0
Processing elevations ang. 49.5
Processing elevations ang. 51.0
Processing elevations ang. 52.5
Processing elevations ang. 54.0
Processing elevations ang. 55.5
Processing elevations ang. 57.0
Processing elevations ang. 58.5
Processing elevations ang. 60.0
Processing elevations ang. 61.5
Processing elevations ang. 63.0
Processing elevations ang. 64.5
Processing elevations ang. 66.0
Processing elevations ang. 67.5
Processing elevations ang. 69.0
Processing elevations ang. 70.5
Processing elevations ang. 72.0
Processing elevations ang. 73.5
Processing elevations ang. 75.0
Processing elevations ang. 76.5
Processing elevations ang. 78.0
Processing elevations ang. 79.5
Processing elevations ang. 81.0
Processing elevations ang. 82.5
Processing elevations ang. 84.0
Processing elevations ang. 85.5
Processing elevations ang. 87.0
Processing elevations ang. 88.5
Processing elevations ang. 90.0
plt.contourf(rhi_T['T'].coordinates['elevation'],rhi_T['T'].coordinates['rrange'],rhi_T['T'],
levels = np.arange(240,295,1),extend='both')
plt.colorbar(label = 'Temp. [K]')
plt.xlabel('Elevation [degree]')
plt.ylabel('Range [m]')<matplotlib.text.Text at 0x7f8576c51750>
It is possible to save extracted profiles to the disk and to load them with the save_hdf5 and the load_hdf5 functions
save_hdf5(rhi_T,'rhi_profile.hdf5')
rhi_T = load_hdf5('rhi_profile.hdf5')