cisma_racmo23_scatter_standalone.py 3.38 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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 -----------------------------------
25
    input_file1 = os.path.join(racmo_path, 'racmo23.ANN_with_latlon_smb.nc')
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    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()