Example 5 - Mapping Curie depth from EMAG2

Using global compilations of the magnetic anomaly, such as EMAG2, means there is no limitation on window sizes. The mapping module in PyCurious has several functions to simplify:

  1. Translating between Coordinate Reference Systems (CRS)

  2. Griding scattered point data

  3. Importing and exporting Geotiffs

These can be accessed from pycurious.mapping

This requires some extra dependencies:

  • pyproj

  • pyepsg

  • cartopy (for visualisation)

which can be installed via pip

pip install [--user] pyproj pyepsg cartopy

Contents

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

import pycurious
from pycurious import mapping

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader

Import EMAG2 data

The EMAG2 (v3) Earth magnetic anomaly grid\(^1\) has been reformatted from a CSV file to a compressed .npz NumPy archive to efficiently load within Python workflows. This dataset can be downloaded from Zenodo. We have added some download tools accessed under pycurious.download to ease the download process and checksum validation.

References

\(^1\) Brian Meyer, Richard Saltus, Arnaud Chulliat (2017): EMAG2: Earth Magnetic Anomaly Grid (2-arc-minute resolution) Version 3. National Centers for Environmental Information, NOAA. Model. doi:10.7289/V5H70CVX

file_resources = [
    # EMAG2 v3
    {
     "local_file":"../../data/EMAG2_V3_20170530.npz",
     "md5":'c0898b6a91efb3f13783873a8b67380c',
     "url":"https://zenodo.org/record/3245551/files/EMAG2_V3_20170530.npz?download=1",
     "expected_size":"500Mb"
    },
    
    # Li et al. 2017
    {
     "local_file":"../../data/Li_et_al_2017.txt",
     "md5":'5f0ea0af3e27c6c21cd7b776c979aa80',
     "url":"https://static-content.springer.com/esm/art%3A10.1038%2Fsrep45129/MediaObjects/41598_2017_BFsrep45129_MOESM71_ESM.txt",
     "expected_size":"15Mb"
    },    
]

from pycurious.download import download_cached_file

for resource in file_resources:
    print ("\nDownloading {:s}".format(resource["local_file"]))
    download_cached_file(resource["url"], resource["local_file"], resource["md5"], resource["expected_size"])

Once loaded, we trim the data to a region of interest - being careful to allow some buffer to accommodate the maximum window size.

EMAG2_V3_20170530.npz contains the following columns:

  1. LON ; Geographic Longitude WGS84 (decimal degrees)

  2. LAT ; Geographic Latitude WGS84 (decimal degrees)

  3. SeaLevel ; Magnetic Anomaly Value at Sea Level(nT)

  4. UpCont ; Magnetic Anomaly Value at continuous 4km altitude (nT)

  5. Code ; Data Source Code (see table below)

  6. Error ; Error estimate (nT)

Code 888 is assigned in certain cells on grid edges where the data source is ambiguous and assigned an error of -888 nT

Code 999 is assigned in cells where no data is reported with the anomaly value assigned 99999 nT and an error of -999 nT

emag2_file = file_resources[0]["local_file"]
with np.load(emag2_file) as f:
    mag_data = f['data']

# filter NaNs
mag_data = mag_data[mag_data[:,3] != 99999.]
mag_data = mag_data[mag_data[:,4] != -888.]
mag_data = mag_data[mag_data[:,4] != -999.]

# print min/max
mincols = mag_data.min(axis=0)
maxcols = mag_data.max(axis=0)

fmt = "min/max {:5.2f} -> {:5.2f}"
for col in range(mag_data.shape[1]):
    print(fmt.format(mincols[col], maxcols[col]))

Interpolate onto grid projection

We choose an arbitrary section of the Earth’s surface and a suitable projection. PyCurious handles transformation between different CRS using a EPSG reference code. EMAG2 is in lons and lats (in decimal degrees), which we want to convert into eastings and northings (in metres). For this example, we transform WGS84 coordinates (EPSG:4326) to the IRENET95 projection (EPSG:2157).

IMPORTANT: The power spectrum must be in eastings & northings to compute the power spectrum in rad/km.

# grid extent - grid of the British Isles in Irenet 95 grid projection
xmin = 0.0
xmax = 1900000.0
ymin = 0.0
ymax = 1700000.0
extent_grid = [xmin, xmax, ymin, ymax]

# extent on the sphere (WGS84)
extent_sphere = mapping.convert_extent(extent_grid, epsg_in=2157, epsg_out=4326)

# map extent - should be narrower than extent_sphere
extent_map = [-10.8, 2, 49.5, 59]

dx, dy = 1e3, 1e3 # 1 km grid resolution
nx, ny = int(round((xmax-xmin)/dx)), int(round((ymax-ymin)/dy))

# interpolate onto grid extent (for computation)
mag_grid = mapping.grid(mag_data[:,:2], mag_data[:,3], extent_grid, shape=(ny,nx), epsg_in=4326, epsg_out=2157)

# interpolate onto extent on the sphere (for mapping)
mag_sphere = mapping.grid(mag_data[:,:2], mag_data[:,3], extent_sphere, shape=(ny,nx))

proj = ccrs.UTM(30)

fig = plt.figure(figsize=(12,10))
ax = plt.axes(projection=proj)
ax.set_extent(extent_map)
ax.coastlines(resolution='50m', linewidth=1.5)
ax.gridlines()

im1 = ax.imshow(mag_sphere, extent=extent_sphere, transform=ccrs.PlateCarree(),
                cmap='Spectral_r', vmin=-200, vmax=200, zorder=0)

fig.colorbar(im1, label='nT')

Compute Curie depth

For this example we use the centroid method of Tanaka et al. (1999), as outlined in the previous notebooks. Both bands are given in rad/km, and check_bands is worth running against a guess at the source thickness before relying on them.

# fitting bands in rad/km -- see Ex2 for how to check these with
# grid.check_bands() before relying on them
zt_range = (1.257, 1.885)
z0_range = (0.0, 0.628)

grid = pycurious.CurieOptimiseTanaka(mag_grid, xmin, xmax, ymin, ymax)

window_size = 400e3

# centroid spacing of 50 km
xc_list, yc_list = grid.create_centroid_list(window_size, spacingX=50e3, spacingY=50e3)
print("number of centroids = {}".format(len(xc_list)))


zt, z0, zt_int, z0_int, zt_stdev, z0_stdev = grid.optimise_routine(window_size, xc_list, yc_list,
                                                                       zt_range, z0_range, taper=np.hamming)

CPD, sigma_CPD = grid.calculate_CPD(zt, z0, zt_stdev, z0_stdev)

centroids = np.column_stack([xc_list, yc_list])

CPD_sphere       = mapping.grid(centroids, CPD, extent_sphere, (ny,nx), epsg_in=2157, epsg_out=4326)
sigma_CPD_sphere = mapping.grid(centroids, sigma_CPD, extent_sphere, (ny,nx), epsg_in=2157, epsg_out=4326)

proj = ccrs.UTM(30)

fig = plt.figure(figsize=(12,10))
ax = plt.axes(projection=proj)
ax.set_extent(extent_map)
ax.coastlines(resolution='50m', linewidth=1.5)
ax.gridlines()

im1 = ax.imshow(CPD_sphere, extent=extent_sphere, transform=ccrs.PlateCarree(),
                cmap='BrBG',  zorder=0)

fig.colorbar(im1, label='km')

Compare to reference models

A global Curie depth reference model has been created by Li et al., 2017, which is a useful resource to compare these results. Download the supplementary material from the open access article to the same directory as EMAG2.

li2017_cpd_file = file_resources[1]['local_file']

li_cpd = np.loadtxt(li2017_cpd_file, skiprows=1)
li_cpd = li_cpd[~np.isnan(li_cpd[:,2])]

li_grid = mapping.grid(li_cpd[:,:2], li_cpd[:,2], extent_sphere, shape=(ny,nx))

proj = ccrs.UTM(30)

fig = plt.figure(figsize=(12,10))
ax = plt.axes(projection=proj)
ax.set_extent(extent_map)
ax.coastlines(resolution='50m', linewidth=1.5)
ax.gridlines()

im1 = ax.imshow(li_grid, extent=extent_sphere, transform=ccrs.PlateCarree(),
                cmap='BrBG', zorder=0, )

fig.colorbar(im1, label='km')

Export GeoTIFF

The GeoTIFF format contains projection information, which makes it ideal for GIS applications such as QGIS.

mapping.export_geotiff("my-curie-depth.tiff", CPD_sphere, extent_sphere, epsg=4326)

We can also import a GeoTIFF using the mapping module. the import_geotiff function prints a wealth of information on the reference spheriod and the projection, including EPSG code.

GTiff, GTextent = mapping.import_geotiff("my-curie-depth.tiff")

print(GTextent)
print(GTiff.shape)