Skip to content

GitLab

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

Fixes to dynamics extension -- now works on Rhea (only)

parent ba4e50ae
Hide whitespace changes
Inline Side-by-side
dynamics/dynamics/cisma_js_standalone.py 0 → 100755
View file @ 95a7ddaf
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()
dynamics/dynamics/cisma_pm_standalone.py 0 → 100755
View file @ 95a7ddaf
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()
dynamics/dynamics/cisma_racmo23_scatter_standalone.py 0 → 100755
View file @ 95a7ddaf
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, 'climos/racmo23_GRN_monthly.smb.1980-1999.ANN.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()
# read remapped CISMA data
ncid2 = Dataset(input_file2, mode='r')
cisma = ncid2.variables['acab'][0, :, :]
ncid2.close()
cisma = cisma * 1000 # unit transfer
# locate the greenland area using the GreenLand mask
ncid3 = Dataset(input_file3, mode='r')
gris_mask = ncid3.variables['GrIS_mask'][:]
gris_mask = ma.masked_equal(gris_mask, 0)
ncid3.close()
racmo_gm = ma.masked_array(racmo, mask=gris_mask.mask)
cisma_gm = ma.masked_array(cisma, mask=gris_mask.mask)
# Reshape the 2d matrix to 1d vector, columnwise
racmo_1d = np.reshape(racmo_gm, (np.product(racmo_gm.shape),))
cisma_1d = np.reshape(cisma_gm, (np.product(cisma_gm.shape),))
# print(np.max(racmo_1d))
# print(np.min(racmo_1d))
# print(np.max(cisma_1d))
# print(np.min(cisma_1d))
# ------- PLOT JJA --------
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111)
ax.plot([-4000, 5000], [-4000, 5000], 'r-')
ax.scatter(x=racmo_1d, y=cisma_1d, color='b', marker='o')
ax.set_aspect(1. / ax.get_data_ratio()) # make axes square
ax.set_xlim(-4000, 5000)
ax.set_ylim(-4000, 5000)
ax.xaxis.set_ticks(np.arange(-4000, 6000, 1000))
ax.yaxis.set_ticks(np.arange(-4000, 6000, 1000))
ax.set_xlabel('RACMO', fontsize=14)
ax.set_ylabel('CISM-A', fontsize=14)
ax.text(-3800, 3000, 'Surface mass balance', fontsize=18)
img_path = os.path.join(out_path, 'scatter_CISMA_RACMO23.png')
plt.savefig(img_path)
plt.close()
img_link = os.path.join(os.path.basename(out_path),
os.path.basename(img_path))
img_elem = el.image('scatter_CISMA_RACMO23',
' '.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()
dynamics/dynamics/cisma_racmo23_standalone.py 0 → 100755
View file @ 95a7ddaf
import os
import Ngl
import numpy.ma as ma
from netCDF4 import Dataset
from livvkit.util import elements as el
describe = """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/'):
# --------------------------------------------------------------
# first use remap_CISMA_to_RACMO.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, 'climos/racmo23_GRN_monthly.smb.1980-1999.ANN.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_file4 = os.path.join(cism_path, 'CISM-Albany/postproc/albany.out.tstep.ANN_acab.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 racmo source
ncid1 = Dataset(input_file1, 'r')
racmo = ncid1.variables['smb'][0, :, :]
lat1 = ncid1.variables['lat'][:]
lon1 = ncid1.variables['lon'][:]
# read remapped CISMA data
ncid2 = Dataset(input_file2, 'r')
cisma = ncid2.variables['acab'][0, :, :]
lat2 = ncid2.variables['y1'][:]
lon2 = ncid2.variables['x1'][:]
cisma = cisma * 1000 # unit transfer
# read original CISMA data
ncid4 = Dataset(input_file4, 'r')
cisma_ori = ncid4.variables['acab'][0, :, :]
cisma_ori = cisma_ori * 1000 # unit transfer
# read lat and lon for original CISMA
ncid6 = Dataset(input_file6, 'r')
lat3 = ncid6.variables['lat'][:]
lon3 = ncid6.variables['lon'][:]
# read simulated thk of CISM
ncid5 = Dataset(input_file5, 'r')
thk = ncid5.variables['thk'][0, :, :]
# locate the greenland area using the GreenLand mask
ncid3 = Dataset(input_file3, 'r')
gris_mask = ncid3.variables['GrIS_mask'][:]
gris_mask = ma.masked_equal(gris_mask, 0)
racmo = ma.masked_array(racmo, mask=gris_mask.mask)
cisma = ma.masked_array(cisma, mask=gris_mask.mask)
thk_mask = ma.masked_less(thk, 0.01)
cisma_ori = ma.masked_array(cisma_ori, mask=thk_mask.mask)
# diff = cisma - racmo
# print(np.max(diff))
# print(np.min(diff))
# ------- 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_RACMO23"))
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 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 = np.arange(-4000,6000,1000)
res1.cnLevels = (-4000, -2000, -1000, -700, -500, -300, -200, -100, -50, -20,
20, 50, 100, 200, 300, 400, 700, 1000, 2000, 5000)
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