Using masks and computing weighted average

This example is based from xarray example http://xarray.pydata.org/en/stable/examples/area_weighted_temperature.html

Import python packages

import xarray as xr
xr.set_options(display_style='html')
import intake
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
%matplotlib inline

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
col = intake.open_esm_datastore(cat_url)
col

pangeo-cmip6 catalog with 7674 dataset(s) from 514818 asset(s):

	unique
activity_id	18
institution_id	36
source_id	88
experiment_id	170
member_id	657
table_id	37
variable_id	700
grid_label	10
zstore	514818
dcpp_init_year	60
version	736
derived_variable_id	0

Search data

cat = col.search(source_id=['NorESM2-LM'], experiment_id=['historical'], table_id=['Amon'], variable_id=['tas'], member_id=['r1i1p1f1'])
cat.df

	activity_id	institution_id	source_id	experiment_id	member_id	table_id	variable_id	grid_label	zstore	dcpp_init_year	version
0	CMIP	NCC	NorESM2-LM	historical	r1i1p1f1	Amon	tas	gn	gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historica...	NaN	20190815

Create dictionary from the list of datasets we found

• This step may take several minutes so be patient!

dset_dict = cat.to_dataset_dict(zarr_kwargs={'use_cftime':True})

--> The keys in the returned dictionary of datasets are constructed as follows:
	'activity_id.institution_id.source_id.experiment_id.table_id.grid_label'

100.00% [1/1 00:12<00:00]

list(dset_dict.keys())

['CMIP.NCC.NorESM2-LM.historical.Amon.gn']

dset = dset_dict[list(dset_dict.keys())[0]]
dset

<xarray.Dataset> Size: 110MB
Dimensions:         (lat: 96, bnds: 2, lon: 144, member_id: 1,
                     dcpp_init_year: 1, time: 1980)
Coordinates:
    height          float64 8B ...
  * lat             (lat) float64 768B -90.0 -88.11 -86.21 ... 86.21 88.11 90.0
    lat_bnds        (lat, bnds) float64 2kB dask.array<chunksize=(96, 2), meta=np.ndarray>
  * lon             (lon) float64 1kB 0.0 2.5 5.0 7.5 ... 352.5 355.0 357.5
    lon_bnds        (lon, bnds) float64 2kB dask.array<chunksize=(144, 2), meta=np.ndarray>
  * time            (time) object 16kB 1850-01-16 12:00:00 ... 2014-12-16 12:...
    time_bnds       (time, bnds) object 32kB dask.array<chunksize=(1980, 2), meta=np.ndarray>
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) float64 8B nan
Dimensions without coordinates: bnds
Data variables:
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 109MB dask.array<chunksize=(1, 1, 990, 96, 144), meta=np.ndarray>
Attributes: (12/65)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    Hybrid-restart from year 1600-01-01 of ...
    branch_time:                      0.0
    branch_time_in_child:             0.0
    branch_time_in_parent:            430335.0
    ...                               ...
    intake_esm_attrs:variable_id:     tas
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/hi...
    intake_esm_attrs:version:         20190815
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.NCC.NorESM2-LM.historical.Amon.gn

xarray.Dataset

• Dimensions:
  • lat: 96
  • bnds: 2
  • lon: 144
  • member_id: 1
  • dcpp_init_year: 1
  • time: 1980
• Coordinates: (9)
  • height
    ()
    float64
    ...

    axis :

    Z

    long_name :

    height

    positive :

    up

    standard_name :

    height

    units :

    m

    [1 values with dtype=float64]

  • lat
    (lat)
    float64
    -90.0 -88.11 -86.21 ... 88.11 90.0

    axis :

    Y

    bounds :

    lat_bnds

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    array([-90.      , -88.105263, -86.210526, -84.315789, -82.421053, -80.526316,
           -78.631579, -76.736842, -74.842105, -72.947368, -71.052632, -69.157895,
           -67.263158, -65.368421, -63.473684, -61.578947, -59.684211, -57.789474,
           -55.894737, -54.      , -52.105263, -50.210526, -48.315789, -46.421053,
           -44.526316, -42.631579, -40.736842, -38.842105, -36.947368, -35.052632,
           -33.157895, -31.263158, -29.368421, -27.473684, -25.578947, -23.684211,
           -21.789474, -19.894737, -18.      , -16.105263, -14.210526, -12.315789,
           -10.421053,  -8.526316,  -6.631579,  -4.736842,  -2.842105,  -0.947368,
             0.947368,   2.842105,   4.736842,   6.631579,   8.526316,  10.421053,
            12.315789,  14.210526,  16.105263,  18.      ,  19.894737,  21.789474,
            23.684211,  25.578947,  27.473684,  29.368421,  31.263158,  33.157895,
            35.052632,  36.947368,  38.842105,  40.736842,  42.631579,  44.526316,
            46.421053,  48.315789,  50.210526,  52.105263,  54.      ,  55.894737,
            57.789474,  59.684211,  61.578947,  63.473684,  65.368421,  67.263158,
            69.157895,  71.052632,  72.947368,  74.842105,  76.736842,  78.631579,
            80.526316,  82.421053,  84.315789,  86.210526,  88.105263,  90.      ])

  • lat_bnds
    (lat, bnds)
    float64
    dask.array<chunksize=(96, 2), meta=np.ndarray>

    Array Chunk Bytes 1.50 kiB 1.50 kiB Shape (96, 2) (96, 2) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • lon
    (lon)
    float64
    0.0 2.5 5.0 ... 352.5 355.0 357.5

    axis :

    X

    bounds :

    lon_bnds

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
            25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
            50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
            75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
           100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
           125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
           150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
           175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
           200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
           325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
           350. , 352.5, 355. , 357.5])

  • lon_bnds
    (lon, bnds)
    float64
    dask.array<chunksize=(144, 2), meta=np.ndarray>

    Array Chunk Bytes 2.25 kiB 2.25 kiB Shape (144, 2) (144, 2) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • time
    (time)
    object
    1850-01-16 12:00:00 ... 2014-12-...

    axis :

    T

    bounds :

    time_bnds

    long_name :

    time

    standard_name :

    time

    array([cftime.DatetimeNoLeap(1850, 1, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 3, 16, 12, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeNoLeap(2014, 10, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2014, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2014, 12, 16, 12, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • time_bnds
    (time, bnds)
    object
    dask.array<chunksize=(1980, 2), meta=np.ndarray>

    Array Chunk Bytes 30.94 kiB 30.94 kiB Shape (1980, 2) (1980, 2) Dask graph 1 chunks in 2 graph layers Data type object numpy.ndarray

  • member_id
    (member_id)
    object
    'r1i1p1f1'

    array(['r1i1p1f1'], dtype=object)

  • dcpp_init_year
    (dcpp_init_year)
    float64
    nan

    array([nan])

• Data variables: (1)
  • tas
    (member_id, dcpp_init_year, time, lat, lon)
    float32
    dask.array<chunksize=(1, 1, 990, 96, 144), meta=np.ndarray>

    cell_measures :

    area: areacella

    cell_methods :

    area: time: mean

    comment :

    near-surface (usually, 2 meter) air temperature

    history :

    2019-08-15T12:42:20Z altered by CMOR: Treated scalar dimension: 'height'. 2019-08-15T12:42:21Z altered by CMOR: Converted type from 'd' to 'f'.

    long_name :

    Near-Surface Air Temperature

    original_name :

    TREFHT

    standard_name :

    air_temperature

    units :

    K

    Array Chunk Bytes 104.41 MiB 52.21 MiB Shape (1, 1, 1980, 96, 144) (1, 1, 990, 96, 144) Dask graph 2 chunks in 3 graph layers Data type float32 numpy.ndarray

• Indexes: (5)
  • lat
    PandasIndex

    PandasIndex(Index([              -90.0,  -88.10526315789474,  -86.21052631578948,
             -84.3157894736842,  -82.42105263157895,  -80.52631578947368,
            -78.63157894736842,  -76.73684210526316,  -74.84210526315789,
            -72.94736842105263,  -71.05263157894737,  -69.15789473684211,
            -67.26315789473685,  -65.36842105263158, -63.473684210526315,
           -61.578947368421055, -59.684210526315795,  -57.78947368421053,
            -55.89473684210527,               -54.0,  -52.10526315789474,
            -50.21052631578947,  -48.31578947368421,  -46.42105263157895,
           -44.526315789473685, -42.631578947368425,  -40.73684210526316,
             -38.8421052631579,  -36.94736842105264,  -35.05263157894737,
            -33.15789473684211, -31.263157894736842, -29.368421052631582,
           -27.473684210526322, -25.578947368421055, -23.684210526315795,
           -21.789473684210535,  -19.89473684210526,               -18.0,
            -16.10526315789474,  -14.21052631578948,  -12.31578947368422,
           -10.421052631578945,  -8.526315789473685,  -6.631578947368425,
            -4.736842105263165, -2.8421052631579045, -0.9473684210526301,
            0.9473684210526301,  2.8421052631578902,    4.73684210526315,
             6.631578947368411,   8.526315789473685,  10.421052631578945,
            12.315789473684205,  14.210526315789465,  16.105263157894726,
                          18.0,   19.89473684210526,   21.78947368421052,
             23.68421052631578,   25.57894736842104,  27.473684210526315,
            29.368421052631575,  31.263157894736835,  33.157894736842096,
            35.052631578947356,   36.94736842105263,   38.84210526315789,
             40.73684210526315,   42.63157894736841,   44.52631578947367,
             46.42105263157893,   48.31578947368419,   50.21052631578948,
             52.10526315789474,                54.0,   55.89473684210526,
             57.78947368421052,   59.68421052631578,   61.57894736842104,
              63.4736842105263,   65.36842105263156,   67.26315789473682,
             69.15789473684211,   71.05263157894737,   72.94736842105263,
             74.84210526315789,   76.73684210526315,   78.63157894736841,
             80.52631578947367,   82.42105263157893,   84.31578947368419,
             86.21052631578945,   88.10526315789474,                90.0],
          dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,  22.5,
           ...
           335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0, 357.5],
          dtype='float64', name='lon', length=144))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1850-01-16 12:00:00, 1850-02-15 00:00:00, 1850-03-16 12:00:00,
                 1850-04-16 00:00:00, 1850-05-16 12:00:00, 1850-06-16 00:00:00,
                 1850-07-16 12:00:00, 1850-08-16 12:00:00, 1850-09-16 00:00:00,
                 1850-10-16 12:00:00,
                 ...
                 2014-03-16 12:00:00, 2014-04-16 00:00:00, 2014-05-16 12:00:00,
                 2014-06-16 00:00:00, 2014-07-16 12:00:00, 2014-08-16 12:00:00,
                 2014-09-16 00:00:00, 2014-10-16 12:00:00, 2014-11-16 00:00:00,
                 2014-12-16 12:00:00],
                dtype='object', length=1980, calendar='noleap', freq=None))

  • member_id
    PandasIndex

    PandasIndex(Index(['r1i1p1f1'], dtype='object', name='member_id'))

  • dcpp_init_year
    PandasIndex

    PandasIndex(Index([nan], dtype='float64', name='dcpp_init_year'))

• Attributes: (65)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  CMIP

  branch_method :

  Hybrid-restart from year 1600-01-01 of piControl

  branch_time :

  0.0

  branch_time_in_child :

  0.0

  branch_time_in_parent :

  430335.0

  cmor_version :

  3.5.0

  contact :

  Please send any requests or bug reports to noresm-ncc@met.no.

  creation_date :

  2019-08-15T12:42:21Z

  data_specs_version :

  01.00.31

  experiment :

  all-forcing simulation of the recent past

  experiment_id :

  historical

  external_variables :

  areacella

  forcing_index :

  1

  frequency :

  mon

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.NCC.NorESM2-LM.historical.none.r1i1p1f1

  grid :

  finite-volume grid with 1.9x2.5 degree lat/lon resolution

  grid_label :

  gn

  history :

  2019-08-15T12:42:21Z ; CMOR rewrote data to be consistent with CMIP6, CF-1.7 CMIP-6.2 and CF standards.

  initialization_index :

  1

  institution :

  NorESM Climate modeling Consortium consisting of CICERO (Center for International Climate and Environmental Research, Oslo 0349), MET-Norway (Norwegian Meteorological Institute, Oslo 0313), NERSC (Nansen Environmental and Remote Sensing Center, Bergen 5006), NILU (Norwegian Institute for Air Research, Kjeller 2027), UiB (University of Bergen, Bergen 5007), UiO (University of Oslo, Oslo 0313) and UNI (Uni Research, Bergen 5008), Norway. Mailing address: NCC, c/o MET-Norway, Henrik Mohns plass 1, Oslo 0313, Norway

  institution_id :

  NCC

  license :

  CMIP6 model data produced by NCC is licensed under a Creative Commons Attribution ShareAlike 4.0 International License (https://creativecommons.org/licenses). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file) and at https:///pcmdi.llnl.gov/. The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  CMIP6

  model_id :

  NorESM2-LM

  nominal_resolution :

  250 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  piControl

  parent_mip_era :

  CMIP6

  parent_source_id :

  NorESM2-LM

  parent_sub_experiment_id :

  none

  parent_time_units :

  days since 0421-01-01

  parent_variant_label :

  r1i1p1f1

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  atmos

  run_variant :

  N/A

  source :

  NorESM2-LM (2017): aerosol: OsloAero atmos: CAM-OSLO (2 degree resolution; 144 x 96; 32 levels; top level 3 mb) atmosChem: OsloChemSimp land: CLM landIce: CISM ocean: MICOM (1 degree resolution; 360 x 384; 70 levels; top grid cell minimum 0-2.5 m [native model uses hybrid density and generic upper-layer coordinate interpolated to z-level for contributed data]) ocnBgchem: HAMOCC seaIce: CICE

  source_id :

  NorESM2-LM

  source_type :

  AOGCM

  status :

  2020-04-30;created; by gcs.cmip6.ldeo@gmail.com

  sub_experiment :

  none

  sub_experiment_id :

  none

  table_id :

  Amon

  table_info :

  Creation Date:(24 July 2019) MD5:08e314340b9dd440c36c1b8e484514f4

  title :

  NorESM2-LM output prepared for CMIP6

  tracking_id :

  hdl:21.14100/2486cf87-033c-4848-ab3e-e828c3b7c759 hdl:21.14100/8d2ec216-30d3-4f90-8664-7cf75e29780b hdl:21.14100/9f798305-d0cf-4eac-89b0-b13dc25574ca hdl:21.14100/e0da9277-4784-446b-b3b5-ba33bb18cc3a hdl:21.14100/ca48487c-955a-4b15-9582-5ce55ce39c73 hdl:21.14100/fc7602aa-fd09-49d8-94c6-e619622c211b hdl:21.14100/cd51a9c2-0a94-45b2-8d9b-59bcb3467214 hdl:21.14100/9b7ad334-1217-4229-98e8-b5a9287c0e5e hdl:21.14100/043ca527-40c7-451f-990b-ccc94720f612 hdl:21.14100/a9613498-ea0d-4be4-97e8-77661872798d hdl:21.14100/52699832-d86c-491b-a0a4-2f27d9ee4d69 hdl:21.14100/358c4a18-33b2-4c44-bd30-f6125a4df63c hdl:21.14100/77254eaf-90ca-44e7-9a3c-763a840e188f hdl:21.14100/62aed49c-6ca5-4027-80fe-4ddd99af575a hdl:21.14100/83ffe172-769b-42b1-b979-b92d768f60fd hdl:21.14100/91a8c4e2-d33b-4b87-a978-0247d8112ef3 hdl:21.14100/e8fb384f-d5a7-41e1-bf12-b8de3e2ffb80

  variable_id :

  tas

  variant_label :

  r1i1p1f1

  netcdf_tracking_ids :

  hdl:21.14100/2486cf87-033c-4848-ab3e-e828c3b7c759 hdl:21.14100/8d2ec216-30d3-4f90-8664-7cf75e29780b hdl:21.14100/9f798305-d0cf-4eac-89b0-b13dc25574ca hdl:21.14100/e0da9277-4784-446b-b3b5-ba33bb18cc3a hdl:21.14100/ca48487c-955a-4b15-9582-5ce55ce39c73 hdl:21.14100/fc7602aa-fd09-49d8-94c6-e619622c211b hdl:21.14100/cd51a9c2-0a94-45b2-8d9b-59bcb3467214 hdl:21.14100/9b7ad334-1217-4229-98e8-b5a9287c0e5e hdl:21.14100/043ca527-40c7-451f-990b-ccc94720f612 hdl:21.14100/a9613498-ea0d-4be4-97e8-77661872798d hdl:21.14100/52699832-d86c-491b-a0a4-2f27d9ee4d69 hdl:21.14100/358c4a18-33b2-4c44-bd30-f6125a4df63c hdl:21.14100/77254eaf-90ca-44e7-9a3c-763a840e188f hdl:21.14100/62aed49c-6ca5-4027-80fe-4ddd99af575a hdl:21.14100/83ffe172-769b-42b1-b979-b92d768f60fd hdl:21.14100/91a8c4e2-d33b-4b87-a978-0247d8112ef3 hdl:21.14100/e8fb384f-d5a7-41e1-bf12-b8de3e2ffb80

  version_id :

  v20190815

  intake_esm_vars :

  ['tas']

  intake_esm_attrs:activity_id :

  CMIP

  intake_esm_attrs:institution_id :

  NCC

  intake_esm_attrs:source_id :

  NorESM2-LM

  intake_esm_attrs:experiment_id :

  historical

  intake_esm_attrs:member_id :

  r1i1p1f1

  intake_esm_attrs:table_id :

  Amon

  intake_esm_attrs:variable_id :

  tas

  intake_esm_attrs:grid_label :

  gn

  intake_esm_attrs:zstore :

  gs://cmip6/CMIP6/CMIP/NCC/NorESM2-LM/historical/r1i1p1f1/Amon/tas/gn/v20190815/

  intake_esm_attrs:version :

  20190815

  intake_esm_attrs:_data_format_ :

  zarr

  intake_esm_dataset_key :

  CMIP.NCC.NorESM2-LM.historical.Amon.gn

Plot the first timestep

projection = ccrs.Mercator(central_longitude=-10)

f, ax = plt.subplots(subplot_kw=dict(projection=projection))

dset['tas'].isel(time=0).plot(transform=ccrs.PlateCarree(), cbar_kwargs=dict(shrink=0.7), cmap='coolwarm')
ax.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f1e479d0c90>

Compute weighted mean

1. Creating weights: for a rectangular grid the cosine of the latitude is proportional to the grid cell area.

2. Compute weighted mean values

def computeWeightedMean(ds):
    # Compute weights based on the xarray you pass
    weights = np.cos(np.deg2rad(ds.lat))
    weights.name = "weights"
    # Compute weighted mean
    air_weighted = ds.weighted(weights)
    weighted_mean = air_weighted.mean(("lon", "lat"))
    return weighted_mean

Compute weighted average over the entire globe

weighted_mean = computeWeightedMean(dset)

Comparison with unweighted mean

• We select a time range

• Note how the weighted mean temperature is higher than the unweighted.

weighted_mean['tas'].sel(time=slice('2000-01-01', '2010-01-01')).plot(label="weighted")
dset['tas'].sel(time=slice('2000-01-01', '2010-01-01')).mean(("lon", "lat")).plot(label="unweighted")

plt.legend()

<matplotlib.legend.Legend at 0x7f1e47a58410>

Compute Weigted arctic average

Let’s try to also take only the data above 60$^\circ$

weighted_mean = computeWeightedMean(dset.where(dset['lat']>60.))

weighted_mean['tas'].sel(time=slice('2000-01-01', '2010-01-01')).plot(label="weighted")
dset['tas'].where(dset['lat']>60.).sel(time=slice('2000-01-01', '2010-01-01')).mean(("lon", "lat")).plot(label="unweighted")

plt.legend()

<matplotlib.legend.Legend at 0x7f1e46139ed0>