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

Energy and clouds now work with both CESM and E3SM, smb still CESM only

parent 5837470d
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
......@@ -26,7 +26,7 @@
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
"""An analysis of CESM's energy balance over Greenland."""
"""An analysis of the Models's energy balance over Greenland."""
from __future__ import absolute_import, division, print_function, unicode_literals
......@@ -40,21 +40,21 @@ from livvkit.util import functions as fn
with fn.temp_sys_path(os.path.dirname(__file__)):
import energy.utils as utils
import energy.model_racmo23_albedo as cesm_albedo
import energy.model_racmo23_latf as cesm_latf
import energy.model_racmo23_lwsd as cesm_lwsd
import energy.model_racmo23_lwsn as cesm_lwsn
import energy.model_racmo23_rnet as cesm_rnet
import energy.model_racmo23_senf as cesm_senf
import energy.model_racmo23_smb as cesm_smb
import energy.model_racmo23_swsd as cesm_swsd
import energy.model_racmo23_swsn as cesm_swsn
import energy.model_racmo23_albedo as albedo
import energy.model_racmo23_latf as latf
import energy.model_racmo23_lwsd as lwsd
import energy.model_racmo23_lwsn as lwsn
import energy.model_racmo23_rnet as rnet
import energy.model_racmo23_senf as senf
import energy.model_racmo23_smb as smb
import energy.model_racmo23_swsd as swsd
import energy.model_racmo23_swsn as swsn
import energy.model_racmo23_t2m_ann as t2m_ann
import energy.model_racmo23_t2m_djf as t2m_djf
import energy.model_racmo23_t2m_jja as t2m_jja
import energy.rad_aavg_model as cesm_rad
import energy.rad_aavg_model as model_rad
import energy.rad_aavg_racmo as racmo_rad
import energy.timeseries_trefht as ts_trefht
import energy.timeseries_tsa as ts_tsa
import energy.timeseries_qice_gt as ts_qice_gt
......@@ -74,36 +74,36 @@ def run(name, config):
fn.mkdir_p(img_dir)
element_list = []
cesm_aavg = cesm_rad.calculate_area_weighted_averages(config)
model_aavg = model_rad.calculate_area_weighted_averages(config)
racmo_aavg = racmo_rad.calculate_area_weighted_averages(config)
table_data = [(key, {'CESM': cesm_aavg[key], 'RACMO 2.3': racmo_aavg[key]}) for key in cesm_aavg.keys()]
table_data = [(key, {'Model': model_aavg[key], 'RACMO 2.3': racmo_aavg[key]}) for key in model_aavg.keys()]
table_dict = OrderedDict(table_data)
# FIXME this is using an undocumented LIVVkit table type, developed for EVE.
table_el = {'Type': 'V-H Table',
'Title': 'Validation',
'TableTitle': 'Area weighted averages of energy balance variables.',
'Headers': ['CESM', 'RACMO 2.3'],
'Headers': ['Model', 'RACMO 2.3'],
'Data': {'': table_dict}
}
element_list.append(table_el)
img_list = []
img_list.extend(cesm_albedo.make_plot(config, out_path=img_dir))
img_list.extend(cesm_latf.make_plot(config, out_path=img_dir))
img_list.extend(cesm_lwsd.make_plot(config, out_path=img_dir))
img_list.extend(cesm_lwsn.make_plot(config, out_path=img_dir))
img_list.extend(cesm_rnet.make_plot(config, out_path=img_dir))
img_list.extend(cesm_senf.make_plot(config, out_path=img_dir))
img_list.extend(cesm_smb.make_plot(config, out_path=img_dir))
img_list.extend(cesm_swsd.make_plot(config, out_path=img_dir))
img_list.extend(cesm_swsn.make_plot(config, out_path=img_dir))
img_list.extend(albedo.make_plot(config, out_path=img_dir))
img_list.extend(latf.make_plot(config, out_path=img_dir))
img_list.extend(lwsd.make_plot(config, out_path=img_dir))
img_list.extend(lwsn.make_plot(config, out_path=img_dir))
img_list.extend(rnet.make_plot(config, out_path=img_dir))
img_list.extend(senf.make_plot(config, out_path=img_dir))
img_list.extend(smb.make_plot(config, out_path=img_dir))
img_list.extend(swsd.make_plot(config, out_path=img_dir))
img_list.extend(swsn.make_plot(config, out_path=img_dir))
img_list.extend(t2m_ann.make_plot(config, out_path=img_dir))
img_list.extend(t2m_djf.make_plot(config, out_path=img_dir))
img_list.extend(t2m_jja.make_plot(config, out_path=img_dir))
img_list.extend(ts_trefht.make_plot(config, out_path=img_dir))
img_list.extend(ts_tsa.make_plot(config, out_path=img_dir))
img_list.extend(ts_qice_gt.make_plot(config, out_path=img_dir))
element_list.append(el.gallery('Figures', img_list))
......
......@@ -12,9 +12,9 @@ from livvkit.util import elements as el
describe = """
Average of the annual average surface albedo over Greenland
for every summer (June-July-August; JJA) from 1980--1999
for E2SM (left) and RACMO 2.3 (middle), and the difference between them (right; E2SM -
for Model (left) and RACMO 2.3 (middle), and the difference between them (right; Model -
RACMO 2.3). The black solid lines denote the 0, 1000, 2000, and 3000 meter
Greenland ice sheet elevation contours as seen by the models (E2SM's contours
Greenland ice sheet elevation contours (Model's contours
shown in the difference plot).
"""
......@@ -22,31 +22,20 @@ title = "Summer surface albedo"
def make_plot(config, out_path='.'):
# ---------------- Data source in TITAN ------------------------
f_cism = os.path.join(config['cism_data'], 'Greenland_5km_v1.1_SacksRev_c110629.nc')
f_e3sm_lnd_climo_jja = os.path.join(config['model_lnd_climos'],
'IGCLM45_MLI_titan_JJA_climo.nc')
f_racmo_alb_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.alb.1980-1999.JJA.nc')
f_racmo_alb_jja_remapped = os.path.join(config['racmo_data'],
'racmo23_GRN_monthly.alb.1980-1999.remap2e3sm.JJA.nc')
f_racmo_mask = os.path.join(config['racmo_data'], 'RACMO23_masks_ZGRN11.nc')
# --------------------------------------------------------------
img_list = []
# f_cism get the following, matrixsize[301,561]
# usrf(0,:,:), lat(0,:,:), lon(0,:,:)
# Goal: plot 3 figures, radiation overlaying on elevation usrf
# read f_cism, the elevation data for E3SM
# adding the [:] transform the variable to type numpy.ndarray
ncid0 = Dataset(f_cism)
ncid0 = Dataset(config['glc_surf'])
usrf = ncid0.variables['usrf'][0, :, :]
lat = ncid0.variables['lat'][0, :, :]
lon = ncid0.variables['lon'][0, :, :]
# read f_e3sm_lnd_climo_jja, variables of E3SM, lat/lon in file1 are 1D vector
ncid1 = Dataset(f_e3sm_lnd_climo_jja)
ncid1 = Dataset(config['lnd_climo'])
fsds = ncid1.variables['FSDS'][0, :, :]
fsr = ncid1.variables['FSR'][0, :, :]
lat1 = ncid1.variables['lat'][:]
......@@ -92,7 +81,7 @@ def make_plot(config, out_path='.'):
wkres = Ngl.Resources()
# wkres.wkOrientation = "portrait" # "portrait" or "landscape"
wks_type = "png"
wks_img = str(os.path.join(out_path, "E2SM_RACMO23_albedo_JJA"))
wks_img = str(os.path.join(out_path, "Model_RACMO23_albedo_JJA"))
wks = Ngl.open_wks(wks_type, wks_img, wkres)
# --- for the map -------
......@@ -115,16 +104,16 @@ def make_plot(config, out_path='.'):
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the albedo contour of E2SM -------
# --- for the albedo contour of the Model -------
res1 = Ngl.Resources()
res1.cnFillPalette = "WhiteYellowOrangeRed"
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.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
res1.lbLabelBarOn = True # Turn on labelbar.
......@@ -138,11 +127,11 @@ def make_plot(config, out_path='.'):
res2.cnFillPalette = "WhiteYellowOrangeRed"
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.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
res2.lbLabelBarOn = True # Turn on labelbar.
......@@ -152,16 +141,16 @@ def make_plot(config, out_path='.'):
res2.sfXArray = lon4
res2.sfYArray = lat4
# --- for the albedo contour of E2SM-RACMO -------
# --- for the albedo contour of Model-RACMO -------
res3 = Ngl.Resources()
res3.cnFillPalette = "BlueWhiteOrangeRed"
res3.nglDraw = False # Don't draw individual plots
res3.nglFrame = False # Don't advance frame.
res3.cnLineLabelsOn = False
res3.cnFillOn = True
res3.cnLinesOn = False
res3.cnLineLabelsOn = False
res3.cnFillMode = "RasterFill"
res3.trGridType = "TriangularMesh"
res3.cnLevelSelectionMode = "ExplicitLevels"
res3.cnLevels = [-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
res3.lbLabelBarOn = True # Turn on labelbar.
......@@ -178,6 +167,7 @@ def make_plot(config, out_path='.'):
sres.cnFillOn = False
sres.cnLinesOn = True
sres.cnLineLabelsOn = False
sres.trGridType = "TriangularMesh"
sres.cnLevelSelectionMode = "ExplicitLevels"
sres.cnLevels = [0, 1000, 2000, 3000]
sres.sfXArray = lon
......@@ -190,6 +180,7 @@ def make_plot(config, out_path='.'):
sres1.cnFillOn = False
sres1.cnLinesOn = True
sres1.cnLineLabelsOn = False
sres1.trGridType = "TriangularMesh"
sres1.cnLevelSelectionMode = "ExplicitLevels"
sres1.cnLevels = [0, 1000, 2000, 3000]
sres1.sfXArray = lon4
......@@ -197,7 +188,6 @@ def make_plot(config, out_path='.'):
# ---- Overlay plots, each one has its own ID
# overlay ice on base map, and then overlay elevation on ice
# usrf is for 5km E2SM; racmo_elev is for RACMO
usrf_plot1 = Ngl.contour(wks, usrf, sres)
usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1)
usrf_plot3 = Ngl.contour(wks, usrf, sres)
......@@ -209,8 +199,7 @@ def make_plot(config, out_path='.'):
# Creat multiple figures and draw, which now contains the elevation and radiation
# "[1,3]" indicates 1 row, 3 columns.
map_title = ["E2SM, Albedo", "RACMO, Albedo", "E2SM-RACMO, Albedo"]
# map_title = ["E2SM, Albedo","RACMO, Albedo"]
map_title = ["Model, Albedo", "RACMO, Albedo", "Model-RACMO, Albedo"]
nmap = 3
plot = []
......
# coding=utf-8
import os
import Ngl
import numpy.ma as ma
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """
Average of the annual average surface albedo over Greenland
for every summer (June-July-August; JJA) from 1980--1999
for CESM (left) and RACMO 2.3 (middle), and the difference between them (right; CESM -
RACMO 2.3). The black solid lines denote the 0, 1000, 2000, and 3000 meter
Greenland ice sheet elevation contours as seen by the models (CESM's contours
shown in the difference plot).
"""
title = "Summer surface albedo"
def make_plot(config, out_path='.'):
# ---------------- Data source in TITAN ------------------------
f_cism = os.path.join(config['cism_data'], 'Greenland_5km_v1.1_SacksRev_c110629.nc')
f_cesm_lnd_climo_jja = os.path.join(config['model_lnd_climos'],
'b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc')
f_racmo_alb_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.alb.1980-1999.JJA.nc')
f_racmo_alb_jja_remapped = os.path.join(config['racmo_data'],
'racmo23_GRN_monthly.alb.1980-1999.remap2cesm.JJA.nc')
f_racmo_mask = os.path.join(config['racmo_data'], 'RACMO23_masks_ZGRN11.nc')
# --------------------------------------------------------------
img_list = []
# f_cism get the following, matrixsize[301,561]
# usrf(0,:,:), lat(0,:,:), lon(0,:,:)
# Goal: plot 3 figures, radiation overlaying on elevation usrf
# read f_cism, the elevation data for CESM
# adding the [:] transform the variable to type numpy.ndarray
ncid0 = Dataset(f_cism)
usrf = ncid0.variables['usrf'][0, :, :]
lat = ncid0.variables['lat'][0, :, :]
lon = ncid0.variables['lon'][0, :, :]
# read f_cesm_lnd_climo_jja, variables of CESM, lat/lon in file1 are 1D vector
ncid1 = Dataset(f_cesm_lnd_climo_jja)
fsds = ncid1.variables['FSDS'][0, :, :]
fsr = ncid1.variables['FSR'][0, :, :]
lat1 = ncid1.variables['lat'][:]
lon1 = ncid1.variables['lon'][:]
albedo = fsr / fsds
# use gris as a mask to mask albedo array
gris_mask = ncid1.variables['gris_mask'][0, :, :]
gris_mask = ma.masked_equal(gris_mask, 0)
albedo_mask = ma.masked_array(albedo, mask=gris_mask.mask)
albedo = albedo_mask
# read f_racmo_alb_jja, variables of RACMO
ncid2 = Dataset(f_racmo_alb_jja)
alb = ncid2.variables['alb'][0, :, :]
# read f_racmo_mask, the lat/lon for RACMO data
ncid4 = Dataset(f_racmo_mask)
lat4 = ncid4.variables['lat'][:]
lon4 = ncid4.variables['lon'][:]
racmo_elev = ncid4.variables['Topography'][0, 0, :, :]
# Use gris as a mask to mask albedo array.
gris_mask = ncid4.variables['GrIS_mask'][:]
gris_mask = ma.masked_equal(gris_mask, 0)
alb_mask = ma.masked_array(alb, mask=gris_mask.mask)
alb = alb_mask
# read input_file3, the remapped RACMO file to calculate difference
ncid3 = Dataset(f_racmo_alb_jja_remapped)
remap_alb = ncid3.variables['alb'][0, :, :]
diff = albedo - remap_alb
# print(np.max(diff))
# print(np.min(diff))
# ------- PLOT --------
# Open a workstation for drawing the plots
wkres = Ngl.Resources()
# wkres.wkOrientation = "portrait" # "portrait" or "landscape"
wks_type = "png"
wks_img = str(os.path.join(out_path, "CESM_RACMO23_albedo_JJA"))
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 = 70. # map area
mres.mpCenterLonF = -44.
mres.mpMinLatF = 57.
mres.mpMaxLatF = 85.
mres.mpMinLonF = -55.
mres.mpMaxLonF = -30.
mres.mpOutlineOn = False
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the albedo contour of CESM -------
res1 = Ngl.Resources()
res1.cnFillPalette = "WhiteYellowOrangeRed"
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.cnLevelSelectionMode = "ExplicitLevels"
res1.cnLevels = [0.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
res1.lbLabelBarOn = True # Turn on labelbar.
res1.lbLabelFontHeightF = 0.04
res1.sfXArray = lon1
res1.sfYArray = lat1
# --- for the albedo contour of RACMO -------
res2 = Ngl.Resources()
res2.cnFillPalette = "WhiteYellowOrangeRed"
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.cnLevelSelectionMode = "ExplicitLevels"
res2.cnLevels = [0.05, 0.25, 0.5, 0.6, 0.7, 0.8, 0.825, 0.85, 0.9]
res2.lbLabelBarOn = True # Turn on labelbar.
res2.lbOrientation = "Vertical" # Verticle labelbar
res2.lbLabelFontHeightF = 0.04 # Make fonts smaller.
res2.sfXArray = lon4
res2.sfYArray = lat4
# --- for the albedo contour of CESM-RACMO -------
res3 = Ngl.Resources()
res3.cnFillPalette = "BlueWhiteOrangeRed"
res3.nglDraw = False # Don't draw individual plots
res3.nglFrame = False # Don't advance frame.
res3.cnLineLabelsOn = False
res3.cnFillOn = True
res3.cnLinesOn = False
res3.cnLineLabelsOn = False
res3.cnFillMode = "RasterFill"
res3.cnLevelSelectionMode = "ExplicitLevels"
res3.cnLevels = [-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
res3.lbLabelBarOn = True # Turn on labelbar.
res3.lbOrientation = "Vertical" # Verticle labelbar
res3.lbLabelFontHeightF = 0.04 # Make fonts smaller.
res3.sfXArray = lon1
res3.sfYArray = lat1
# ---- for the elevation -------
sres = Ngl.Resources()
sres.nglDraw = False # Don't draw individual plots
sres.nglFrame = False # Don't advance frame.
sres.cnFillOn = False
sres.cnLinesOn = True
sres.cnLineLabelsOn = False
sres.cnLevelSelectionMode = "ExplicitLevels"
sres.cnLevels = [0, 1000, 2000, 3000]
sres.sfXArray = lon
sres.sfYArray = lat
# ---- for the elevation of RACMO -------
sres1 = Ngl.Resources()
sres1.nglDraw = False # Don't draw individual plots
sres1.nglFrame = False # Don't advance frame.
sres1.cnFillOn = False
sres1.cnLinesOn = True
sres1.cnLineLabelsOn = False
sres1.cnLevelSelectionMode = "ExplicitLevels"
sres1.cnLevels = [0, 1000, 2000, 3000]
sres1.sfXArray = lon4
sres1.sfYArray = lat4
# ---- Overlay plots, each one has its own ID
# overlay ice on base map, and then overlay elevation on ice
# usrf is for 5km CESM; racmo_elev is for RACMO
usrf_plot1 = Ngl.contour(wks, usrf, sres)
usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1)
usrf_plot3 = Ngl.contour(wks, usrf, sres)
albedo_plot = Ngl.contour(wks, albedo, res1)
alb_plot = Ngl.contour(wks, alb, res2)
# diff_plot = Ngl.contour(wks,Remap_alb,res1) # Remapped RACMO
diff_plot = Ngl.contour(wks, diff, res3)
# Creat multiple figures and draw, which now contains the elevation and radiation
# "[1,3]" indicates 1 row, 3 columns.
map_title = ["CESM, Albedo", "RACMO, Albedo", "CESM-RACMO, Albedo"]
# map_title = ["CESM, Albedo","RACMO, Albedo"]
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], albedo_plot)
Ngl.overlay(plot[0], usrf_plot1)
Ngl.overlay(plot[1], alb_plot)
Ngl.overlay(plot[1], usrf_plot2)
Ngl.overlay(plot[2], diff_plot)
Ngl.overlay(plot[2], usrf_plot3)
Ngl.panel(wks, plot, [1, nmap])
img_link = os.path.join(os.path.basename(out_path),
os.path.basename(wks_img + '.' + wks_type))
img_elem = el.image(title,
' '.join(describe.split()),
img_link)
img_elem['Height'] = config['image_height']
img_list.append(img_elem)
return img_list
# coding=utf-8
import os
import Ngl
import numpy.ma as ma
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """
Average of the annual average surface albedo over Greenland
for every summer (June-July-August; JJA) from 1980--1999
for E2SM (left) and RACMO 2.3 (middle), and the difference between them (right; E2SM -
RACMO 2.3). The black solid lines denote the 0, 1000, 2000, and 3000 meter
Greenland ice sheet elevation contours as seen by the models (E2SM's contours
shown in the difference plot).
"""
title = "Summer surface albedo"
def make_plot(config, out_path='.'):
# ---------------- Data source in TITAN ------------------------
f_cism = os.path.join(config['cism_data'], 'Greenland_5km_v1.1_SacksRev_c110629.nc')
f_e3sm_lnd_climo_jja = os.path.join(config['model_lnd_climos'],
'IGCLM45_MLI_titan_JJA_climo.nc')
f_racmo_alb_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.alb.1980-1999.JJA.nc')
f_racmo_alb_jja_remapped = os.path.join(config['racmo_data'],
'racmo23_GRN_monthly.alb.1980-1999.remap2e3sm.JJA.nc')
f_racmo_mask = os.path.join(config['racmo_data'], 'RACMO23_masks_ZGRN11.nc')
# --------------------------------------------------------------
img_list = []
# f_cism get the following, matrixsize[301,561]
# usrf(0,:,:), lat(0,:,:), lon(0,:,:)
# Goal: plot 3 figures, radiation overlaying on elevation usrf
# read f_cism, the elevation data for E3SM
# adding the [:] transform the variable to type numpy.ndarray
ncid0 = Dataset(f_cism)
usrf = ncid0.variables['usrf'][0, :, :]
lat = ncid0.variables['lat'][0, :, :]
lon = ncid0.variables['lon'][0, :, :]
# read f_e3sm_lnd_climo_jja, variables of E3SM, lat/lon in file1 are 1D vector
ncid1 = Dataset(f_e3sm_lnd_climo_jja)
fsds = ncid1.variables['FSDS'][0, :, :]
fsr = ncid1.variables['FSR'][0, :, :]
lat1 = ncid1.variables['lat'][:]
lon1 = ncid1.variables['lon'][:]
albedo = fsr / fsds
# use gris as a mask to mask albedo array
gris_mask = ncid1.variables['gris_mask'][0, :, :]
gris_mask = ma.masked_equal(gris_mask, 0)
albedo_mask = ma.masked_array(albedo, mask=gris_mask.mask)
albedo = albedo_mask
# read f_racmo_alb_jja, variables of RACMO
ncid2 = Dataset(f_racmo_alb_jja)
alb = ncid2.variables['alb'][0, :, :]
# read f_racmo_mask, the lat/lon for RACMO data
ncid4 = Dataset(f_racmo_mask)
lat4 = ncid4.variables['lat'][:]
lon4 = ncid4.variables['lon'][:]
racmo_elev = ncid4.variables['Topography'][0, 0, :, :]
# Use gris as a mask to mask albedo array.
gris_mask = ncid4.variables['GrIS_mask'][:]
gris_mask = ma.masked_equal(gris_mask, 0)
alb_mask = ma.masked_array(alb, mask=gris_mask.mask)
alb = alb_mask
# read input_file3, the remapped RACMO file to calculate difference
ncid3 = Dataset(f_racmo_alb_jja_remapped)
remap_alb = ncid3.variables['alb'][0, :, :]
diff = albedo - remap_alb
# print(np.max(diff))
# print(np.min(diff))
# ------- PLOT --------
# Open a workstation for drawing the plots
wkres = Ngl.Resources()
# wkres.wkOrientation = "portrait" # "portrait" or "landscape"
wks_type = "png"
wks_img = str(os.path.join(out_path, "E2SM_RACMO23_albedo_JJA"))
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 = 70. # map area
mres.mpCenterLonF = -44.
mres.mpMinLatF = 57.
mres.mpMaxLatF = 85.
mres.mpMinLonF = -55.
mres.mpMaxLonF = -30.
mres.mpOutlineOn = False
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the albedo contour of E2SM -------
res1 = Ngl.Resources()
res1.cnFillPalette = "WhiteYellowOrangeRed"
res1.nglDraw = False # Don't draw individual plots
res1.nglFrame = False # Don't advance frame.