Use case: writing a binary input file for MITgcm

You may want to write binary files to create forcing data, initial condition,... for your MITgcm configuration. Here we show how xmitgcm can help.

Simple case: a regular grid



In [1]:

    
import numpy as np
import xmitgcm
import matplotlib.pylab as plt

Let's build a regular lat/lon grid with one degree resolution and create a pseudo-field on this regular grid:



In [2]:

    
lon = np.arange(-180,180,1)
lat = np.arange(-90,90,1)

lon2, lat2 = np.meshgrid(lon,lat)

pseudo = np.cos(2*np.pi*lat2/360) * np.cos(4*np.pi*np.pi*lon2*lon2/360/360)



In [3]:

    
plt.contourf(lon2, lat2, pseudo)
plt.colorbar()









    Out[3]:





<matplotlib.colorbar.Colorbar at 0x7ff5add1c668>

We can write the field as a binary file, to be used as an input file for the model with xmitgcm.utils.write_to_binary. Default is single precision, but double precision can be written with corresponding numpy.dtype. Note that here we use a numpy.array but one can use xarray as well using the DataArray.values



In [4]:

    
xmitgcm.utils.write_to_binary?









    





Signature: xmitgcm.utils.write_to_binary(flatdata, fileout, dtype=dtype('float32'))
Docstring:
write data in binary file

PARAMETERS:

flatdata: numpy.array
    vector of data to write
fileout: str
    output file name
dtype: np.dtype
    single/double precision
File:      ~/.conda/envs/production2/lib/python3.6/site-packages/xmitgcm-0.2.2-py3.6.egg/xmitgcm/utils.py
Type:      function



In [5]:

    
xmitgcm.utils.write_to_binary(pseudo.flatten(), 'file1.bin')

More complicated case: a LLC grid

In this case, let's assume we have a xarray dataarray or dataset well formatted on the llc grid. This dataset can be the result of a regridding onto the LLC grid that we want to use as an initial condition for the model (for example). We need to generate the binary file that MITgcm can read. Here's what we can do:



In [6]:

    
# First let's download a sample dataset
! wget https://ndownloader.figshare.com/files/14066567
! tar -xf 14066567









    



--2019-05-01 21:16:09--  https://ndownloader.figshare.com/files/14066567
Resolving ndownloader.figshare.com (ndownloader.figshare.com)... 52.48.232.64, 63.32.41.137, 54.229.125.140, ...
Connecting to ndownloader.figshare.com (ndownloader.figshare.com)|52.48.232.64|:443... connected.
HTTP request sent, awaiting response... 302 Found
Location: https://s3-eu-west-1.amazonaws.com/pfigshare-u-files/14066567/global_oce_llc90.tar.gz [following]
--2019-05-01 21:16:09--  https://s3-eu-west-1.amazonaws.com/pfigshare-u-files/14066567/global_oce_llc90.tar.gz
Resolving s3-eu-west-1.amazonaws.com (s3-eu-west-1.amazonaws.com)... 52.218.36.2
Connecting to s3-eu-west-1.amazonaws.com (s3-eu-west-1.amazonaws.com)|52.218.36.2|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 84966646 (81M) [application/gzip]
Saving to: ‘14066567.1’

100%[======================================>] 84,966,646  12.9MB/s   in 6.6s   

2019-05-01 21:16:16 (12.3 MB/s) - ‘14066567.1’ saved [84966646/84966646]

We can load this dataset with xmitgcm:



In [7]:

    
extra_metadata = xmitgcm.utils.get_extra_metadata(domain='llc', nx=90)

ds = xmitgcm.open_mdsdataset('./global_oce_llc90/', iters= [8], geometry='llc', extra_metadata=extra_metadata)



In [8]:

    
ds









    Out[8]:





<xarray.Dataset>
Dimensions:   (face: 13, i: 90, i_g: 90, j: 90, j_g: 90, k: 50, k_l: 50, k_p1: 51, k_u: 50, time: 1)
Coordinates:
  * i         (i) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
  * i_g       (i_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * j         (j) int64 0 1 2 3 4 5 6 7 8 9 10 ... 80 81 82 83 84 85 86 87 88 89
  * j_g       (j_g) int64 0 1 2 3 4 5 6 7 8 9 ... 80 81 82 83 84 85 86 87 88 89
  * k         (k) int64 0 1 2 3 4 5 6 7 8 9 10 ... 40 41 42 43 44 45 46 47 48 49
  * k_u       (k_u) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_l       (k_l) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
  * k_p1      (k_p1) int64 0 1 2 3 4 5 6 7 8 9 ... 41 42 43 44 45 46 47 48 49 50
  * face      (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
    iter      (time) int64 8
  * time      (time) int64 8
    XC        (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    YC        (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    XG        (face, j_g, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    YG        (face, j_g, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    CS        (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    SN        (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    Z         (k) >f4 -5.0 -15.0 -25.0 -35.0 ... -5039.25 -5461.25 -5906.25
    Zp1       (k_p1) >f4 0.0 -10.0 -20.0 -30.0 ... -5244.5 -5678.0 -6134.5
    Zu        (k_u) >f4 -10.0 -20.0 -30.0 -40.0 ... -5244.5 -5678.0 -6134.5
    Zl        (k_l) >f4 0.0 -10.0 -20.0 -30.0 ... -4834.0 -5244.5 -5678.0
    rA        (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    dxG       (face, j_g, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    dyG       (face, j, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    Depth     (face, j, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    rAz       (face, j_g, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    dxC       (face, j, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    dyC       (face, j_g, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    rAw       (face, j, i_g) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    rAs       (face, j_g, i) >f4 dask.array<shape=(13, 90, 90), chunksize=(1, 90, 90)>
    drC       (k_p1) >f4 5.0 10.0 10.0 10.0 10.0 ... 399.0 422.0 445.0 228.25
    drF       (k) >f4 10.0 10.0 10.0 10.0 10.0 ... 364.5 387.5 410.5 433.5 456.5
    PHrefC    (k) >f4 49.05 147.15 245.25 ... 49435.043 53574.863 57940.312
    PHrefF    (k_p1) >f4 0.0 98.1 196.2 294.3 ... 51448.547 55701.18 60179.445
    hFacC     (k, face, j, i) >f4 dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    hFacW     (k, face, j, i_g) >f4 dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    hFacS     (k, face, j_g, i) >f4 dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    maskC     (k, face, j, i) bool dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    maskW     (k, face, j, i_g) bool dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    maskS     (k, face, j_g, i) bool dask.array<shape=(50, 13, 90, 90), chunksize=(1, 1, 90, 90)>
Data variables:
    W         (time, k_l, face, j, i) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
    PHL       (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    PH        (time, k, face, j, i) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
    V         (time, k, face, j_g, i) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
    Eta       (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    S         (time, k, face, j, i) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
    ETAN      (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIarea    (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIheff    (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIhsnow   (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    DETADT2   (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    PHIBOT    (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    sIceLoad  (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    MXLDEPTH  (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    oceSPDep  (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIatmQnt  (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIatmFW   (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    oceQnet   (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    oceFWflx  (time, face, j, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    oceTAUX   (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    oceTAUY   (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    ADVxHEFF  (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    ADVyHEFF  (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    DFxEHEFF  (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    DFyEHEFF  (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    ADVxSNOW  (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    ADVySNOW  (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    DFxESNOW  (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    DFyESNOW  (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIuice    (time, face, j, i_g) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    SIvice    (time, face, j_g, i) >f4 dask.array<shape=(1, 13, 90, 90), chunksize=(1, 1, 90, 90)>
    U         (time, k, face, j, i_g) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
    T         (time, k, face, j, i) >f4 dask.array<shape=(1, 50, 13, 90, 90), chunksize=(1, 1, 1, 90, 90)>
Attributes:
    Conventions:  CF-1.6
    title:        netCDF wrapper of MITgcm MDS binary data
    source:       MITgcm
    history:      Created by calling `open_mdsdataset(extra_metadata={'has_fa...

Now let's say we want to use the temperature T and make it the initial condition for another simulation. First we need to re-build the facets, concatenate them into the compact form that MITgcm reads/writes and then write the compact to a binary file. We would do this as follows:



In [9]:

    
# temperature
facets = xmitgcm.utils.rebuild_llc_facets(ds['T'], extra_metadata)
compact = xmitgcm.utils.llc_facets_3d_spatial_to_compact(facets, 'Z', extra_metadata)
xmitgcm.utils.write_to_binary(compact, 'T_initial_condition.bin')

In this case, we already had the binary file to read from so we can compare checksums:



In [10]:

    
!md5sum T_initial_condition.bin









    



045086cc85a66aa59a713ab7d716539e  T_initial_condition.bin



In [11]:

    
!md5sum ./global_oce_llc90/T.0000000008.data









    



045086cc85a66aa59a713ab7d716539e  ./global_oce_llc90/T.0000000008.data

The file generated is identical to the original file! A similar function exist for 2d files,...