Remap geothermal dataset (Davies, 2013) to tx2_3

%matplotlib inline

import xarray as xr
import xesmf, intake
import numpy as np
import datetime
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

today = datetime.date.today().strftime("%y%m%d")
print(today)

260305

fname = '../mesh/tx2_3v3_grid.nc'
ds_out = xr.open_dataset(fname).rename({'tlon': 'lon','tlat': 'lat', 'qlon': 'lon_b',
                                        'qlat': 'lat_b', 'nx' : 'xh', 'ny' : 'yh',
                                        'depth' : 'z_l'})
ds_out

<xarray.Dataset> Size: 42MB
Dimensions:  (yh: 480, xh: 540, nxp: 541, nyp: 481)
Dimensions without coordinates: yh, xh, nxp, nyp
Data variables: (12/20)
    lon      (yh, xh) float64 2MB ...
    lat      (yh, xh) float64 2MB ...
    ulon     (yh, nxp) float64 2MB ...
    ulat     (yh, nxp) float64 2MB ...
    vlon     (nyp, xh) float64 2MB ...
    vlat     (nyp, xh) float64 2MB ...
    ...       ...
    tarea    (yh, xh) float64 2MB ...
    tmask    (yh, xh) float64 2MB ...
    angle    (yh, xh) float64 2MB ...
    z_l      (yh, xh) float64 2MB ...
    ar       (yh, xh) float64 2MB ...
    egs      (yh, xh) float64 2MB ...
Attributes:
    Description:  CESM MOM6 2/3 degree grid
    Author:       Frank, Fred, Gustavo (gmarques@ucar.edu)
    Created:      2026-03-05T14:49:28.971877
    type:         Glogal 2/3 degree grid file

xarray.Dataset

• Dimensions:
  • yh: 480
  • xh: 540
  • nxp: 541
  • nyp: 481
• Data variables: (20)
  • lon
    (yh, xh)
    float64
    ...

    long_name :

    array of t-grid longitudes

    units :

    degrees_east

    [259200 values with dtype=float64]

  • lat
    (yh, xh)
    float64
    ...

    long_name :

    array of t-grid latitudes

    units :

    degrees_north

    [259200 values with dtype=float64]

  • ulon
    (yh, nxp)
    float64
    ...

    long_name :

    array of u-grid longitudes

    units :

    degrees_east

    [259680 values with dtype=float64]

  • ulat
    (yh, nxp)
    float64
    ...

    long_name :

    array of u-grid latitudes

    units :

    degrees_north

    [259680 values with dtype=float64]

  • vlon
    (nyp, xh)
    float64
    ...

    long_name :

    array of v-grid longitudes

    units :

    degrees_east

    [259740 values with dtype=float64]

  • vlat
    (nyp, xh)
    float64
    ...

    long_name :

    array of v-grid latitudes

    units :

    degrees_north

    [259740 values with dtype=float64]

  • lon_b
    (nyp, nxp)
    float64
    ...

    long_name :

    array of q-grid longitudes

    units :

    degrees_east

    [260221 values with dtype=float64]

  • lat_b
    (nyp, nxp)
    float64
    ...

    long_name :

    array of q-grid latitudes

    units :

    degrees_north

    [260221 values with dtype=float64]

  • dxt
    (yh, xh)
    float64
    ...

    long_name :

    x-distance between u-points, centered at t

    units :

    meters

    [259200 values with dtype=float64]

  • dyt
    (yh, xh)
    float64
    ...

    long_name :

    y-distance between v-points, centered at t

    units :

    meters

    [259200 values with dtype=float64]

  • dxCv
    (yh, xh)
    float64
    ...

    long_name :

    x-distance between q-points, centered at v

    units :

    meters

    [259200 values with dtype=float64]

  • dyCu
    (yh, xh)
    float64
    ...

    long_name :

    y-distance between q-points, centered at u

    units :

    meters

    [259200 values with dtype=float64]

  • dxCu
    (yh, xh)
    float64
    ...

    long_name :

    x-distance between t-points, centered at u

    units :

    meters

    [259200 values with dtype=float64]

  • dyCv
    (yh, xh)
    float64
    ...

    long_name :

    y-distance between t-points, centered at v

    units :

    meters

    [259200 values with dtype=float64]

  • tarea
    (yh, xh)
    float64
    ...

    long_name :

    area of t-cells

    units :

    meters^2

    [259200 values with dtype=float64]

  • tmask
    (yh, xh)
    float64
    ...

    long_name :

    ocean fraction at t-cell centers

    units :

    none

    [259200 values with dtype=float64]

  • angle
    (yh, xh)
    float64
    ...

    long_name :

    angle grid makes with latitude line

    units :

    degrees

    [259200 values with dtype=float64]

  • z_l
    (yh, xh)
    float64
    ...

    long_name :

    depth at h points

    units :

    meters

    [259200 values with dtype=float64]

  • ar
    (yh, xh)
    float64
    ...

    long_name :

    grid aspect ratio (dyt/dxt)

    units :

    none

    [259200 values with dtype=float64]

  • egs
    (yh, xh)
    float64
    ...

    long_name :

    grid effective grid spacing

    units :

    degrees

    [259200 values with dtype=float64]

• Attributes: (4)

  Description :

  CESM MOM6 2/3 degree grid

  Author :

  Frank, Fred, Gustavo (gmarques@ucar.edu)

  Created :

  2026-03-05T14:49:28.971877

  type :

  Glogal 2/3 degree grid file

infile = '/glade/campaign/cgd/oce/datasets/obs/geothermal/geothermal_davies2013_ncar.nc'
ds_in = xr.open_dataset(infile).rename({'tlon': 'lon','tlat': 'lat', 'qlon': 'lon_b',
                                        'qlat': 'lat_b'})
ds_in

<xarray.Dataset> Size: 56MB
Dimensions:        (yh: 1080, xh: 1440, yq: 1081, xq: 1441)
Coordinates:
  * yh             (yh) float32 4kB 0.5 1.5 2.5 ... 1.078e+03 1.078e+03 1.08e+03
  * xh             (xh) float32 6kB 0.5 1.5 2.5 ... 1.438e+03 1.438e+03 1.44e+03
Dimensions without coordinates: yq, xq
Data variables:
    geothermal_hf  (yh, xh) float32 6MB ...
    lat            (yh, xh) float64 12MB ...
    lon            (yh, xh) float64 12MB ...
    lat_b          (yq, xq) float64 12MB ...
    lon_b          (yq, xq) float64 12MB ...
Attributes:
    title:          Geothermal heat flow from Davies, 2013, re-gridded to OM4...
    reference:      Davies, J. Huw, 2013: Global map of solid Earth surface h...
    reference_url:  http://dx.doi.org/10.1002/ggge.20271

xarray.Dataset

• Dimensions:
  • yh: 1080
  • xh: 1440
  • yq: 1081
  • xq: 1441
• Coordinates: (2)
  • yh
    (yh)
    float32
    0.5 1.5 2.5 ... 1.078e+03 1.08e+03

    long_name :

    Grid j-index

    array([5.0000e-01, 1.5000e+00, 2.5000e+00, ..., 1.0775e+03, 1.0785e+03,
           1.0795e+03], shape=(1080,), dtype=float32)

  • xh
    (xh)
    float32
    0.5 1.5 2.5 ... 1.438e+03 1.44e+03

    long_name :

    Grid i-index

    array([5.0000e-01, 1.5000e+00, 2.5000e+00, ..., 1.4375e+03, 1.4385e+03,
           1.4395e+03], shape=(1440,), dtype=float32)

• Data variables: (5)
  • geothermal_hf
    (yh, xh)
    float32
    ...

    long_name :

    Geothermal heat flow

    standard_name :

    upward_geothermal_heat_flux_at_sea_floor

    units :

    W m-2

    cell_methods :

    area: mean

    sha1 :

    519321c160a5495abde42aebe83c7aa525eac40f

    [1555200 values with dtype=float32]

  • lat
    (yh, xh)
    float64
    ...

    description :

    Latitude of tracer (T) points

    units :

    degrees_north

    [1555200 values with dtype=float64]

  • lon
    (yh, xh)
    float64
    ...

    description :

    Longitude of tracer (T) points

    units :

    degrees_east

    [1555200 values with dtype=float64]

  • lat_b
    (yq, xq)
    float64
    ...

    description :

    Latitude of corner (Bu) points

    units :

    degrees_north

    [1557721 values with dtype=float64]

  • lon_b
    (yq, xq)
    float64
    ...

    description :

    Longitude of corner (Bu) points

    units :

    degrees_east

    [1557721 values with dtype=float64]

• Attributes: (3)

  title :

  Geothermal heat flow from Davies, 2013, re-gridded to OM4_025

  reference :

  Davies, J. Huw, 2013: Global map of solid Earth surface heat flow. Geochemistry, Geophysics, Geosystems, 14 (10), pp 4608--4622. doi:10.1002/ggge.20271

  reference_url :

  http://dx.doi.org/10.1002/ggge.20271

def regrid_tracer(fld, ds_in, ds_out, method='conservative'):

    regrid = xesmf.Regridder(
        ds_in,
        ds_out,
        method=method,
        periodic=True,
    )
    fld_out = regrid(ds_in[fld])
    return fld_out

Potential temperature

geothermal_hf = regrid_tracer('geothermal_hf', ds_in, ds_out)
geothermal_hf

<xarray.DataArray (yh: 480, xh: 540)> Size: 1MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
        0.05329704],
       [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
        0.05336878],
       [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
        0.05337462]], shape=(480, 540), dtype=float32)
Dimensions without coordinates: yh, xh
Attributes:
    regrid_method:  conservative

xarray.DataArray

• yh: 480
• xh: 540

• 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.05587 0.05508 0.05432 0.05364 0.05337

  array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         ...,
         [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
          0.05329704],
         [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
          0.05336878],
         [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
          0.05337462]], shape=(480, 540), dtype=float32)

• Attributes: (1)

  regrid_method :

  conservative

geothermal_hf = geothermal_hf.rename('geothermal_hf')

Visual inspection

Make sure original and remapped plots look similar.

# visual inspection. Make sure original and remapped plots look similar
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(18,4))
geothermal_hf.plot.pcolormesh(ax=axes[0], vmin=0,vmax=2)
ds_in['geothermal_hf'].plot.pcolormesh(ax=axes[1], vmin=0,vmax=2)
axes[0].set_title('Remapped')
axes[1].set_title('Original grid')
plt.suptitle('geothermal_hf')

Text(0.5, 0.98, 'geothermal_hf')

data = xr.DataArray(geothermal_hf.data, dims=('yh', 'xh'), name='geothermal_hf')
data

<xarray.DataArray 'geothermal_hf' (yh: 480, xh: 540)> Size: 1MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
        0.05329704],
       [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
        0.05336878],
       [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
        0.05337462]], shape=(480, 540), dtype=float32)
Dimensions without coordinates: yh, xh

xarray.DataArray

'geothermal_hf'

• yh: 480
• xh: 540

• 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.05587 0.05508 0.05432 0.05364 0.05337

  array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         ...,
         [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
          0.05329704],
         [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
          0.05336878],
         [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
          0.05337462]], shape=(480, 540), dtype=float32)

data

<xarray.DataArray 'geothermal_hf' (yh: 480, xh: 540)> Size: 1MB
array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
        0.        ],
       ...,
       [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
        0.05329704],
       [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
        0.05336878],
       [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
        0.05337462]], shape=(480, 540), dtype=float32)
Dimensions without coordinates: yh, xh

xarray.DataArray

'geothermal_hf'

• yh: 480
• xh: 540

• 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.05587 0.05508 0.05432 0.05364 0.05337

  array([[0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         [0.        , 0.        , 0.        , ..., 0.        , 0.        ,
          0.        ],
         ...,
         [0.05298106, 0.05332395, 0.05414148, ..., 0.05428809, 0.05358988,
          0.05329704],
         [0.05306507, 0.05335549, 0.05416296, ..., 0.05430644, 0.05361681,
          0.05336878],
         [0.0530719 , 0.05337847, 0.05417501, ..., 0.05431673, 0.05363644,
          0.05337462]], shape=(480, 540), dtype=float32)

data = xr.DataArray(geothermal_hf.data, dims=('yh', 'xh'), name='geothermal_hf')

# Global attrs
data.attrs['long_name'] = 'Geothermal heat flow'
data.attrs['standard_name'] = 'upward_geothermal_heat_flux_at_sea_floor'
data.attrs['units'] = 'W m-2'
data.attrs['cell_methods'] = 'area: mean'

# Create a dataset and add the DataArray to it
dataset = xr.Dataset({'geothermal_hf': data})

dataset.attrs['description'] = 'Geothermal heat flow from Davies, 2013, re-gridded to tx2_3'
dataset.attrs['author'] = 'Gustavo Marques (gmarques@ucar.edu)'
dataset.attrs['date'] = today
dataset.attrs['infile'] = infile
dataset.attrs['url'] = 'https://github.com/NCAR/tx2_3/geothermal'
# save
fname = 'geothermal_davies2013_tx2_3v3_{}.nc'.format(today)

dataset.to_netcdf(fname, 
                  engine="netcdf4",
                  format="NETCDF3_64BIT_DATA",
                  encoding={'geothermal_hf': {'_FillValue': None}})