<jupyter_start><jupyter_text>MIT LicenseCopyright (c) 2019 <NAME>, https://orcid.org/0000-0001-9626-8615 (ORCID)<jupyter_code>import xarray as xr
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Define functions<jupyter_code>from scipy.ndimage.filters import gaussian_filter

# band filter
def raster_gamma_range(raster0, g1, g2):
    raster = raster0.copy()
    raster.values = raster.values.astype(np.float32)
    raster.values = gaussian_filter(raster.values,g1) \
            - gaussian_filter(raster.values,g2)
    return raster
# dem, srtm
def correlogram(raster1, raster2, gammas):
    # spatial filtering
    rasters1 = []
    rasters2 = []
    for g in gammas:
        print (g,". ", end = '')

        _raster1 = raster_gamma_range(raster1, g-.5, g+.5)
        rasters1.append(_raster1)

        _raster2 = raster_gamma_range(raster2, g-.5, g+.5)
        rasters2.append(_raster2)
    print ()

    corrs = []
    for ridx in range(len(gammas)):
        print (ridx+1,". ", end = '')
        _raster2 = rasters2[ridx]
        for didx in range(len(gammas)):
            _raster1 = rasters1[didx]
            df = pd.DataFrame({'raster1': _raster1.values.flatten(), 'raster2': _raster2.values.flatten()})
            corr = round((df.corr()).iloc[0,1],2)
            corrs.append(corr)

    da_corr = xr.DataArray(np.array(corrs).reshape([len(gammas),len(gammas)]),
                          coords=[resolution*gammas,resolution*gammas],
                          dims=['raster2','raster1'])

    return (rasters1, rasters2, da_corr)<jupyter_output><empty_output><jupyter_text>Define parameters<jupyter_code># to load source data
GEBCO="GEBCO_2019.Africa30x30.tif"
GRAVITY="WGM2012_Bouguer_ponc_2min.Africa30x30.tif"

# rasters below defined in decimal degrees
# this coefficient [km/pixel] for pixel-based filtering and plotting
resolution = 3.7
GAMMA = 28
DGAMMA= 2<jupyter_output><empty_output><jupyter_text>Load datasets<jupyter_code>dem = xr.open_rasterio(GEBCO).rename({'x':'lon','y':'lat'})
dem
grav = xr.open_rasterio(GRAVITY).rename({'x':'lon','y':'lat'})
grav<jupyter_output><empty_output><jupyter_text>Compare source datasets<jupyter_code>fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10.5,5))

grav.plot(ax=ax1, cmap='terrain')
ax1.set_title('WGM2012 Bouguer',fontsize=16)

dem.plot(ax=ax2, cmap='terrain')
ax2.set_title('GEBCO_2019',fontsize=16)

fig.tight_layout(rect=[0.03, 0.0, 1, 0.9])
plt.suptitle('GEBCO_2019 and WGM2012 Bouguer Gravity Anomalies',fontsize=16)
#plt.savefig('GEBCO_2019 and WGM2012 Bouguer Gravity Anomalies.jpg', dpi=150)
plt.show()<jupyter_output><empty_output><jupyter_text>Make correlogram<jupyter_code>gammas = np.arange(1,GAMMA+DGAMMA/2,DGAMMA)
(dems,gravs,da_corr) = correlogram(dem, grav, gammas)
float(da_corr.min()),float(da_corr.max())
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10.5,5))

da_corr.plot(cmap='RdBu_r',ax=ax1, vmin=-1,vmax=1)
ax1.set_xlabel('GEBCO_2019 Wavelength, km',fontsize=12)
ax1.set_ylabel('WGM2012 Gravity Wavelength, km',fontsize=12)

da_corr.plot.contour(levels=np.linspace(-1,1,41),cmap='RdBu_r',add_colorbar=True, ax=ax2)
ax2.set_xlabel('GEBCO_2019 Wavelength, km',fontsize=12)
ax2.set_ylabel('WGM2012 Gravity Wavelength, km',fontsize=12)

plt.suptitle('Pearson Correlation Coefficient:\nGEBCO_2019 and WGM2012 Bouguer Gravity Anomalies',fontsize=16)
fig.tight_layout(rect=[0.03, 0.0, 1, 0.9])
#plt.savefig('Correlation GEBCO_2019 and WGM2012 Bouguer Gravity Anomalies.jpg', dpi=150)
plt.show()
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10.5,5))

# define wavelength [km] and index for it
wavelength = 80
gidx = np.argmin((gammas*resolution-wavelength)**2)

gravs[gidx].plot(cmap='RdBu',ax=ax1)
ax1.set_title('WGM2012 Bouguer',fontsize=16)

dems[gidx].plot(cmap='RdBu_r',ax=ax2)
ax2.set_title('GEBCO_2019',fontsize=16)

fig.tight_layout(rect=[0.03, 0.0, 1, 0.9])
plt.suptitle('Wavelength %dkm:\nGEBCO_2019 and WGM2012 Bouguer Gravity Anomalies\n' % wavelength,fontsize=16)
plt.show()
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(10.5,5))

# define wavelength [km] and index for it
wavelength = 100
gidx = np.argmin((gammas*resolution-wavelength)**2)

gravs[gidx].plot(cmap='RdBu',ax=ax1)
ax1.set_title('WGM2012 Bouguer',fontsize=16)

dems[gidx].plot(cmap='RdBu_r',ax=ax2)
ax2.set_title('GEBCO_2019',fontsize=16)

fig.tight_layout(rect=[0.03, 0.0, 1, 0.9])
plt.suptitle('Wavelength %dkm:\nGEBCO_2019 and WGM2012 Bouguer Gravity Anomalies\n' % wavelength,fontsize=16)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1554290718038529, 'fraction_numerical': 0.06589316987740805, 'mean_word_length': 5.536480686695279, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 26, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '199489', 'license': 'MIT', 'n_tokens_mistral': 2099, 'n_tokens_neox': 1895, 'n_words': 278, 'path': 'Africa/GEBCO_Bouguer_correlogram.Africa30x30.ipynb', 'repo_name': 'RichardScottOZ/gis-snippets', 'size': '4578'}