Commit 1fd9ca28 authored by Vladimir Onoprienko's avatar Vladimir Onoprienko

Initial commit

parents
#Copmare observations meteo data and model output
This project provides tools for operating and comparing observed meteo conditions and gridded fields of meteo quantities produced by an atmospheric model.
Tools are based on [NCAR Command language](https://www.ncl.ucar.edu) (ncl). The source data are converted into the netCDF format to be operated.
## Data sources
### List of WMO stations {#wmo-list}
The official list of meteo stations provided by WMO is [here](https://www.wmo.int/cpdb/volume_a_observing_stations/list_stations). It is hard to be combined into a single file, so the tool uses the list downloaded here [here](ftp://ftp-cdc.dwd.de/pub/CDC/help/stations_list_CLIMAT_data.txt). It is not ultimately full and may become outdated. If so you may edit it manually.
### Meteo observations
Currently the only supported format of meteo observations is isd-lite. You can download data the from [NCAR database](https://www1.ncdc.noaa.gov/pub/data/noaa/isd-lite). The files have format XXXXX0-99999-YYYY where XXXXX is WMO id of a stations and YYYY is a year of reported observations.
### Atmospheric models
Currently the only model which output is supported is [Arctic System Reanalysis NCAR](https://rda.ucar.edu/datasets/ds631.0/) of both 15km and 30km resolutions. You can download data there. Register If you have no NCAR account. For the script to operate it needs at least full one-year set of model data.
## Dependencies
You need UNIX and ncl-ncarg package (v6.4.0) installed .
## How to use
This section gives step-by-step instructions how to compare observations and model results using tools provided here. There are two ways to operate the scripts: using command line arguments or editing the configuring script *workspace/config.ncl*
1. Update the file *data/wmo_ids.txt* with list of WMO stations if needed. See [wmo list](#wmo-list) section
2. Enter working directory
cd workspace
3. Convert wmo stations list *data/wmo_ids.txt* from text format to a netCDF format file *netcdf/workspace/wmo_ids.nc*.
ncl wmoid_to_netcdf.ncl
To change default input and output files edit values of `wmo_ids_list_txt`, `wmo_ids_list_nc` attributes in configurating script *workspace/config.ncl* or pass command line arguments when running the script:
ncl wmoid_to_netcdf.ncl 'wmo_ids_list_txt="../data/wmo_ids.txt"' 'wmo_ids_list_nc="./netcdf/wmo_ids.nc"'
4. Download *isd-lite* into *data/isd*. Convert them to a netCDF format. List wmo identifiers in the attribute `wmo_ids` of configuring script or pass as cmd option. For setting periods of time list years with 'years'. If you need only one station `wmo_id` cmd argument is recognized, for one year only you may pass `year` argument.
ncl isd_to_obs.ncl 'wmo_ids=(/ "23724", "01025", "71916"/)' year=2004
This converts files *../data/isd/XXXXX0-99999-2004* to *./workspace/netcdf/obs_XXXXX-2004.nc* adding station metadata from *./netcdf/wmo_ids.nc*. You may change defaults using `wmo_ids_list_nc`, `isd_dir`, and `obs_nc_prefix` options.
5. Interpolate ASR gridded data to station observations. Pass `isd_nc_dir` to set directory where isd files are searched, to set stations pass `wmo_id` or `wmo_ids`, to set years pass `year` or `years`, to set directory where asr file are searched pass `asr_dir`, pass `asr_prefix`: source asr files are *{asr_dir}/asr_prefix-YYYYMMDD.nc*, pass `output_prefix` to produce output files {output_prefix}_XXXXX-YYYY
ncl asr_to_obs.ncl 'isd_nc_dir="./netcdf"' 'asr_dir="~/workspace/data/asr30"' 'asr_prefix="asr30km.anl.2D."' 'wmo_ids=(/"23724", "01025", "71916"/)' year=2004 'output_prefix="netcdf/asr30"'
5. To plot observed values and corresponding model values run at one graph run
ncl compare.ncl
## For developers
- To print sea level pressure provided by reanalisys run
ncl plot.ncl
- *station.ncl* seems to be not working
- *check_re.ncl* is not ready yet
This diff is collapsed.
;Interpolate ASR NCAR gridded model data to meteostation
; load configuring script
if(.not.isdefined("config_src_ncl")) then
config_src_ncl = "./config.ncl"
end if
loadscript(config_src_ncl)
begin
;set global constants
msg_read_cmd = "Checking input command line variables: obs_nc_prefix, asr_dir, asr_basename, wmo_id, wmo_ids, year, years, model_obs_prefix"
msg_print_init = "Following values are set:"
nstation_dimlen = 1
Celsius_zero_in_Kelvin = 273.15
linmsg_opt = -1; to extapolate missing values at the end the intervals
;get configuration
if(.not.isdefined("use_config")) then
use_config = 1
end if
config = get_default_config(use_config)
; check cmd variables
print(""+msg_read_cmd)
if ( .not.isdefined("obs_nc_prefix") ) then
obs_nc_prefix = config@obs_nc_prefix
end if
if ( .not.isdefined("asr_dir") ) then
asr_dir = config@asr_dir
end if
if ( .not.isdefined("asr_basename") ) then
asr_basename = config@asr_basename
end if
if ( (.not.isdefined("wmo_id")).and.(.not.isdefined("wmo_ids")) ) then
wmo_ids = config@wmo_ids
else
if ( (isdefined("wmo_id")).and.(.not.isdefined("wmo_ids")) ) then
wmo_ids = (/tostring(wmo_id)/)
end if
end if
if ( (.not.isdefined("year")).and.(.not.isdefined("years")) ) then
years = config@years
else
if ( (isdefined("year")).and.(.not.isdefined("years")) ) then
years = (/ year /)
end if
end if
if ( .not.isdefined("model_obs_prefix") ) then
model_obs_prefix = config@model_obs_prefix
end if
print(""+msg_print_init)
print("obs_nc_prefix : " + obs_nc_prefix)
print("asr_dir : " + asr_dir)
print("asr_basename : " + asr_basename)
print("years : " + years)
print("wmo_ids : " + wmo_ids)
print("model_obs_prefix: " + model_obs_prefix)
;get number of stations and number of years
num_ids = dimsizes(wmo_ids)
num_years = dimsizes(years)
;main loop
do year_num = 0, num_years - 1
year = stringtointeger(years(year_num))
time_zero = 0.0
time_zero@units = "days since " + year + "-01-01 00:00:0.0"
time_current = time_zero
if(isleapyear(year)) then
num_days = 366
else
num_days = 365
end if
model_files = new((/num_days/), string)
;;set names of model files
do day_num = 0, num_days-1
time_current = time_zero + day_num
date_str = cd_calendar(time_current, 2)
filename = asr_basename + date_str + ".nc"
model_files(day_num) = asr_dir + "/" + filename
end do
f_model = addfiles(model_files, "r")
ListSetType (f_model, "cat")
time_model_read = f_model[:]->Time
lat2d = f_model[0]->XLAT
lon2d = f_model[0]->XLONG
do id_num = 0, num_ids-1
isd_nc_file = obs_nc_prefix + wmo_ids(id_num) + "-" + year + ".nc"
output_file = model_obs_prefix + wmo_ids(id_num) + "-" + year + ".nc"
f = addfile(isd_nc_file,"r")
time = f->time
lat = (/ f->latitude /)
lon = (/ f->longitude /)
alt = (/ f->altitude /)
id_obs = (/ f->wmo_id /)
pressure_obs = f->psl
temperature_obs = f->t
print("year, id, lat, lon:" + " " + year(year_num) + " " + wmo_ids(id_num) + " " + lat + " " + lon)
nm = getind_latlon2d (lat2d, lon2d, lat, lon)
print(" lat, lon of neighbour model grid cell and station : (" + lat2d(nm(0,0), nm(0,1)) + ", " + lon2d(nm(0,0), nm(0,1)) + "), (" + lat(0) + ", " + lon(0) + ")")
time_model = cd_convert(time_model_read, time@units)
;;get model pressure
pressure_model_read = f_model[:]->PMSL(:, nm(0,0), nm(0,1))
pressure_model_read = pressure_model_read * 0.01
pressure_model = new((/dimsizes(time)/), float)
pressure_model@long_name = "Model Sea-level Pressure interpolated to station observations"
pressure_model@standard_name = "pressure"
pressure_model@units = "hPa"
pressure_model@source_prefix = asr_basename
pressure_model = linint1(time_model, pressure_model_read, False, time, 0)
pressure_model = linmsg(pressure_model, linmsg_opt)
pressure_model!0 = "nrecord"
;;get model temperature
temperature_model_read = (/ f_model[:]->T2M(:, nm(0,0), nm(0,1)) /)
temperature_model_read = temperature_model_read - Celsius_zero_in_Kelvin
temperature_model = new((/dimsizes(time)/), float)
temperature_model_read@long_name = "Model air temperature at 2m height intepolated to station observation"
temperature_model_read@standard_name = "temperature"
temperature_model_read@units = "Celsius"
temperature_model_read@source_prefix = asr_basename
temperature_model = linint1(time_model, temperature_model_read, False, time, 0)
temperature_model = linmsg(temperature_model, linmsg_opt)
temperature_model!0 = "nrecord"
;;write output netCDF file
setfileoption("nc","preFill",False)
setfileoption("nc","defineMode",True)
system("rm -f " + output_file)
fout = addfile(output_file,"c")
;system("cp " + isd_nc_file + " " + output_file )
dim_names = getvardims(f)
dim_sizes = getfiledimsizes(f)
dim_unlimited = conform (dim_names,(/False/),0)
filedimdef(fout,dim_names,dim_sizes,dim_unlimited)
filevardef(fout, "latitude" , typeof(lat) , "nstation")
filevardef(fout, "longitude", typeof(lon) , "nstation")
filevardef(fout, "altitude" , typeof(alt) , "nstation")
filevardef(fout, "wmo_id" , typeof(id_obs), "ncharacter")
filevardef(fout, "time" , typeof(time) , "nrecord")
filevardef(fout, "psl" , typeof(pressure_model) , "nrecord")
filevardef(fout, "t" , typeof(temperature_model), "nrecord")
filevardef(fout, "psl_obs", typeof(pressure_obs) , "nrecord")
filevardef(fout, "t_obs" , typeof(temperature_obs) , "nrecord")
filevarattdef(fout,"psl" , pressure_model_read)
filevarattdef(fout,"t" , temperature_model_read)
;filevarattdef(fout,"psl_obs", pressure_obs)
;filevarattdef(fout,"t_obs" , temperature_obs)
lat!0 = "nstation"
lon!0 = "nstation"
alt!0 = "nstation"
lat!0 = "nstation"
lat!0 = "nstation"
id_obs!0 = "ncharacter"
pressure_model!0 = "nrecord"
temperature_model!0 = "nrecord"
pressure_obs!0 = "nrecord"
temperature_obs!0 = "nrecord"
fout->latitude = lat
fout->longitude = lon
fout->altitude = alt
fout->wmo_id = id_obs
fout->time = time
fout->psl = pressure_model
fout->t = temperature_model
fout->psl_obs = pressure_obs
fout->t_obs = temperature_obs
delete(fout)
delete(time)
delete(pressure_obs)
delete(temperature_obs)
delete(pressure_model)
delete(pressure_model_read)
delete(temperature_model)
delete(temperature_model_read)
end do
end do
delete(f)
delete(f_model)
end
; Convert isd data to a netCDF format
; load configuring script
if(.not.isdefined("config_src_ncl")) then
config_src_ncl = "./config.ncl"
end if
loadscript(config_src_ncl)
begin
;set global constants
wmo_id_len = 5
num_isd_pressure_field = 7
num_isd_temperature_field = 5
nstation_dimlen = 1
msg_read_cmd = "Checking input command line variables: wmo_ids_list_nc, wmo_id, wmo_ids, year, years, isd_dir, obs_nc_prefix"
msg_print_init = "Following values are set:"
;get configuration
if(.not.isdefined("use_config")) then
use_config = 1
end if
config = get_default_config(use_config)
;check cmd variables
print(""+msg_read_cmd)
if ( .not.isdefined("wmo_ids_list_nc") ) then
wmo_ids_list_nc = config@wmo_ids_list_nc
end if
if ( (.not.isdefined("wmo_id")).and.(.not.isdefined("wmo_ids")) ) then
wmo_ids = config@wmo_ids
else
if ( (isdefined("wmo_id")).and.(.not.isdefined("wmo_ids")) ) then
wmo_ids = (/tostring(wmo_id)/)
end if
end if
if ( (.not.isdefined("year")).and.(.not.isdefined("years")) ) then
years = config@years
else
if ( (isdefined("year")).and.(.not.isdefined("years")) ) then
years = (/ year /)
end if
end if
if ( .not.isdefined("isd_dir") ) then
isd_dir = config@isd_dir
end if
if ( .not.isdefined("obs_nc_prefix") ) then
obs_nc_prefix = config@obs_nc_prefix
end if
print(""+msg_print_init)
print("wmo_ids_list_nc :" + wmo_ids_list_nc)
print("wmo_ids :" + wmo_ids)
print("years :" + years)
print("isd_dir : " + isd_dir)
print("obs_nc_prefix : " + obs_nc_prefix)
; get nuber of stations and years; set an array
num_ids = dimsizes(wmo_ids)
num_years = dimsizes(years)
wmoid_index = new((/num_ids/), integer)
;read list of wmo stations
f1 = addfile(wmo_ids_list_nc,"r")
ids_wmo = tostring(f1->id(:,0:wmo_id_len-1))
lat_wmo = f1->latitude
lon_wmo = f1->longitude
alt_wmo = f1->altitude
delete(f1)
num_stations_wmo = dimsizes(lat_wmo)
;find given stations in the list
do id_num = 0, num_ids-1
do station_num = 0, num_stations_wmo-1
if (wmo_ids(id_num).eq.ids_wmo(station_num)) then
wmoid_index(id_num) = station_num
print("station " + wmo_ids(id_num) + " found, number: " + wmoid_index(id_num))
break
end if
end do
end do
;read metadata of inerested stations in the list
lat = new((/num_ids/), float)
lon = new((/num_ids/), float)
alt = new((/num_ids/), float)
id = new((/num_ids/), string)
copy_VarAtts(ids_wmo, id)
copy_VarAtts(lat_wmo, lat)
copy_VarAtts(lon_wmo, lon)
copy_VarAtts(alt_wmo, alt)
lat!0 = "nstation"
lon!0 = "nstation"
alt!0 = "nstation"
do id_num = 0, num_ids-1
id(id_num) = ids_wmo(wmoid_index(id_num))
lat(id_num) = tofloat(lat_wmo(wmoid_index(id_num)))
lon(id_num) = tofloat(lon_wmo(wmoid_index(id_num)))
alt(id_num) = tofloat(alt_wmo(wmoid_index(id_num)))
print("id, lat, lon, alt: " + id(id_num) + " " + lat(id_num) + " " + lon(id_num) + " " + alt(id_num))
end do
;read isd file and convert to netCDF
do id_num = 0, num_ids-1
do year_num = 0, num_years - 1
year = years(year_num)
filename = wmo_ids(id_num) + "0-99999-" + year
file_isd = isd_dir + "/" + filename
print("file_isd: " + file_isd)
nrow = numAsciiRow(file_isd)
ncol = numAsciiCol(file_isd)
isd_data = asciiread(file_isd, (/nrow, ncol/) , "integer")
year_arr = isd_data(:,0)
month = isd_data(:,1)
day = isd_data(:,2)
hour = isd_data(:,3)
minute = hour*0
second = hour*0
time_units = "days since " + year + "-01-01 00:00:0.0"
opt = 0
opt@return_type = "float"
opt@calendar = "proleptic_gregorian"
time = cd_inv_calendar(year_arr, month, day, hour,minute,second,time_units, opt)
time!0 = "nrecord"
time@standard_name = "time"
time@long_name = "Time of Observation at Station"
;read pressure
pressure = tofloat(isd_data(:,num_isd_pressure_field-1))
pressure = pressure*0.1
pressure!0 = "nrecord"
pressure@standard_name = "pressure"
pressure@long_name = "Observed Sea Level Pressure at screen height (~2m)"
pressure@units = "hPa"
temperature = tofloat(isd_data(:,num_isd_temperature_field-1))
temperature = temperature*0.1
temperature!0 = "nrecord"
temperature@standard_name = "temperature"
temperature@long_name = "Observed temperature of the air"
temperature@units = "Celsius"
print("pressure :" + pressure(1:3))
print("temperature :" + temperature(1:3))
out_filename = wmo_ids(id_num) + "-" + year + ".nc"
out_name = obs_nc_prefix + out_filename
print("out_name: " + out_name)
setfileoption("nc","preFill",False)
setfileoption("nc","defineMode",True)
system("rm -f " + out_name)
fout = addfile(out_name,"c")
id_char = stringtochar(id(id_num))
dim_names = (/"nrecord", "nstation" , "ncharacter"/)
dim_sizes = (/nrow , nstation_dimlen, wmo_id_len /)
dim_unlimited = (/False , False , False /)
filedimdef(fout,dim_names,dim_sizes,dim_unlimited)
filevardef(fout, "latitude" , typeof(lat), "nstation")
filevardef(fout, "longitude", typeof(lon), "nstation")
filevardef(fout, "altitude" , typeof(alt), "nstation")
filevardef(fout, "time", typeof(time) , "nrecord")
filevardef(fout, "psl" , typeof(pressure) , "nrecord")
filevardef(fout, "t" , typeof(temperature), "nrecord")
lat_out = lat(id_num)
lon_out = lon(id_num)
alt_out = alt(id_num)
lat_out!0 = "nstation"
lon_out!0 = "nstation"
alt_out!0 = "nstation"
id_out = id_char(0:wmo_id_len-1)
id_out!0 = "ncharacter"
fout->latitude = lat_out
fout->longitude = lon_out
fout->altitude = alt_out
fout->wmo_id = id_out
fout->time = time
fout->psl = pressure
fout->t = temperature
delete(fout)
delete(isd_data)
delete(pressure)
delete(temperature)
delete(time)
delete(year_arr)
delete(month)
delete(day)
delete(hour)
delete(minute)
delete(second)
delete(id_char)
delete(id_out)
end do
end do
end
;Convert wmo ids to a netCDF format
; load configuring script
if(.not.isdefined("config_src_ncl")) then
config_src_ncl = "./config.ncl"
end if
loadscript(config_src_ncl)
begin
;set global constants
config_var_names = (/"wmo_ids_list_txt", "wmo_ids_list_nc"/)
msg_input_vars = "Checking input variables: "
msg_print_init = "Following values are set: "
;get configuration
if(.not.isdefined("use_config")) then
use_config = 1
end if
config = get_default_config(use_config)
;print list of input variables
print("" +msg_input_vars)
print(" "+config_var_names)
; check cmd variables
if ( .not.isdefined("wmo_ids_list_txt").and.(use_config.gt.0) ) then
wmo_ids_list_txt = config@wmo_ids_list_txt
end if
if ( .not.isdefined("wmo_ids_list_nc").and.(use_config.gt.0) ) then
wmo_ids_list_nc = config@wmo_ids_list_nc
end if
if (any(.not.isdefined(config_var_names) )) then
print("ERROR: undefined input variable")
print(" Following input variables are not defined, but screept need them.")
print(" " + mask(config_var_names, .not.isdefined(config_var_names), 1))
print("Try to edit config file or define variable via command line.")
print("EXIT.")
exit()
end if
print(""+msg_print_init)
print(" wmo_ids_list_txt: " + wmo_ids_list_txt)
print(" wmo_ids_list_nc : " + wmo_ids_list_nc)
; read wmo ids text file
nrow = numAsciiRow(wmo_ids_list_txt)
ncol = numAsciiCol(wmo_ids_list_txt)
print("nrow, ncol: " + nrow + ", " + ncol)
data = asciiread(wmo_ids_list_txt, (/nrow/) , "string")
fields = str_split_csv(data, ";", 0)
row_index = new((/nrow-1/), integer)
rec_num = 0
mask_arr = .not.ismissing(fields(:, :)).and.(strlen(str_strip(fields(:,:))).gt.0); check if bad field
mask_arr(0,:) = False
do row_num = 1, nrow-1
if ( all(mask_arr(row_num,:)) ) then
row_index(rec_num) = row_num
rec_num = rec_num + 1
end if
end do
num_rec = rec_num
print("num_rec: " + num_rec)
id = new((/num_rec/), string)
name = new((/num_rec/), string)
lat = new((/num_rec/), float)
lon = new((/num_rec/), float)
alt = new((/num_rec/), float)
country = new((/num_rec/), string)
do rec_num = 0, num_rec-1
id(rec_num) = str_strip(fields(row_index(rec_num), 0))
name(rec_num) = str_strip(fields(row_index(rec_num), 1))
lat(rec_num) = stringtofloat(fields(row_index(rec_num), 2))
lon(rec_num) = stringtofloat(fields(row_index(rec_num), 3))
alt(rec_num) = stringtofloat(fields(row_index(rec_num), 4))
country(rec_num) = str_strip(fields(row_index(rec_num), 5))
end do