Kenya_AGBpot/plot_maps.py at master · GCEL/Kenya_AGBpot · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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
"""
27/07/2018 - JFE
this script plots maps of potential biomass in Kenya


"""

import matplotlib.pyplot as plt
import numpy as np
import os
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from netCDF4 import Dataset
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from mpl_toolkits.axes_grid1 import AxesGrid
from cartopy.mpl.geoaxes import GeoAxes
from osgeo import gdal

import sys

path = '/disk/scratch/local.2/jexbraya/kenya_ODA/'

nc_med = Dataset(path+'/output/Kenya_ODA_v31_AGBpot_mean_WC2_SOTWIS_GridSearch.nc')

# load observed and potential AGB
obs = nc_med.variables['AGB_mean'][:]
pot = nc_med.variables['AGBpot_mean'][:]

# JFE added forest mask to only plot forests
frst = nc_med.variables['training'][:] == 2
#replace places outside forests in obs as 0
obs[~frst] = 0.

#set extent and instantiate mappable items
ext = [33.5,42.5,-5,6]

brd= cfeat.BORDERS;brd.scale='10m'
lk = cfeat.LAKES;lk.scale='10m'
co = cfeat.COASTLINE; co.scale = '110m'
oc = cfeat.OCEAN; co.scale = '110m'
la = cfeat.LAND;la.scale='110m'

titles = ['a) AGB','b) Potential Forest Biomass (PFB)','c) PFB-AGB']

#define projection
prj=ccrs.PlateCarree()
axes_class = (GeoAxes,dict(map_projection=prj))

#create map and plot using axesgrid
figmaps = plt.figure('maps',figsize=(12,8));figmaps.clf()

axgr = AxesGrid(figmaps,111,nrows_ncols=(1,3),axes_class=axes_class,label_mode='',cbar_mode='each',cbar_pad = 0.04,cbar_size="5%",axes_pad=1.)

vmn = [0,0,0]
vmx = [200,200,100]
cmaps = ['viridis','viridis','plasma']

obs[(nc_med.variables['training'][:]==1)] = 0.
obs.data[(obs.mask)*(~pot.mask)] = 0.
obs.mask[(obs.mask)*(~pot.mask)] = False

titles = ['a) AGB$_{2015}$','b) AGB$_{pot}$','c) AGB$_{pot}$ - AGB$_{2015}$']

#iterate maps on the grid
for mm,map2plot in enumerate([obs,pot,pot-obs]):

    #plot
    ax = axgr[mm]
    im = ax.imshow(map2plot,origin='upper',vmin = vmn[mm],vmax=vmx[mm],extent=ext,interpolation='nearest',cmap=cmaps[mm])
    #add colorbar and label on rightmost one
    cb = axgr.cbar_axes[mm].colorbar(im)
    if mm ==2:
        cb.set_label_text('Mg ha$^{-1}$')
    #add features: ocean and land to have map on a light grey background
    ax.add_feature(oc,facecolor='silver',zorder=-1)
    ax.add_feature(la,facecolor='silver',zorder = -1)


    #ax.text(0.99,0.99,titles[mm],transform = ax.transAxes,va='top',ha='right',weight='bold',size='medium')
    ax.set_title(titles[mm])

    #set the lat/lon
    ax.set_xticks(np.arange(34,42.1,2))
    ax.xaxis.set_major_formatter(LongitudeFormatter())

    ax.set_yticks(np.arange(-4,5.1,2))
    ax.yaxis.set_major_formatter(LatitudeFormatter())

    print(mm, (map2plot*nc_med.variables['areas'][:]).sum()*1e-13)

figmaps.show()
#figmaps.savefig('figures/compare_maps_V2_WC2_SOTWIS.png', bbox_inches='tight')