Commit b6958ea4 authored by Kennedy, Joseph H's avatar Kennedy, Joseph H
Browse files

LEX-0.1-beta

A beta version of the LIVVkit Extensions repository to demonstrate our
extensions functionality. It lincludes the CESM Greenland analyses for
cloud cover, energy balance, and surface mass balance, as well the
CISM-Albany dynamics analyses, which will be documented in the Evans et
al. (2018) LIVVkit validation paper.

This pre-release:
  Makes the Zwally basin numbers easier to see
  Fixes rhea setup script to not cause a logout
  Fixes bad regridding in dynamics extension
  Removes the CESM2 table from energy extension
  Fixes the yearly cylce plots (were only showing out to 10 months)
  Clean up some whitespace in the energy extension
  Include the CISM-Albany dynamics extension
parents 4e99b93b e788d578
...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.', ...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.',
img_list = [] img_list = []
months = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] months = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
percent_vals = [] percent_vals = []
model_vals = [] model_vals = []
cldsat_vals = [] cldsat_vals = []
...@@ -25,17 +25,17 @@ def make_plot(config=None, out_path='.', ...@@ -25,17 +25,17 @@ def make_plot(config=None, out_path='.',
for month in months: for month in months:
# CESM1 # CESM1
f_cesm = os.path.join(cesm_path, "postproc", "atm", "climos", f_cesm = os.path.join(cesm_path, "postproc", "atm", "climos",
"b.e10.BG20TRCN.f09_g16.002_{}_aavg_climo.nc".format(month)) "b.e10.BG20TRCN.f09_g16.002_{:02d}_aavg_climo.nc".format(month))
ncid1 = Dataset(f_cesm, mode='r') ncid1 = Dataset(f_cesm, mode='r')
model_cld = ncid1.variables['CLDHGH'][0] model_cld = ncid1.variables['CLDHGH'][0]
# CLDSAT # CLDSAT
f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{}_aavg_climo.nc".format(month)) f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{:02d}_aavg_climo.nc".format(month))
ncid2 = Dataset(f_cloudsat, mode='r') ncid2 = Dataset(f_cloudsat, mode='r')
cldsat_cld = ncid2.variables['CLDHGH'][0] cldsat_cld = ncid2.variables['CLDHGH'][0]
# ISCCP # ISCCP
f_isccp = os.path.join(cloud_path, "ISCCP_{}_aavg_climo.nc".format(month)) f_isccp = os.path.join(cloud_path, "ISCCP_{:02d}_aavg_climo.nc".format(month))
ncid3 = Dataset(f_isccp, mode='r') ncid3 = Dataset(f_isccp, mode='r')
isccp_cld = ncid3.variables['CLDHGH'][0] isccp_cld = ncid3.variables['CLDHGH'][0]
...@@ -58,8 +58,10 @@ def make_plot(config=None, out_path='.', ...@@ -58,8 +58,10 @@ def make_plot(config=None, out_path='.',
plt.plot(months, isccp_vals, 'c-.') plt.plot(months, isccp_vals, 'c-.')
plt.xlabel('Months of climatology') plt.xlabel('Months of climatology')
plt.ylabel('Percent total cloud') plt.ylabel('Percent total cloud')
plt.tight_layout()
img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDHGH.png') img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDHGH.png')
plt.savefig(img_path, bbox_inches='tight') plt.savefig(img_path)
plt.close() plt.close()
img_link = os.path.join(os.path.basename(out_path), img_link = os.path.join(os.path.basename(out_path),
......
...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.', ...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.',
img_list = [] img_list = []
months = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] months = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
percent_vals = [] percent_vals = []
model_vals = [] model_vals = []
cldsat_vals = [] cldsat_vals = []
...@@ -26,17 +26,17 @@ def make_plot(config=None, out_path='.', ...@@ -26,17 +26,17 @@ def make_plot(config=None, out_path='.',
# ---------------- Data source at OLCF ------------------------ # ---------------- Data source at OLCF ------------------------
# CESM1 # CESM1
input_file1 = os.path.join(cesm_path, "postproc", "atm", "climos", input_file1 = os.path.join(cesm_path, "postproc", "atm", "climos",
"b.e10.BG20TRCN.f09_g16.002_{}_aavg_climo.nc".format(month)) "b.e10.BG20TRCN.f09_g16.002_{:02d}_aavg_climo.nc".format(month))
ncid1 = Dataset(input_file1, mode='r') ncid1 = Dataset(input_file1, mode='r')
model_cld = ncid1.variables['CLDLOW'][0] model_cld = ncid1.variables['CLDLOW'][0]
# CLDSAT # CLDSAT
f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{}_aavg_climo.nc".format(month)) f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{:02d}_aavg_climo.nc".format(month))
ncid2 = Dataset(f_cloudsat, mode='r') ncid2 = Dataset(f_cloudsat, mode='r')
cldsat_cld = ncid2.variables['CLDLOW'][0] cldsat_cld = ncid2.variables['CLDLOW'][0]
# ISCCP # ISCCP
f_isccp = os.path.join(cloud_path, "ISCCP_{}_aavg_climo.nc".format(month)) f_isccp = os.path.join(cloud_path, "ISCCP_{:02d}_aavg_climo.nc".format(month))
ncid3 = Dataset(f_isccp, mode='r') ncid3 = Dataset(f_isccp, mode='r')
isccp_cld = ncid3.variables['CLDLOW'][0] isccp_cld = ncid3.variables['CLDLOW'][0]
...@@ -59,8 +59,10 @@ def make_plot(config=None, out_path='.', ...@@ -59,8 +59,10 @@ def make_plot(config=None, out_path='.',
plt.plot(months, isccp_vals, 'c-.') plt.plot(months, isccp_vals, 'c-.')
plt.xlabel('Months of climatology') plt.xlabel('Months of climatology')
plt.ylabel('Percent total cloud') plt.ylabel('Percent total cloud')
plt.tight_layout()
img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDLOW.png') img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDLOW.png')
plt.savefig(img_path, bbox_inches='tight') plt.savefig(img_path)
plt.close() plt.close()
img_link = os.path.join(os.path.basename(out_path), img_link = os.path.join(os.path.basename(out_path),
......
...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.', ...@@ -16,7 +16,7 @@ def make_plot(config=None, out_path='.',
img_list = [] img_list = []
months = ['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'] months = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]
percent_vals = [] percent_vals = []
model_vals = [] model_vals = []
cldsat_vals = [] cldsat_vals = []
...@@ -26,17 +26,17 @@ def make_plot(config=None, out_path='.', ...@@ -26,17 +26,17 @@ def make_plot(config=None, out_path='.',
# ---------------- Data source at OLCF ------------------------ # ---------------- Data source at OLCF ------------------------
# CESM1 # CESM1
f_cesm = os.path.join(cesm_path, "postproc", "atm", "climos", f_cesm = os.path.join(cesm_path, "postproc", "atm", "climos",
"b.e10.BG20TRCN.f09_g16.002_{}_aavg_climo.nc".format(month)) "b.e10.BG20TRCN.f09_g16.002_{:02d}_aavg_climo.nc".format(month))
ncid1 = Dataset(f_cesm, mode='r') ncid1 = Dataset(f_cesm, mode='r')
model_cld = ncid1.variables['CLDTOT'][0] model_cld = ncid1.variables['CLDTOT'][0]
# CLDSAT # CLDSAT
f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{}_aavg_climo.nc".format(month)) f_cloudsat = os.path.join(cloud_path, "CLOUDSAT_{:02d}_aavg_climo.nc".format(month))
ncid2 = Dataset(f_cloudsat, mode='r') ncid2 = Dataset(f_cloudsat, mode='r')
cldsat_cld = ncid2.variables['CLDTOT'][0] cldsat_cld = ncid2.variables['CLDTOT'][0]
# ISCCP # ISCCP
f_isccp = os.path.join(cloud_path, "ISCCP_{}_aavg_climo.nc".format(month)) f_isccp = os.path.join(cloud_path, "ISCCP_{:02d}_aavg_climo.nc".format(month))
ncid3 = Dataset(f_isccp, mode='r') ncid3 = Dataset(f_isccp, mode='r')
isccp_cld = ncid3.variables['CLDTOT'][0] isccp_cld = ncid3.variables['CLDTOT'][0]
...@@ -59,8 +59,10 @@ def make_plot(config=None, out_path='.', ...@@ -59,8 +59,10 @@ def make_plot(config=None, out_path='.',
plt.plot(months, isccp_vals, 'c-.') plt.plot(months, isccp_vals, 'c-.')
plt.xlabel('Months of climatology') plt.xlabel('Months of climatology')
plt.ylabel('Percent total cloud') plt.ylabel('Percent total cloud')
plt.tight_layout()
img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDTOT.png') img_path = os.path.join(out_path, 'CESM_yearly_cycle_CLDTOT.png')
plt.savefig(img_path, bbox_inches='tight') plt.savefig(img_path)
plt.close() plt.close()
img_link = os.path.join(os.path.basename(out_path), img_link = os.path.join(os.path.basename(out_path),
......
import os
import Ngl
import numpy.ma as ma
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """CISMA_js plot"""
# plots velnorm, zooms in on Jakobshavn 69.2N 51.1W
def make_plot(config=None, out_path='.',
cism_path='/lustre/atlas1/cli115/world-shared/4ue/'):
# ---------------- Read Data -----------------------------------
input_file1 = os.path.join(cism_path, 'Gris_init/greenland_1km_2016_12_01.mcb_velnorm.nc')
input_file4 = os.path.join(cism_path, 'CISM-Albany/postproc/albany.out.tstep.ANN_velnorm.nc')
input_file5 = os.path.join(cism_path, 'CISM-Albany/postproc/albany.out.tstep.ANN_thk.nc')
input_file6 = os.path.join(cism_path, 'CISM-grids/climo_greenland_1km_2017_05_06_notime.nc')
# --------------------------------------------------------------
img_list = []
# read input source, x1, y1
ncid1 = Dataset(input_file1, 'r')
init = ncid1.variables['velnorm'][0, :, :]
# read lat and lon for original CISMA, x1 y1
ncid6 = Dataset(input_file6, 'r')
lat1 = ncid6.variables['lat'][:]
lon1 = ncid6.variables['lon'][:]
# read CISMA data output x0, y0
ncid4 = Dataset(input_file4, 'r')
cisma = ncid4.variables['velnorm'][0, 0, :, :]
# read lat and lon for CISMA output, velnorm x0, y0
# read simulated thk of CISM x1, y1
ncid5 = Dataset(input_file5, 'r')
thk = ncid5.variables['thk'][0, :, :]
# read lat and lon for CISMA output, thk x1 y1
thk_mask = ma.masked_less(thk, 0.001)
vel_mask = ma.masked_less(cisma, 0.001)
cisma = ma.masked_array(cisma, mask=vel_mask.mask)
init = ma.masked_array(init, mask=thk_mask.mask)
init = init[:-1, :-1]
diff = cisma - init
# ------- PLOT --------
# Open a workstation for drawing the plots
wkres = Ngl.Resources()
wkres.wkColorMap = "WhiteBlueGreenYellowRed"
# wkres.wkOrientation = "portrait" # "portrait" or "landscape"
wks_type = "png"
wks_img = str(os.path.join(out_path, "CISMA_js"))
wks = Ngl.open_wks(wks_type, wks_img, wkres)
# --- for the map -------
# Define plotting area, Greenland
mres = Ngl.Resources()
mres.nglDraw = False # Don't draw individual plots
mres.nglFrame = False # Don't advance frame.
mres.pmTickMarkDisplayMode = "Never" # Turn off map tickmarks.
mres.mpGridAndLimbOn = False # Turn off grid and limb lines.
mres.mpProjection = "Aitoff"
mres.mpLimitMode = "LatLon" # limit map via lat/lon, to zoom in
mres.mpCenterLatF = 72. # map area
mres.mpCenterLonF = -48.
mres.mpMinLatF = 68.
mres.mpMaxLatF = 78.
mres.mpMinLonF = -55.
mres.mpMaxLonF = -42.
mres.mpOutlineOn = False
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the file1 contour -------
res1 = Ngl.Resources()
res1.nglDraw = False # Don't draw individual plots
res1.nglFrame = False # Don't advance frame.
res1.cnLineLabelsOn = False
res1.cnFillOn = True
res1.cnLinesOn = False
res1.cnLineLabelsOn = False
res1.cnFillMode = "RasterFill"
res1.trGridType = "TriangularMesh"
res1.cnLevelSelectionMode = "ExplicitLevels"
res1.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res1.lbLabelBarOn = True # Turn on labelbar.
res1.lbOrientation = "Horizontal" # Verticle labelbar
res1.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res1.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res1.sfXArray = lon1
res1.sfYArray = lat1
# --- for the file2 contour -------
res2 = Ngl.Resources()
res2.nglDraw = False # Don't draw individual plots
res2.nglFrame = False # Don't advance frame.
res2.cnLineLabelsOn = False
res2.cnFillOn = True
res2.cnLinesOn = False
res2.cnLineLabelsOn = False
res2.cnFillMode = "RasterFill"
res2.trGridType = "TriangularMesh"
res2.cnLevelSelectionMode = "ExplicitLevels"
res2.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res2.lbLabelBarOn = True # Turn on labelbar.
res2.lbOrientation = "Horizontal" # Verticle labelbar
res2.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res2.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res2.sfXArray = lon1
res2.sfYArray = lat1
# --- for diff=file1-file2 contour -------
res22 = Ngl.Resources()
res22.nglDraw = False # Don't draw individual plots
res22.nglFrame = False # Don't advance frame.
res22.cnLineLabelsOn = False
res22.cnFillOn = True
res22.cnLinesOn = False
res22.cnLineLabelsOn = False
res22.cnFillMode = "RasterFill"
res22.trGridType = "TriangularMesh"
res22.cnLevelSelectionMode = "ExplicitLevels"
res22.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res22.lbLabelBarOn = True # Turn on labelbar.
res22.lbOrientation = "Horizontal" # Verticle labelbar
res22.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res22.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res22.sfXArray = lon1
res22.sfYArray = lat1
init_plot = Ngl.contour(wks, init, res1)
cisma_plot = Ngl.contour(wks, cisma, res2)
diff_plot = Ngl.contour(wks, diff, res22)
# Creat multiple figures and draw, which now contains the elevation and temperature
# "[1,3]" indicates 1 row, 3 columns.
map_title = ["Init", "CISM-A", "CISM-A - Init"]
nmap = 3
plot = []
for i in range(nmap):
mres.tiMainString = map_title[i]
plot.append(Ngl.map(wks, mres))
# Overlay everything on the map plot.
Ngl.overlay(plot[0], cisma_plot)
Ngl.overlay(plot[1], init_plot)
Ngl.overlay(plot[2], diff_plot)
Ngl.panel(wks, plot, [1, 3])
img_link = os.path.join(os.path.basename(out_path),
os.path.basename(wks_img + '.' + wks_type))
img_elem = el.image('CISMA_js',
' '.join(describe.split()),
img_link)
if config:
img_elem['Height'] = config['image_height']
img_list.append(img_elem)
return img_list
if __name__ == '__main__':
make_plot()
import os
import Ngl
import numpy.ma as ma
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """CISMA_pm plot"""
# Plots velnorm, zooms in on Petermann glacier, 80.5N -59.5
def make_plot(config=None, out_path='.',
cism_path='/lustre/atlas1/cli115/world-shared/4ue/'):
# ---------------- Read Data -----------------------------------
input_file1 = os.path.join(cism_path, 'Gris_init/greenland_1km_2016_12_01.mcb_velnorm.nc')
input_file4 = os.path.join(cism_path, 'CISM-Albany/postproc/albany.out.tstep.ANN_velnorm.nc')
input_file5 = os.path.join(cism_path, 'CISM-Albany/postproc/albany.out.tstep.ANN_thk.nc')
input_file6 = os.path.join(cism_path, 'CISM-grids/climo_greenland_1km_2017_05_06_notime.nc')
# input_file7 = '/lustre/atlas1/cli115/world-shared/4ue/CISM-grids/zurich_mask.nc' # x1 y1 mask for simulated CISMA
# --------------------------------------------------------------
img_list = []
# read input source, x1, y1
ncid1 = Dataset(input_file1, 'r')
init = ncid1.variables['velnorm'][0, :, :]
# read lat and lon for original CISMA, x1 y1
ncid6 = Dataset(input_file6, 'r')
lat1 = ncid6.variables['lat'][:]
lon1 = ncid6.variables['lon'][:]
# read CISMA data output x0, y0
ncid4 = Dataset(input_file4, 'r')
cisma = ncid4.variables['velnorm'][0, 0, :, :]
# read lat and lon for CISMA output, velnorm x0, y0
# read simulated thk of CISM x1, y1
ncid5 = Dataset(input_file5, 'r')
thk = ncid5.variables['thk'][0, :, :]
# read lat and lon for CISMA output, thk x1 y1
thk_mask = ma.masked_less(thk, 0.001)
vel_mask = ma.masked_less(cisma, 0.001)
cisma = ma.masked_array(cisma, mask=vel_mask.mask)
init = ma.masked_array(init, mask=thk_mask.mask)
init = init[:-1, :-1]
diff = cisma - init
# ------- PLOT --------
# Open a workstation for drawing the plots
wkres = Ngl.Resources()
wkres.wkColorMap = "WhiteBlueGreenYellowRed"
# wkres.wkOrientation = "portrait" # "portrait" or "landscape"
wks_type = "png"
wks_img = str(os.path.join(out_path, "CISMA_pm"))
wks = Ngl.open_wks(wks_type, wks_img, wkres)
# --- for the map -------
# Define plotting area, Greenland
mres = Ngl.Resources()
mres.nglDraw = False # Don't draw individual plots
mres.nglFrame = False # Don't advance frame.
mres.pmTickMarkDisplayMode = "Never" # Turn off map tickmarks.
mres.mpGridAndLimbOn = False # Turn off grid and limb lines.
mres.mpProjection = "Aitoff"
mres.mpLimitMode = "LatLon" # limit map via lat/lon, to zoom in
mres.mpCenterLatF = 80. # map area
mres.mpCenterLonF = -48.
mres.mpMinLatF = 76.
mres.mpMaxLatF = 84.
mres.mpMinLonF = -55.
mres.mpMaxLonF = -42.
mres.mpOutlineOn = False
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the file1 contour -------
res1 = Ngl.Resources()
res1.nglDraw = False # Don't draw individual plots
res1.nglFrame = False # Don't advance frame.
res1.cnLineLabelsOn = False
res1.cnFillOn = True
res1.cnLinesOn = False
res1.cnLineLabelsOn = False
res1.cnFillMode = "RasterFill"
res1.trGridType = "TriangularMesh"
res1.cnLevelSelectionMode = "ExplicitLevels"
res1.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res1.lbLabelBarOn = True # Turn on labelbar.
res1.lbOrientation = "Horizontal" # Verticle labelbar
res1.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res1.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res1.sfXArray = lon1
res1.sfYArray = lat1
# --- for the file2 contour -------
res2 = Ngl.Resources()
res2.nglDraw = False # Don't draw individual plots
res2.nglFrame = False # Don't advance frame.
res2.cnLineLabelsOn = False
res2.cnFillOn = True
res2.cnLinesOn = False
res2.cnLineLabelsOn = False
res2.cnFillMode = "RasterFill"
res2.trGridType = "TriangularMesh"
res2.cnLevelSelectionMode = "ExplicitLevels"
res2.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res2.lbLabelBarOn = True # Turn on labelbar.
res2.lbOrientation = "Horizontal" # Verticle labelbar
res2.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res2.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res2.sfXArray = lon1
res2.sfYArray = lat1
# --- for diff=file1-file2 contour -------
res22 = Ngl.Resources()
res22.nglDraw = False # Don't draw individual plots
res22.nglFrame = False # Don't advance frame.
res22.cnLineLabelsOn = False
res22.cnFillOn = True
res22.cnLinesOn = False
res22.cnLineLabelsOn = False
res22.cnFillMode = "RasterFill"
res22.trGridType = "TriangularMesh"
res22.cnLevelSelectionMode = "ExplicitLevels"
res22.cnLevels = (0, 1, 2, 10, 30, 70, 100, 200, 500, 1000, 2000)
res22.lbLabelBarOn = True # Turn on labelbar.
res22.lbOrientation = "Horizontal" # Verticle labelbar
res22.lbLabelFontHeightF = 0.02 # Make fonts smaller.
res22.pmLabelBarOrthogonalPosF = -0.04 # move label bar closer
res22.sfXArray = lon1
res22.sfYArray = lat1
init_plot = Ngl.contour(wks, init, res1)
cisma_plot = Ngl.contour(wks, cisma, res2)
diff_plot = Ngl.contour(wks, diff, res22)
# Creat multiple figures and draw, which now contains the elevation and temperature
# "[1,3]" indicates 1 row, 3 columns.
map_title = ["Init", "CISM-A", "CISM-A - Init"]
nmap = 3
plot = []
for i in range(nmap):
mres.tiMainString = map_title[i]
plot.append(Ngl.map(wks, mres))
# Overlay everything on the map plot.
Ngl.overlay(plot[0], cisma_plot)
Ngl.overlay(plot[1], init_plot)
Ngl.overlay(plot[2], diff_plot)
Ngl.panel(wks, plot, [1, 3])
img_link = os.path.join(os.path.basename(out_path),
os.path.basename(wks_img + '.' + wks_type))
img_elem = el.image('CISMA_pm',
' '.join(describe.split()),
img_link)
if config:
img_elem['Height'] = config['image_height']
img_list.append(img_elem)
return img_list
if __name__ == '__main__':
make_plot()
import os
import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """scatter_CISMA_RACMO23 plot."""
def make_plot(config=None, out_path='.',
racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/',
cism_path='/lustre/atlas1/cli115/world-shared/4ue/'):
img_list = []
# --------------------------------------------------------------
# first use Remap_CISM2RACMO.sh to get remapped file and source data file
# which were saved in my local directory NCLplot
# ---------------- Read Data -----------------------------------
input_file1 = os.path.join(racmo_path, 'racmo23.ANN_with_latlon_smb.nc')
input_file2 = os.path.join(cism_path, 'CISM-Albany/postproc/remapped_cism/CISM-Albany.acab.remap2racmo.ANN.nc')
input_file3 = os.path.join(racmo_path, 'RACMO23_masks_ZGRN11.nc')
# --------------------------------------------------------------
# input_file1 get the following, lat(rlat,rlon),lat=312,lon=306
# smb(0,:,:), smb(time, lat, lon) with unit"mmWE", have some values= 0.005072316
# input_file2 get y1(y1,x1), y1=312,x1=306, here lat_name is y1
# acab(0,:,:), acab(time, y1, x1) with unit"meter/year" have some values = 0
# Goal: plot scatter figure of surface mass balance between RACMO and CISMA
# read racmo source data
ncid1 = Dataset(input_file1, mode='r')
racmo = ncid1.variables['smb'][0, :, :]
ncid1.close()