Commit ec1f19f5 authored by Evans, Katherine's avatar Evans, Katherine
Browse files

change scripts to point to E3SM data, update E3SM data. CESM now broken....

change scripts to point to E3SM data, update E3SM data. CESM now broken. timeseries producing nans and still needs a fix. Still needs updated E3SM elevation
parent c454342a
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.
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.
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.
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.
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.
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.
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.
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.
......@@ -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 CESM (left) and RACMO 2.3 (middle), and the difference between them (right; CESM -
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 (CESM's contours
Greenland ice sheet elevation contours as seen by the models (E2SM's contours
shown in the difference plot).
"""
......@@ -24,11 +24,11 @@ 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_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.remap2cesm.JJA.nc')
'racmo23_GRN_monthly.alb.1980-1999.remap2e3sm.JJA.nc')
f_racmo_mask = os.path.join(config['racmo_data'], 'RACMO23_masks_ZGRN11.nc')
# --------------------------------------------------------------
......@@ -38,15 +38,15 @@ def make_plot(config, out_path='.'):
# usrf(0,:,:), lat(0,:,:), lon(0,:,:)
# Goal: plot 3 figures, radiation overlaying on elevation usrf
# read f_cism, the elevation data for CESM
# 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_cesm_lnd_climo_jja, variables of CESM, lat/lon in file1 are 1D vector
ncid1 = Dataset(f_cesm_lnd_climo_jja)
# 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'][:]
......@@ -92,7 +92,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, "CESM_RACMO23_albedo_JJA"))
wks_img = str(os.path.join(out_path, "E2SM_RACMO23_albedo_JJA"))
wks = Ngl.open_wks(wks_type, wks_img, wkres)
# --- for the map -------
......@@ -115,7 +115,7 @@ def make_plot(config, out_path='.'):
mres.mpFillOn = False
mres.mpPerimOn = True # add box around map
# --- for the albedo contour of CESM -------
# --- for the albedo contour of E2SM -------
res1 = Ngl.Resources()
res1.cnFillPalette = "WhiteYellowOrangeRed"
res1.nglDraw = False # Don't draw individual plots
......@@ -152,7 +152,7 @@ def make_plot(config, out_path='.'):
res2.sfXArray = lon4
res2.sfYArray = lat4
# --- for the albedo contour of CESM-RACMO -------
# --- for the albedo contour of E2SM-RACMO -------
res3 = Ngl.Resources()
res3.cnFillPalette = "BlueWhiteOrangeRed"
res3.nglDraw = False # Don't draw individual plots
......@@ -197,7 +197,7 @@ 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 CESM; racmo_elev is for RACMO
# 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 +209,8 @@ 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 = ["CESM, Albedo", "RACMO, Albedo", "CESM-RACMO, Albedo"]
# map_title = ["CESM, Albedo","RACMO, Albedo"]
map_title = ["E2SM, Albedo", "RACMO, Albedo", "E2SM-RACMO, Albedo"]
# map_title = ["E2SM, Albedo","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.
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 E2SM-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 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)
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 = ["E2SM, Albedo", "RACMO, Albedo", "E2SM-RACMO, Albedo"]
# map_title = ["E2SM, 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, <