Skip to content

Commit

Permalink
doc/example/pre-processing updates for v2.5.0 Part1 (#477)
Browse files Browse the repository at this point in the history 
* improve tc track pre-processing script for obs

* update IM release code

* convert IBTrACS track density from 3 hourly data to 6 hourly

* preprocessing to use ibtracs

* docs updates

* add new files

* address review
  • Loading branch information
@chengzhuzhang
chengzhuzhang authored Jun 17, 2021
1 parent 91e1376 commit b036014
Show file tree
Hide file tree
Showing 4 changed files with 97 additions and 96 deletions.
12 changes: 12 additions & 0 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,3 +73,15 @@ In addition to default model versus observation comparison, the package also pro

<img src="misc/example_fig7.png" alt="Figure7" style="width: 280px;"/>
<h5 align="center">Figure 7: An example of Taylor diagram summarizing metrics calculated based on lat-lon contour plots diagnostics of several key variables</h5>

## License

Copyright (c) 2018-2021, Energy Exascale Earth System Model Project
All rights reserved

SPDX-License-Identifier: (BSD-3-Clause)

See [LICENSE](./LICENSE) for details

Unlimited Open Source - BSD 3-clause Distribution
`LLNL-CODE-819717`
10 changes: 6 additions & 4 deletions acme_diags/plot/cartopy/tc_analysis_plot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,22 +24,25 @@
border = (-0.06, -0.03, 0.13, 0.03)

plot_info = {}
# Each key gives a list with ax extent, x ticks , y ticks, title, clevs, reference
# Each key gives a list with ax extent, x ticks , y ticks, title, clevs, reference and time resolution ratio (convert 3hrly to 6hrly data, density needs to be devided by 2)
# TODO flexible to apply to 3hrly model output when compare track density.
plot_info["aew"] = [
[182, 359, 0, 35],
[240, 300],
[0, 15, 30],
"African Easterly Wave Density",
np.arange(0, 15.1, 1),
"EAR5 (2000-2014)",
1,
]
plot_info["cyclone"] = [
[0, 359, -60, 60],
[0, 60, 120, 180, 240, 300, 359.99],
[-60, -30, 0, 30, 60],
"TC Tracks Density",
np.arange(0, 1, 0.2),
np.arange(0, 0.3, 0.05),
"IBTrACS (1979-2018)",
2,
]


Expand All @@ -61,12 +64,11 @@ def plot_panel(n, fig, proj, var, var_num_years, region, title):
ax = fig.add_axes(panel[n], projection=proj)
ax.set_extent(plot_info[region][0], ccrs.PlateCarree())

clevs = np.arange(0, 15.1, 1)
clevs = plot_info[region][4]
p1 = ax.contourf(
var.getLongitude(),
var.getLatitude(),
var / var_num_years,
var / var_num_years / plot_info[region][6],
transform=ccrs.PlateCarree(),
levels=clevs,
extend="both",
Expand Down
156 changes: 66 additions & 90 deletions analysis_data_preprocess/create_global_tc_track_density.py
View file
Original file line number Diff line number Diff line change
@@ -1,107 +1,79 @@
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 20 11:06:41 2019
The observational TC data are obtained from Prof. Kerry Emanuel’s website (https://emanuel.mit.edu/products). Click “Global Tropical Cyclone Data in NETCDF format (updated through 2018)” the files will be downloaded. The TC locations are tracked 6 hourly.
@author: bala635
The observational TC data are obtained from Prof. Kerry Emanuel’s website (https://emanuel.mit.edu/products). Click “Global Tropical Cyclone Data in NETCDF format (updated through 2018)” the files will be downloaded.
Another data source is IBTrACS from NOAA, which has 3hourly track data.
Note, The tc density maps created from both data source have identical distribution, but values are twice from IBTrACS compared to MIT due to higher time frequency. For now the TC analysis has been using 6hourly E3SM output, therefore the MIT data is being used for consistency.
"""

import os
import sys
import warnings
from datetime import datetime, timedelta

if not sys.warnoptions:
warnings.simplefilter("ignore")
import numpy as np
from netCDF4 import Dataset as netcdffile

###################################################
###################################################

all_lon = []
all_lat = []

###################################################

nc = netcdffile('./attracks.nc')
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] > 1980 and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 25:

all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])

###################################################

nc = netcdffile('./eptracks.nc')
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] > 1980 and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 25:

all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])

###################################################

nc = netcdffile('./wptracks.nc')
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] > 1980 and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 25:

all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])

###################################################

nc = netcdffile('./iotracks.nc')
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] > 1980 and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 25:

all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])

###################################################

nc = netcdffile('./shtracks.nc')
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] > 1980 and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 25:

all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])

###################################################
origin_path = '/Users/zhang40/Documents/ACME/e3sm_tc_diags'
start_yr = 1979
end_yr = 2018

ibtracs = False #MIT data
ibtracs = True


if ibtracs:
data_name = 'IBTrACS'
print('Use IBTrACS 3hourly')
basins = ["NA", "WP", "EP", "NI", "SI", "SP"]
for basin in basins:
nc = netcdffile(
os.path.join(origin_path, "IBTrACS.{}.v04r00.nc".format(basin))
)
latmc = np.squeeze(nc["lat"][:, :])
longmc = np.squeeze(nc["lon"][:, :])
time = np.squeeze(nc["time"][:, :])
yearic = np.zeros((np.shape(latmc)[0]))

for i in range(0, np.shape(time)[0]):
for j in range(0, np.shape(time)[1]):
if time[i, j] is not np.ma.masked:
day_hurr = datetime(1858, 11, 17, 0, 0, 0) + timedelta(time[i, j])
yearic[i] = day_hurr.year
if yearic[i] >= start_yr and yearic[i] <= end_yr and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 0:
if(longmc[i, j]<0):
longmc[i, j] = 360 + longmc[i, j]
all_lat.append(latmc[i, j])
all_lon.append(longmc[i, j])
print(len(all_lat))


else:
print('MIT 6hrly')
data_name = 'MIT'
basins = ['at', 'ep','wp','io','sh']
for basin in basins:
nc = netcdffile('{}/{}tracks.nc'.format(origin_path, basin))
latmc = np.squeeze(nc['latmc'][:,:])
longmc = np.squeeze(nc['longmc'][:,:])
vsmc = np.squeeze(nc['vsmc'][:,:])
yearic = np.squeeze(nc['yearic'][:])

for i in range(0,np.shape(latmc)[0]):
for j in range(0,np.shape(latmc)[1]):

if yearic[j] >=start_yr and yearic[j]<=end_yr and np.abs(latmc[i,j]) > 0 and np.abs(longmc[i,j]) > 0:
all_lon.append(longmc[i,j])
all_lat.append(latmc[i,j])
print(len(all_lon))

all_lat = np.asarray(all_lat)
all_lon = np.asarray(all_lon)
Expand All @@ -111,13 +83,17 @@

import pandas as pd

raw_data = {'lon': all_lon_new,
raw_data = {
'lon': all_lon_new,
'lat': all_lat_new}

df = pd.DataFrame.from_dict(raw_data)
# columns = ['SST', 'STRAT','TCHP','TDY_SALT','TDY_NOSALT','INT','LAT','LON','SHR','DIV','PER','DELTA_INT']
df.to_csv('cyclones_all_obs.csv', sep='\t')
print(data_name)
out_file = '/Users/zhang40/Documents/ACME/e3sm_tc_diags/cyclones_hist_{}_{}_{}.csv'.format(data_name,start_yr, end_yr)
with open(out_file, 'w') as file:
file.write('start {}\n'.format(len(all_lon)))
df.to_csv(file, header=False, sep='\t')

###################################################
# Convert the .csv file to .nc by calling tempest-extremes from command line:
#tempestextremes/bin/HistogramNodes --in cyclones_all_obs.csv --iloncol 3 --ilatcol 4 --out cyclones_all_obs.nc
#tempestextremes/bin/HistogramNodes --in cyclones_hist_IBTrACS_1979_2018.csv --iloncol 2 --ilatcol 3 --out cyclones_hist_IBTrACS_1979_2018.nc
15 changes: 13 additions & 2 deletions docs/source/index.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ Algorithm and visualization codes for **latitude-longitude contour maps**,
**polar contour maps**, the accompanying **summarizing table** and **Taylor diagram plots**, **pressure-latitude zonal mean contour plots**,
**zonal mean line plots**, **pressure-longitude meridional mean contour plots**, **area mean time series plots**, and **Cloud Top Height-Tau** joint histograms
from COSP cloud simulator output. Plots can be created for annual
and seasonal climatologies, and monthly mean time series.
and seasonal climatologies, and monthly mean time series. In additional to the core sets being released in v1, **ENSO diags**, **QBO diags**, **Diurnal cycle phase plot**, **Streamflow evaluation**, **ARM diags**, and **TC analysis** are implemented in v2 release.

The package also supports custom user diagnostics, by specifying
plot type, desired region (global, ocean, land, etc.),
Expand Down Expand Up @@ -139,10 +139,21 @@ Additional back-ends could be implemented if the need arose.
| | |
| ARM diagnostics monthly diurnal cycle of cloud plot | ARM diagnostics convection onset statistics plot |
+--------------------------------------------------------+------------------------------------------------------+
| .. figure:: _static/index/fig21.png | .. figure:: _static/index/fig22.png |
| :align: center | :align: center |
| :target: _static/index/fig21.png | :target: _static/index/fig22.png |
| | |
| Tropical Cyclone Track Density | Annual Cycle Zonel Mean plot |
+--------------------------------------------------------+------------------------------------------------------+
| .. figure:: _static/index/fig23.png | .. figure:: _static/index/fig24.png |
| :align: center | :align: center |
| :target: _static/index/fig23.png | :target: _static/index/fig24.png |
| | |
| Tropical Cyclone frequency per basin | Per-basin Tropical Cyclone frac seasonal cycle |
+--------------------------------------------------------+------------------------------------------------------+

The above plots and more can be found
`here <https://portal.nersc.gov/cfs/e3sm/zhang40/tutorials/run_v230_allsets/viewer/>`_.
ARM diagnostics plots can be found `here <https://web.lcrc.anl.gov/public/e3sm/e3sm_diags_test_data/unit_test_complete_run/expected/all_sets/viewer/arm_diags/>`_.

Feature availability for each backend
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Expand Down

0 comments on commit b036014

Please sign in to comment.