In [1]:

import numpy as np
import dask
import xarray as xr
import matplotlib.pyplot as plt
import os
from os.path import isfile
import datetime as dt
from salishsea_tools import evaltools as et, viz_tools
import glob
import dask
from dask.distributed import Client, progress
import dask.array as da
%matplotlib inline

In [2]:

meshPath='/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702_noLPE.nc'
maskName='tmask'

In [3]:

j0=230;j1=470;i0=0;i1=200;
ct=240;cz=40;cy=50;cx=200;

In [4]:

def prepro(ds):
    return ds.loc[dict(y=slice(j0,j1),x=slice(i0,i1))]

In [5]:

mesh=xr.open_mfdataset(meshPath,chunks={'t':1, 'z':cz, 'y':cy, 'x':cx},preprocess=prepro,parallel=True)

/home/eolson/anaconda3/envs/py37/lib/python3.7/site-packages/ipykernel/__main__.py:1: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`
will change. To retain the existing behavior, pass
combine='nested'. To use future default behavior, pass
combine='by_coords'. See
http://xarray.pydata.org/en/stable/combining.html#combining-multi

  if __name__ == '__main__':
/home/eolson/.local/lib/python3.7/site-packages/xarray/backends/api.py:933: FutureWarning: The datasets supplied have global dimension coordinates. You may want
to use the new `combine_by_coords` function (or the
`combine='by_coords'` option to `open_mfdataset`) to order the datasets
before concatenation. Alternatively, to continue concatenating based
on the order the datasets are supplied in future, please use the new
`combine_nested` function (or the `combine='nested'` option to
open_mfdataset).
  from_openmfds=True,

In [6]:

mesh2=mesh.loc[{'z':0}]

In [7]:

plt.pcolormesh(mesh2['tmask'].squeeze())

Out[7]:

<matplotlib.collections.QuadMesh at 0x7f01f83154a8>

In [8]:

flistmuZ=et.index_model_files(dt.datetime(2015,6,1),dt.datetime(2015,9,1),
                             '/data/eolson/results/MEOPAR/SS36runs/CedarRuns/testmuZ/',
                             'long',10,'ptrc_T',1)

first file starts on  2015-05-31 00:00:00

In [9]:

flistHC=et.index_model_files(flistmuZ['t_0'][0],flistmuZ['t_n'][len(flistmuZ)-1]-dt.timedelta(days=1),
                             '/results/SalishSea/hindcast.201905/',
                             'nowcast',1,'ptrc_T',1)

In [10]:

client = dask.distributed.Client(
    n_workers=2, threads_per_worker=2, processes=True)
client

Out[10]:

Client

Scheduler: tcp://127.0.0.1:45792
Dashboard: http://127.0.0.1:8787/status

Cluster

Workers: 2
Cores: 4
Memory: 16.77 GB

In [11]:

chunks = {"time_counter": 3, "deptht": 40*3, "y": 898*3, "x": 398*3}

In [12]:

t0=dt.datetime.now()
dataHC=xr.open_mfdataset(flistHC['paths'],data_vars=('diatoms','microzooplankton'),chunks=chunks,
            drop_variables=('nitrate','ammonium','silicon','flagellates','ciliates','dissolved_organic_nitrogen',
                            'particulate_organic_nitrogen','biogenic_silicon','mesozooplankton'),preprocess=prepro,parallel=True)
t1=dt.datetime.now()
print(t1-t0)

/home/eolson/anaconda3/envs/py37/lib/python3.7/site-packages/ipykernel/__main__.py:4: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`
will change. To retain the existing behavior, pass
combine='nested'. To use future default behavior, pass
combine='by_coords'. See
http://xarray.pydata.org/en/stable/combining.html#combining-multi

/home/eolson/.local/lib/python3.7/site-packages/xarray/backends/api.py:933: FutureWarning: The datasets supplied have global dimension coordinates. You may want
to use the new `combine_by_coords` function (or the
`combine='by_coords'` option to `open_mfdataset`) to order the datasets
before concatenation. Alternatively, to continue concatenating based
on the order the datasets are supplied in future, please use the new
`combine_nested` function (or the `combine='nested'` option to
open_mfdataset).
  from_openmfds=True,

0:01:29.028540

In [13]:

dataHC

Out[13]:

<xarray.Dataset>
Dimensions:               (axis_nbounds: 2, deptht: 40, nvertex: 4, time_counter: 2400, x: 200, y: 240)
Coordinates:
    nav_lon               (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
    nav_lat               (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
  * deptht                (deptht) float32 0.5000003 1.5000031 ... 441.4661
    time_centered         (time_counter) datetime64[ns] dask.array<chunksize=(3,), meta=np.ndarray>
  * time_counter          (time_counter) datetime64[ns] 2015-05-31T00:30:00 ... 2015-09-07T23:30:00
Dimensions without coordinates: axis_nbounds, nvertex, x, y
Data variables:
    bounds_lat            (y, x, nvertex) float32 dask.array<chunksize=(240, 200, 4), meta=np.ndarray>
    deptht_bounds         (deptht, axis_nbounds) float32 dask.array<chunksize=(40, 2), meta=np.ndarray>
    area                  (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
    bounds_lon            (y, x, nvertex) float32 dask.array<chunksize=(240, 200, 4), meta=np.ndarray>
    time_centered_bounds  (time_counter, axis_nbounds) datetime64[ns] dask.array<chunksize=(3, 2), meta=np.ndarray>
    time_counter_bounds   (time_counter, axis_nbounds) datetime64[ns] dask.array<chunksize=(3, 2), meta=np.ndarray>
    diatoms               (time_counter, deptht, y, x) float32 dask.array<chunksize=(3, 40, 240, 200), meta=np.ndarray>
    microzooplankton      (time_counter, deptht, y, x) float32 dask.array<chunksize=(3, 40, 240, 200), meta=np.ndarray>
Attributes:
    name:         SalishSea_1h_20150526_20150531
    description:  biogeochemical variables
    title:        biogeochemical variables
    Conventions:  CF-1.6
    timeStamp:    2019-Sep-25 01:38:45 GMT
    uuid:         7abb4fe2-2cf5-475e-b3c2-b5200adfd953

In [14]:

t0=dt.datetime.now()
datamuZ=xr.open_mfdataset(flistmuZ['paths'],data_vars=('diatoms','microzooplankton'),chunks=chunks,
            drop_variables=('nitrate','ammonium','silicon','flagellates','ciliates','dissolved_organic_nitrogen',
                            'particulate_organic_nitrogen','biogenic_silicon','mesozooplankton'),preprocess=prepro,parallel=True)
t1=dt.datetime.now()
print(t1-t0)

/home/eolson/anaconda3/envs/py37/lib/python3.7/site-packages/ipykernel/__main__.py:4: FutureWarning: In xarray version 0.15 the default behaviour of `open_mfdataset`
will change. To retain the existing behavior, pass
combine='nested'. To use future default behavior, pass
combine='by_coords'. See
http://xarray.pydata.org/en/stable/combining.html#combining-multi

/home/eolson/.local/lib/python3.7/site-packages/xarray/backends/api.py:933: FutureWarning: The datasets supplied have global dimension coordinates. You may want
to use the new `combine_by_coords` function (or the
`combine='by_coords'` option to `open_mfdataset`) to order the datasets
before concatenation. Alternatively, to continue concatenating based
on the order the datasets are supplied in future, please use the new
`combine_nested` function (or the `combine='nested'` option to
open_mfdataset).
  from_openmfds=True,

0:00:07.497164

In [15]:

datamuZ=datamuZ.rename({'time_counter':'t','deptht':'z'})

In [16]:

datamuZ

Out[16]:

<xarray.Dataset>
Dimensions:               (axis_nbounds: 2, nvertex: 4, t: 2400, x: 200, y: 240, z: 40)
Coordinates:
    nav_lon               (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
    nav_lat               (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
  * z                     (z) float32 0.5000003 1.5000031 ... 414.5341 441.4661
    time_centered         (t) datetime64[ns] dask.array<chunksize=(3,), meta=np.ndarray>
  * t                     (t) datetime64[ns] 2015-05-31T00:30:00 ... 2015-09-07T23:30:00
Dimensions without coordinates: axis_nbounds, nvertex, x, y
Data variables:
    deptht_bounds         (z, axis_nbounds) float32 dask.array<chunksize=(40, 2), meta=np.ndarray>
    area                  (y, x) float32 dask.array<chunksize=(240, 200), meta=np.ndarray>
    bounds_nav_lon        (y, x, nvertex) float32 dask.array<chunksize=(240, 200, 4), meta=np.ndarray>
    bounds_nav_lat        (y, x, nvertex) float32 dask.array<chunksize=(240, 200, 4), meta=np.ndarray>
    time_centered_bounds  (t, axis_nbounds) datetime64[ns] dask.array<chunksize=(3, 2), meta=np.ndarray>
    time_counter_bounds   (t, axis_nbounds) datetime64[ns] dask.array<chunksize=(3, 2), meta=np.ndarray>
    diatoms               (t, z, y, x) float32 dask.array<chunksize=(3, 40, 240, 200), meta=np.ndarray>
    microzooplankton      (t, z, y, x) float32 dask.array<chunksize=(3, 40, 240, 200), meta=np.ndarray>
Attributes:
    name:         SalishSea_1h_20150531_20150709_ptrc_T
    description:  biogeochemical variables
    title:        biogeochemical variables
    Conventions:  CF-1.6
    timeStamp:    2019-Oct-30 15:36:14 GMT
    uuid:         396ea30f-075b-4210-9e66-197509446563

In [ ]:

In [17]:

tmask=da.asanyarray(mesh['tmask'])
e3t=da.asanyarray(mesh['e3t_0'])
diatomsmuZ=da.asanyarray(datamuZ['diatoms'])
diatomsHC=da.asanyarray(dataHC['diatoms'])

In [18]:

tmask

Out[18]:

	Array	Chunk
Bytes	1.92 MB	400.00 kB
Shape	(1, 40, 240, 200)	(1, 40, 50, 200)
Count	43 Tasks	6 Chunks
Type	int8	numpy.ndarray

In [19]:

e3t

Out[19]:

	Array	Chunk
Bytes	15.36 MB	3.20 MB
Shape	(1, 40, 240, 200)	(1, 40, 50, 200)
Count	43 Tasks	6 Chunks
Type	float64	numpy.ndarray

In [20]:

diatomsmuZ

Out[20]:

	Array	Chunk
Bytes	18.43 GB	23.04 MB
Shape	(2400, 40, 240, 200)	(3, 40, 240, 200)
Count	2410 Tasks	800 Chunks
Type	float32	numpy.ndarray

In [21]:

prod1=tmask*e3t*diatomsmuZ

In [22]:

prod1

Out[22]:

	Array	Chunk
Bytes	36.86 GB	9.60 MB
Shape	(2400, 40, 240, 200)	(3, 40, 50, 200)
Count	17702 Tasks	4800 Chunks
Type	float64	numpy.ndarray

In [23]:

sum1=dask.array.sum(prod1,axis=1)

In [24]:

sum1

Out[24]:

	Array	Chunk
Bytes	921.60 MB	240.00 kB
Shape	(2400, 240, 200)	(3, 50, 200)
Count	27302 Tasks	4800 Chunks
Type	float64	numpy.ndarray

In [25]:

diatomsmuZInt=da.mean(sum1,0)

In [26]:

prod2=tmask*e3t*diatomsHC
sum2=dask.array.sum(prod2,axis=1)
diatomsHCInt=da.mean(sum2,0)

In [27]:

#diatomsmuZInt.visualize()

In [28]:

t0=dt.datetime.now()
diatomsmuZInt=diatomsmuZInt.compute()
t1=dt.datetime.now()
print(t1-t0)

0:22:01.200207

In [29]:

t0=dt.datetime.now()
diatomsHCInt=diatomsHCInt.compute()
t1=dt.datetime.now()
print(t1-t0)

0:16:50.440003

In [30]:

t0=dt.datetime.now()
fig,ax=plt.subplots(1,3,figsize=(15,4))
m0=ax[0].pcolormesh(diatomsmuZInt)
plt.colorbar(m0,ax=ax[0])
ax[0].set_title('muZ Diatoms')
m1=ax[1].pcolormesh(diatomsHCInt)
plt.colorbar(m1,ax=ax[1])
ax[1].set_title('HC Diatoms')
m2=ax[2].pcolormesh(diatomsmuZInt-diatomsHCInt)
plt.colorbar(m2,ax=ax[2])
ax[2].set_title('Diff')
t1=dt.datetime.now()
print(t1-t0)

0:00:00.313700

In [31]:

client.close()