In [42]:

import xarray as xr
import numpy as np
import pandas as pd

import seaborn as sns
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
%matplotlib inline

In [43]:

#%config InlineBackend.close_figures=False # tell pyplot not to close figures
#np.warnings.filterwarnings('ignore') # silence numpy warning

In [44]:

# path to netcdf file
f = r'\\...\precip.comb.v2018to2016-v6monitorafter.total.nc'

In [45]:

# open netcdf file
ds = xr.open_dataset(f)
ds.dims

Out[45]:

Frozen({'lat': 180, 'lon': 360, 'time': 1548})

In [46]:

# explore variables
ds.var

Out[46]:

<bound method ImplementsDatasetReduce._reduce_method.<locals>.wrapped_func of <xarray.Dataset>
Dimensions:  (lat: 180, lon: 360, time: 1548)
Coordinates:
  * lat      (lat) float32 89.5 88.5 87.5 86.5 85.5 ... -86.5 -87.5 -88.5 -89.5
  * lon      (lon) float32 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * time     (time) datetime64[ns] 1891-01-01 1891-02-01 ... 2019-12-01
Data variables:
    precip   (time, lat, lon) float32 ...
Attributes:
    Original_Source:  http://www.dwd.de/en/FundE/Klima/KLIS/int/GPCC/GPCC.htm...
    Reference:        Users of the data sets are kindly requested to give fee...
    original_source:  ftp://ftp-anon.dwd.de/pub/data/gpcc/html/download_gate....
    Conventions:      CF 1.0
    dataset_title:    Global Precipitation Climatology Centre (GPCC)
    title:            GPCC Full Data Reanalysis Version 2018 1.0x1.0 Monthly ...
    history:          Created 09/2018  based on V2018 data obtained via ftp
    References:       https://www.psl.noaa.gov/data/gridded/data.gpcc.html>

In [47]:

# ewxplore precipitation variable
prec = ds.precip
prec

Out[47]:

<xarray.DataArray 'precip' (time: 1548, lat: 180, lon: 360)>
[100310400 values with dtype=float32]
Coordinates:
  * lat      (lat) float32 89.5 88.5 87.5 86.5 85.5 ... -86.5 -87.5 -88.5 -89.5
  * lon      (lon) float32 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
  * time     (time) datetime64[ns] 1891-01-01 1891-02-01 ... 2019-12-01
Attributes:
    long_name:     GPCC Monthly total of precipitation
    statistic:     Total
    valid_range:   [   0. 8000.]
    parent_stat:   Observations
    var_desc:      Precipitation
    units:         mm
    level:         Surface
    actual_range:  [   0.   3349.61]
    dataset:       GPCC Precipitation 1.0degree V2018 Full Reanalysis

xarray.DataArray

'precip'

time: 1548
lat: 180
lon: 360

...
```
[100310400 values with dtype=float32]
```

Coordinates: (3)

lat

(lat)

float32

89.5 88.5 87.5 ... -88.5 -89.5

units :: degrees_north
actual_range :: [-89.5 89.5]
long_name :: Latitude
standard_name :: latitude
axis :: Y
coordinate_defines :: point

array([ 89.5,  88.5,  87.5,  86.5,  85.5,  84.5,  83.5,  82.5,  81.5,  80.5,
        79.5,  78.5,  77.5,  76.5,  75.5,  74.5,  73.5,  72.5,  71.5,  70.5,
        69.5,  68.5,  67.5,  66.5,  65.5,  64.5,  63.5,  62.5,  61.5,  60.5,
        59.5,  58.5,  57.5,  56.5,  55.5,  54.5,  53.5,  52.5,  51.5,  50.5,
        49.5,  48.5,  47.5,  46.5,  45.5,  44.5,  43.5,  42.5,  41.5,  40.5,
        39.5,  38.5,  37.5,  36.5,  35.5,  34.5,  33.5,  32.5,  31.5,  30.5,
        29.5,  28.5,  27.5,  26.5,  25.5,  24.5,  23.5,  22.5,  21.5,  20.5,
        19.5,  18.5,  17.5,  16.5,  15.5,  14.5,  13.5,  12.5,  11.5,  10.5,
         9.5,   8.5,   7.5,   6.5,   5.5,   4.5,   3.5,   2.5,   1.5,   0.5,
        -0.5,  -1.5,  -2.5,  -3.5,  -4.5,  -5.5,  -6.5,  -7.5,  -8.5,  -9.5,
       -10.5, -11.5, -12.5, -13.5, -14.5, -15.5, -16.5, -17.5, -18.5, -19.5,
       -20.5, -21.5, -22.5, -23.5, -24.5, -25.5, -26.5, -27.5, -28.5, -29.5,
       -30.5, -31.5, -32.5, -33.5, -34.5, -35.5, -36.5, -37.5, -38.5, -39.5,
       -40.5, -41.5, -42.5, -43.5, -44.5, -45.5, -46.5, -47.5, -48.5, -49.5,
       -50.5, -51.5, -52.5, -53.5, -54.5, -55.5, -56.5, -57.5, -58.5, -59.5,
       -60.5, -61.5, -62.5, -63.5, -64.5, -65.5, -66.5, -67.5, -68.5, -69.5,
       -70.5, -71.5, -72.5, -73.5, -74.5, -75.5, -76.5, -77.5, -78.5, -79.5,
       -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5],
      dtype=float32)

lon
(lon)
float32
0.5 1.5 2.5 ... 357.5 358.5 359.5
long_name :
Longitude
units :
degrees_east
standard_name :
longitude
actual_range :
[ 0.5 359.5]
axis :
X
coordinate_defines :
point
```
array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5], dtype=float32)
```

time

(time)

datetime64[ns]

1891-01-01 ... 2019-12-01

long_name :: Time
delta_t :: 0000-01-00 00:00:00
avg_period :: 0000-01-00 00:00:00
standard_name :: time
axis :: T
coordinate_defines :: start
actual_range :: [33237. 80322.]

array(['1891-01-01T00:00:00.000000000', '1891-02-01T00:00:00.000000000',
       '1891-03-01T00:00:00.000000000', ..., '2019-10-01T00:00:00.000000000',
       '2019-11-01T00:00:00.000000000', '2019-12-01T00:00:00.000000000'],
      dtype='datetime64[ns]')

Attributes: (9)
long_name :
GPCC Monthly total of precipitation
statistic :
Total
valid_range :
[ 0. 8000.]
parent_stat :
Observations
var_desc :
Precipitation
units :
mm
level :
Surface
actual_range :
[ 0. 3349.61]
dataset :
GPCC Precipitation 1.0degree V2018 Full Reanalysis

In [48]:

# extract data for specific date
oct2018 = prec.sel(time='2018-10-01')

In [49]:

# plot data
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines() 
oct2018.plot()
plt.show()

In [50]:

# extract year and month precipitation for a specific location: vanocuver BC
prec_slice = prec.sel(time=slice("2010-01-01", "2019-12-31"))

pre_monthly_loc = prec_slice.sel(lat = 49, 
                           lon = -123+360, # add 360 to convert negative longs from (-180, 180) to (0, 360)
                           method='nearest')

df = pre_monthly_loc.to_dataframe()
df.reset_index(drop=False, inplace=True)
df = df[['time', 'precip']]
df['year'] = pd.DatetimeIndex(df['time']).year
df['Month'] = pd.DatetimeIndex(df['time']).month
df['month'] = pd.DatetimeIndex(df['time']).month_name().str.slice(stop=3)
df

Out[50]:

	time	precip	year	Month	month
0	2010-01-01	312.390015	2010	1	Jan
1	2010-02-01	200.419998	2010	2	Feb
2	2010-03-01	205.570007	2010	3	Mar
3	2010-04-01	177.970001	2010	4	Apr
4	2010-05-01	131.910004	2010	5	May
...	...	...	...	...	...
115	2019-08-01	58.080002	2019	8	Aug
116	2019-09-01	189.190002	2019	9	Sep
117	2019-10-01	231.259995	2019	10	Oct
118	2019-11-01	186.910004	2019	11	Nov
119	2019-12-01	347.899994	2019	12	Dec

120 rows × 5 columns

In [51]:

# Create a Heatmap for year/month precipitation in Vancouver, BC
data = df.pivot(["Month","month"], "year", "precip")

data = data.sort_values(by='Month',ascending=False)

fig = plt.gcf()
fig.set_size_inches(18, 6)

ax = sns.heatmap(data, cmap="YlGnBu").set(title='Monthly Precipitation in Vancouver BC from 2010 - 2019')

In [52]:

# extract oct 2018 precipitation in British Columbia - this time using .where instead of .sel
bc_oct2018_prec = oct2018.where((oct2018.lat > 45) & 
                                (oct2018.lat < 70) & 
                                (oct2018.lon > -137+360) & # add 360 to convert negative longs from (-180, 180) to (0, 360)
                                (oct2018.lon < -114+360),
                                drop=True)

In [53]:

# export extracted data to dataframe
df = bc_oct2018_prec.to_dataframe()
df.tail()

Out[53]:

		time	precip
lat	lon
45.5	241.5	2018-10-01	56.549999
	242.5	2018-10-01	52.779999
	243.5	2018-10-01	51.759998
	244.5	2018-10-01	73.239998
	245.5	2018-10-01	47.490002

In [54]:

# calculate annual precpitation from monthly data
annual_prec = ds.groupby("time.year").sum("time")
annual_prec

In [55]:

# plot annual precipitation for selected year: 2018
prec = annual_prec.precip
year2018 = prec.sel(year=2018)

fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines() 
year2018.plot()
plt.show()

In [56]:

# extract annual precipitation for a specific location: vanocuver BC
annual_prec_loc = annual_prec.sel(lat = 49, 
                                  lon = -123+360,
                                  year = [x for x in range(1990, 2020)],
                                  
                                   method='nearest')

annual_prec_loc

In [57]:

# plot results
annual_prec_loc.precip.plot()

mean_val = annual_prec_loc.mean('year')['precip']
plt.axhline (mean_val, color='g', linestyle='dashed') # plot the mean value

fig = plt.gcf()
fig.set_size_inches(18, 6)

plt.show

Out[57]:

<function matplotlib.pyplot.show(close=None, block=None)>

In [58]:

# calculate Annual Precipitation anomaly (z_score) for each year
mean = annual_prec.mean('year')['precip']
stdev = annual_prec.std('year')['precip']

annual_prec['z_score'] = (annual_prec.precip - mean) / stdev
annual_prec

In [59]:

# plot the z_score for selected year: 2018
z_score = annual_prec.z_score
year2018 = z_score.sel(year=2018)

fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
ax = plt.axes(projection=ccrs.PlateCarree())
year2018.plot()
plt.show()

In [60]:

# Histogam of z_score
annual_prec.z_score.plot()

fig = plt.gcf()
fig.set_size_inches(18, 6)

In [61]:

# save the new xarray dataset to netCDF
output = r'\\...\precip_annual_anomalies.nc' 
annual_prec.to_netcdf(path=output, mode='w')

In [ ]: