Skip to content

Commit

Permalink
check bug in profile output
Browse files Browse the repository at this point in the history
  • Loading branch information
@lewisblake
lewisblake committed Jul 22, 2023
1 parent 7fb61b0 commit 51be00e
Show file tree
Hide file tree
Showing 2 changed files with 123 additions and 106 deletions.
226 changes: 120 additions & 106 deletions pyaerocom/aeroval/coldatatojson_engine.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
import logging
import os
from time import time
import shutil

from pyaerocom import ColocatedData, TsType
from pyaerocom._lowlevel_helpers import write_json
Expand Down Expand Up @@ -49,6 +50,12 @@ def run(self, files):
list of files that have been converted.
"""
# LB: Hacky approach to make sure writting new output
out_dirs = self.cfg.path_manager.get_json_output_dirs(True)
for idir in out_dirs:
if os.path.exists(out_dirs[idir]):
shutil.rmtree(out_dirs[idir])

converted = []
for file in files:
logger.info(f"Processing: {file}")
Expand Down Expand Up @@ -179,130 +186,137 @@ def process_coldata(self, coldata: ColocatedData):
if annual_stats_constrained:
data = _apply_annual_constraint(data) # LB: maybe this is setting everything to nans

if not diurnal_only:
logger.info("Processing statistics timeseries for all regions")
input_freq = self.cfg.statistics_opts.stats_tseries_base_freq

for reg in regnames:
try:
stats_ts = _process_statistics_timeseries(
data=data,
freq=main_freq,
region_ids={reg: regnames[reg]},
use_weights=use_weights,
use_country=use_country,
data_freq=input_freq,
if not "just_for_viz" in coldata.data.attrs: # make the regular json output
if not diurnal_only:
logger.info("Processing statistics timeseries for all regions")
input_freq = self.cfg.statistics_opts.stats_tseries_base_freq

for reg in regnames:
try:
stats_ts = _process_statistics_timeseries(
data=data,
freq=main_freq,
region_ids={reg: regnames[reg]},
use_weights=use_weights,
use_country=use_country,
data_freq=input_freq,
)

except TemporalResolutionError:
stats_ts = {}
fname = get_timeseries_file_name(
regnames[reg], obs_name, var_name_web, vert_code
)
ts_file = os.path.join(out_dirs["hm/ts"], fname)
_add_heatmap_entry_json(
ts_file, stats_ts, obs_name, var_name_web, vert_code, model_name, model_var
)

except TemporalResolutionError:
stats_ts = {}
fname = get_timeseries_file_name(regnames[reg], obs_name, var_name_web, vert_code)
ts_file = os.path.join(out_dirs["hm/ts"], fname)
_add_heatmap_entry_json(
ts_file, stats_ts, obs_name, var_name_web, vert_code, model_name, model_var
)

# breakpoint() # LB : here we need to do something for the different vertical layers.
logger.info("Processing heatmap data for all regions")

hm_all = _process_heatmap_data(
data,
regnames,
use_weights,
use_country,
meta_glob,
periods,
seasons,
add_trends,
trends_min_yrs,
)

for freq, hm_data in hm_all.items():
fname = get_heatmap_filename(freq)

hm_file = os.path.join(out_dirs["hm"], fname)
# breakpoint() # LB : here we need to do something for the different vertical layers.
logger.info("Processing heatmap data for all regions")

_add_heatmap_entry_json(
hm_file, hm_data, obs_name, var_name_web, vert_code, model_name, model_var
hm_all = _process_heatmap_data(
data,
regnames,
use_weights,
use_country,
meta_glob,
periods,
seasons,
add_trends,
trends_min_yrs,
)

logger.info("Processing regional timeseries for all regions")
ts_objs_regional = _process_regional_timeseries(data, regnames, regions_how, meta_glob)

_write_site_data(ts_objs_regional, out_dirs["ts"])
if coldata.has_latlon_dims:
for cd in data.values():
if cd is not None:
cd.data = cd.flatten_latlondim_station_name().data
for freq, hm_data in hm_all.items():
fname = get_heatmap_filename(freq)

logger.info("Processing individual site timeseries data")
(ts_objs, map_meta, site_indices) = _process_sites(data, regs, regions_how, meta_glob)
hm_file = os.path.join(out_dirs["hm"], fname)

_write_site_data(ts_objs, out_dirs["ts"])

scatter_freq = min(TsType(fq) for fq in self.cfg.time_cfg.freqs)
scatter_freq = min(scatter_freq, main_freq)
_add_heatmap_entry_json(
hm_file, hm_data, obs_name, var_name_web, vert_code, model_name, model_var
)

logger.info("Processing map and scat data by period")
for period in periods:
# compute map_data and scat_data just for this period
map_data, scat_data = _process_map_and_scat(
data,
map_meta,
site_indices,
[period],
str(scatter_freq),
stats_min_num,
seasons,
add_trends,
trends_min_yrs,
use_fairmode,
obs_var,
logger.info("Processing regional timeseries for all regions")
ts_objs_regional = _process_regional_timeseries(
data, regnames, regions_how, meta_glob
)

# the files in /map and /scat will be split up according to their time period as well
map_name = get_json_mapname(
obs_name, var_name_web, model_name, model_var, vert_code, period
_write_site_data(ts_objs_regional, out_dirs["ts"])
if coldata.has_latlon_dims:
for cd in data.values():
if cd is not None:
cd.data = cd.flatten_latlondim_station_name().data

logger.info("Processing individual site timeseries data")
(ts_objs, map_meta, site_indices) = _process_sites(
data, regs, regions_how, meta_glob
)
# breakpoint() # need format for output now. currently rewriting over previous .json files
outfile_map = os.path.join(out_dirs["map"], map_name)
write_json(map_data, outfile_map, ignore_nan=True)

outfile_scat = os.path.join(out_dirs["scat"], map_name)
write_json(scat_data, outfile_scat, ignore_nan=True)
_write_site_data(ts_objs, out_dirs["ts"])

scatter_freq = min(TsType(fq) for fq in self.cfg.time_cfg.freqs)
scatter_freq = min(scatter_freq, main_freq)

logger.info("Processing map and scat data by period")
for period in periods:
# compute map_data and scat_data just for this period
map_data, scat_data = _process_map_and_scat(
data,
map_meta,
site_indices,
[period],
str(scatter_freq),
stats_min_num,
seasons,
add_trends,
trends_min_yrs,
use_fairmode,
obs_var,
)

if coldata.ts_type == "hourly" and use_diurnal:
logger.info("Processing diurnal profiles")
(ts_objs_weekly, ts_objs_weekly_reg) = _process_sites_weekly_ts(
coldata, regions_how, regnames, meta_glob
)
outdir = os.path.join(out_dirs["ts/diurnal"])
for ts_data_weekly in ts_objs_weekly:
# writes json file
_write_stationdata_json(ts_data_weekly, outdir)
if ts_objs_weekly_reg != None:
for ts_data_weekly_reg in ts_objs_weekly_reg:
# writes json file
_write_stationdata_json(ts_data_weekly_reg, outdir)
# the files in /map and /scat will be split up according to their time period as well
map_name = get_json_mapname(
obs_name, var_name_web, model_name, model_var, vert_code, period
)
# breakpoint() # need format for output now. currently rewriting over previous .json files
outfile_map = os.path.join(out_dirs["map"], map_name)
write_json(map_data, outfile_map, ignore_nan=True)

if (
"vertical_layer" in coldata.data.attrs
): # LB: Will need some sort of additional flag to deal with the two colocation level types
logger.info("Processing profile data for vizualization")
outfile_scat = os.path.join(out_dirs["scat"], map_name)
write_json(scat_data, outfile_scat, ignore_nan=True)

for regid in regnames:
profile_viz = process_profile_data(
data,
regid,
use_country,
periods,
seasons,
if coldata.ts_type == "hourly" and use_diurnal:
logger.info("Processing diurnal profiles")
(ts_objs_weekly, ts_objs_weekly_reg) = _process_sites_weekly_ts(
coldata, regions_how, regnames, meta_glob
)
outdir = os.path.join(out_dirs["ts/diurnal"])
for ts_data_weekly in ts_objs_weekly:
# writes json file
_write_stationdata_json(ts_data_weekly, outdir)
if ts_objs_weekly_reg != None:
for ts_data_weekly_reg in ts_objs_weekly_reg:
# writes json file
_write_stationdata_json(ts_data_weekly_reg, outdir)
else:
if (
"vertical_layer" in coldata.data.attrs
): # LB: Will need some sort of additional flag to deal with the two colocation level types
logger.info("Processing profile data for vizualization")

for regid in regnames:
profile_viz = process_profile_data(
data,
regid,
use_country,
periods,
seasons,
)

fname = get_profile_filename(regid, obs_name, var_name_web)
fname = get_profile_filename(regid, obs_name, var_name_web)

outfile_profile = os.path.join(out_dirs["profile"], fname)
add_profile_entry_json(outfile_profile, data, profile_viz, periods, seasons)
outfile_profile = os.path.join(out_dirs["profile"], fname)
add_profile_entry_json(outfile_profile, data, profile_viz, periods, seasons)

# for reg in regions:
# fname = get_profile_filename(reg, obs_name, var_name_web)
Expand Down
3 changes: 3 additions & 0 deletions pyaerocom/colocation_3d.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -483,6 +483,9 @@ def colocate_vertical_profile_gridded(
# Each element in the tuple is a list of ColocatedData objects.
# The length of these lists is the same as the number of colocation layers

for coldata in output_prep[1]:
coldata.data.attrs["just_for_viz"] = 1

colocated_data_lists = ColocatedDataLists(
output_prep[0], output_prep[1]
) # put the list of prepared output into namedtuple object s.t. both position and named arguments can be used
Expand Down

0 comments on commit 51be00e

Please sign in to comment.