08 Geostatistics demo

UW Geospatial Data Analysis
CEE467/CEWA567
David Shean, Eric Gagliano, Quinn Brencher

Overview

This notebook explores fundamental geostatistical concepts using the GLAS elevation data from lab02

We’ll explore key geostatistical techniques:

1. Spatial autocorrelation (Global and Local Moran’s I)

2. Variogram analysis

3. Spatial interpolation (deterministic vs. kriging)

Imports

# !pip install pysal
# !pip install scikit-gstat

# Core libraries
import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
from matplotlib import colors
import os
from pathlib import Path
import warnings

# Geospatial libraries
import geopandas as gpd
import xarray as xr
import rasterio
import rioxarray as rxr
from pyproj import CRS
import contextily as ctx
from shapely.geometry import Point, Polygon

# Geostatistical libraries
import pysal
from pysal.explore import esda
from pysal.lib import weights
from splot.esda import moran_scatterplot
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
import skgstat as skg

import easysnowdata

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/spaghetti/network.py:41: FutureWarning: The next major release of pysal/spaghetti (2.0.0) will drop support for all ``libpysal.cg`` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ``libpysal.cg`` geometries should prepare for this simply by converting to ``shapely`` geometries.
  warnings.warn(dep_msg, FutureWarning, stacklevel=1)

Geostatistics with vector data

WUS_aea_crs = "+proj=aea +lon_0=-111.2695312 +lat_1=33.5538453 +lat_2=46.581155 +lat_0=40.0675002 +datum=WGS84 +units=m +no_defs"
glas_df = pd.read_csv('../02_NumPy_Pandas_Matplotlib/data/GLAH14_tllz_conus_lulcfilt_demfilt.csv')
glas_gdf = gpd.GeoDataFrame(glas_df, geometry=gpd.points_from_xy(glas_df.lon, glas_df.lat)).set_crs('EPSG:4326')
glas_gdf = glas_gdf.to_crs(WUS_aea_crs)
glas_gdf

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry
0	2003.139571	731266.943345	44.157897	-105.356562	1398.51	1400.52	0.33	31	POINT (470693.831 472504.893)
1	2003.139571	731266.943346	44.150175	-105.358116	1387.11	1384.64	0.43	31	POINT (470627.12 471636.996)
2	2003.139571	731266.943347	44.148632	-105.358427	1392.83	1383.49	0.28	31	POINT (470613.749 471463.57)
3	2003.139571	731266.943347	44.147087	-105.358738	1384.24	1382.85	0.84	31	POINT (470600.39 471289.921)
4	2003.139571	731266.943347	44.145542	-105.359048	1369.21	1380.24	1.73	31	POINT (470587.11 471116.278)
...	...	...	...	...	...	...	...	...	...
65231	2009.775995	733691.238340	37.896222	-117.044399	1556.16	1556.43	0.00	31	POINT (-504775.82 -226245.044)
65232	2009.775995	733691.238340	37.897769	-117.044675	1556.02	1556.43	0.00	31	POINT (-504788.788 -226071.184)
65233	2009.775995	733691.238340	37.899319	-117.044952	1556.19	1556.44	0.00	31	POINT (-504801.821 -225896.984)
65234	2009.775995	733691.238340	37.900869	-117.045230	1556.18	1556.44	0.00	31	POINT (-504814.94 -225722.778)
65235	2009.775995	733691.238341	37.902420	-117.045508	1556.32	1556.44	0.00	31	POINT (-504828.05 -225548.46)

65236 rows × 9 columns

f,ax=plt.subplots(figsize=(10,8))

glas_gdf.plot(ax=ax,column='glas_z',cmap='inferno',legend=True,markersize=1,legend_kwds={'label': "Elevation (m)"})
ctx.add_basemap(ax,crs=glas_gdf.crs.to_string())

ax.set_aspect('equal')

ax.set_title('GLAS Elevation')

f.tight_layout()

Regular correlation - comparing glas_z and dem_z values

corr = glas_gdf.corr(numeric_only=True)
corr

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc
decyear	1.000000	1.000000	-0.030542	0.022455	-0.010960	-0.011376	-0.023535	-0.006947
ordinal	1.000000	1.000000	-0.030537	0.022457	-0.010963	-0.011379	-0.023538	-0.006953
lat	-0.030542	-0.030537	1.000000	-0.159387	-0.112952	-0.113471	0.103647	-0.111976
lon	0.022455	0.022457	-0.159387	1.000000	0.507210	0.507571	0.141909	-0.073070
glas_z	-0.010960	-0.010963	-0.112952	0.507210	1.000000	0.999929	0.588513	-0.206267
dem_z	-0.011376	-0.011379	-0.113471	0.507571	0.999929	1.000000	0.589164	-0.206353
dem_z_std	-0.023535	-0.023538	0.103647	0.141909	0.588513	0.589164	1.000000	-0.106402
lulc	-0.006947	-0.006953	-0.111976	-0.073070	-0.206267	-0.206353	-0.106402	1.000000

m,b = np.polyfit(glas_gdf['glas_z'], glas_gdf['dem_z'], 1)

f,ax=plt.subplots(figsize=(10,10))
glas_gdf.plot(x='glas_z', y='dem_z', kind='scatter', ax=ax)
ax.set_aspect('equal')
ax.plot([0, 4500], [0, 4500], color='black', linestyle='--')
ax.set_title(f'glas_z vs dem_z\nbest fit line: y={m:.2f}+{b:.2f}x\npearson r={corr.loc["glas_z","dem_z"]:.5f}')
ax.plot(glas_gdf['glas_z'], m*glas_gdf['glas_z']+b, color='red',linewidth=0.5)

[<matplotlib.lines.Line2D at 0x7f30846fed20>]

Spatial autocorrelation - Global Moran’s I

First, we’ll take a sample of our points and create a spatial weights matrix. Experiment with some of the spatial weight / connectivity definitions!

sample_size = min(10000, len(glas_gdf))  # Adjust based on computational resources
glas_sample = glas_gdf.sample(sample_size, random_state=42)

# Create a spatial weights matrix
# We'll use k-nearest neighbors (KNN) weights
k = 8  # Number of neighbors
knn_w = weights.KNN.from_dataframe(glas_sample, k=k)
knn_w = weights.distance.DistanceBand.from_dataframe(glas_sample, threshold=5000)
knn_w.transform = 'r'  

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/libpysal/weights/distance.py:153: UserWarning: The weights matrix is not fully connected: 
 There are 70 disconnected components.
  W.__init__(self, neighbors, id_order=ids, **kwargs)

('WARNING: ', 14865, ' is an island (no neighbors)')
('WARNING: ', 8557, ' is an island (no neighbors)')
('WARNING: ', 6769, ' is an island (no neighbors)')
('WARNING: ', 15460, ' is an island (no neighbors)')
('WARNING: ', 24589, ' is an island (no neighbors)')
('WARNING: ', 24581, ' is an island (no neighbors)')
('WARNING: ', 18310, ' is an island (no neighbors)')
('WARNING: ', 5382, ' is an island (no neighbors)')
('WARNING: ', 45376, ' is an island (no neighbors)')
('WARNING: ', 21821, ' is an island (no neighbors)')
('WARNING: ', 6288, ' is an island (no neighbors)')
('WARNING: ', 53823, ' is an island (no neighbors)')
('WARNING: ', 56887, ' is an island (no neighbors)')
('WARNING: ', 3308, ' is an island (no neighbors)')
('WARNING: ', 39216, ' is an island (no neighbors)')
('WARNING: ', 26231, ' is an island (no neighbors)')
('WARNING: ', 60581, ' is an island (no neighbors)')
('WARNING: ', 34683, ' is an island (no neighbors)')
('WARNING: ', 22120, ' is an island (no neighbors)')
('WARNING: ', 23900, ' is an island (no neighbors)')
('WARNING: ', 21844, ' is an island (no neighbors)')
('WARNING: ', 40189, ' is an island (no neighbors)')
('WARNING: ', 2791, ' is an island (no neighbors)')
('WARNING: ', 4448, ' is an island (no neighbors)')
('WARNING: ', 21588, ' is an island (no neighbors)')
('WARNING: ', 62377, ' is an island (no neighbors)')
('WARNING: ', 34004, ' is an island (no neighbors)')
('WARNING: ', 24983, ' is an island (no neighbors)')
('WARNING: ', 53791, ' is an island (no neighbors)')
('WARNING: ', 55983, ' is an island (no neighbors)')
('WARNING: ', 12892, ' is an island (no neighbors)')
('WARNING: ', 55989, ' is an island (no neighbors)')
('WARNING: ', 59932, ' is an island (no neighbors)')
('WARNING: ', 6234, ' is an island (no neighbors)')
('WARNING: ', 2527, ' is an island (no neighbors)')
('WARNING: ', 15318, ' is an island (no neighbors)')
('WARNING: ', 8492, ' is an island (no neighbors)')
('WARNING: ', 38604, ' is an island (no neighbors)')
('WARNING: ', 7107, ' is an island (no neighbors)')
('WARNING: ', 5517, ' is an island (no neighbors)')
('WARNING: ', 41014, ' is an island (no neighbors)')
('WARNING: ', 16140, ' is an island (no neighbors)')
('WARNING: ', 14624, ' is an island (no neighbors)')
('WARNING: ', 11310, ' is an island (no neighbors)')
('WARNING: ', 27248, ' is an island (no neighbors)')
('WARNING: ', 36391, ' is an island (no neighbors)')
('WARNING: ', 45846, ' is an island (no neighbors)')
('WARNING: ', 37112, ' is an island (no neighbors)')
('WARNING: ', 21808, ' is an island (no neighbors)')
('WARNING: ', 31079, ' is an island (no neighbors)')
('WARNING: ', 1736, ' is an island (no neighbors)')
('WARNING: ', 48743, ' is an island (no neighbors)')
('WARNING: ', 2541, ' is an island (no neighbors)')
('WARNING: ', 15969, ' is an island (no neighbors)')
('WARNING: ', 47256, ' is an island (no neighbors)')
('WARNING: ', 5540, ' is an island (no neighbors)')
('WARNING: ', 38807, ' is an island (no neighbors)')
('WARNING: ', 33648, ' is an island (no neighbors)')
('WARNING: ', 36209, ' is an island (no neighbors)')
('WARNING: ', 5476, ' is an island (no neighbors)')
('WARNING: ', 2182, ' is an island (no neighbors)')
('WARNING: ', 18307, ' is an island (no neighbors)')
('WARNING: ', 61370, ' is an island (no neighbors)')
('WARNING: ', 37151, ' is an island (no neighbors)')
('WARNING: ', 12900, ' is an island (no neighbors)')
('WARNING: ', 60893, ' is an island (no neighbors)')
('WARNING: ', 9670, ' is an island (no neighbors)')
('WARNING: ', 30717, ' is an island (no neighbors)')
('WARNING: ', 21018, ' is an island (no neighbors)')
('WARNING: ', 54196, ' is an island (no neighbors)')
('WARNING: ', 21838, ' is an island (no neighbors)')
('WARNING: ', 31522, ' is an island (no neighbors)')
('WARNING: ', 55603, ' is an island (no neighbors)')
('WARNING: ', 13010, ' is an island (no neighbors)')
('WARNING: ', 39485, ' is an island (no neighbors)')
('WARNING: ', 9093, ' is an island (no neighbors)')
('WARNING: ', 26791, ' is an island (no neighbors)')
('WARNING: ', 15302, ' is an island (no neighbors)')
('WARNING: ', 5547, ' is an island (no neighbors)')
('WARNING: ', 55728, ' is an island (no neighbors)')
('WARNING: ', 53157, ' is an island (no neighbors)')
('WARNING: ', 17599, ' is an island (no neighbors)')
('WARNING: ', 19321, ' is an island (no neighbors)')
('WARNING: ', 51062, ' is an island (no neighbors)')
('WARNING: ', 33652, ' is an island (no neighbors)')
('WARNING: ', 13835, ' is an island (no neighbors)')
('WARNING: ', 63592, ' is an island (no neighbors)')
('WARNING: ', 13459, ' is an island (no neighbors)')
('WARNING: ', 41015, ' is an island (no neighbors)')
('WARNING: ', 42118, ' is an island (no neighbors)')
('WARNING: ', 33636, ' is an island (no neighbors)')
('WARNING: ', 8375, ' is an island (no neighbors)')
('WARNING: ', 55941, ' is an island (no neighbors)')
('WARNING: ', 6268, ' is an island (no neighbors)')
('WARNING: ', 8925, ' is an island (no neighbors)')
('WARNING: ', 21841, ' is an island (no neighbors)')
('WARNING: ', 47813, ' is an island (no neighbors)')
('WARNING: ', 18441, ' is an island (no neighbors)')
('WARNING: ', 32163, ' is an island (no neighbors)')
('WARNING: ', 1805, ' is an island (no neighbors)')
('WARNING: ', 56839, ' is an island (no neighbors)')
('WARNING: ', 61971, ' is an island (no neighbors)')
('WARNING: ', 60492, ' is an island (no neighbors)')
('WARNING: ', 5528, ' is an island (no neighbors)')
('WARNING: ', 56855, ' is an island (no neighbors)')
('WARNING: ', 14898, ' is an island (no neighbors)')
('WARNING: ', 22130, ' is an island (no neighbors)')
('WARNING: ', 32159, ' is an island (no neighbors)')
('WARNING: ', 3452, ' is an island (no neighbors)')
('WARNING: ', 64506, ' is an island (no neighbors)')
('WARNING: ', 3431, ' is an island (no neighbors)')
('WARNING: ', 61338, ' is an island (no neighbors)')
('WARNING: ', 40568, ' is an island (no neighbors)')
('WARNING: ', 55725, ' is an island (no neighbors)')
('WARNING: ', 45990, ' is an island (no neighbors)')
('WARNING: ', 21597, ' is an island (no neighbors)')
('WARNING: ', 23945, ' is an island (no neighbors)')
('WARNING: ', 62597, ' is an island (no neighbors)')
('WARNING: ', 54184, ' is an island (no neighbors)')
('WARNING: ', 49169, ' is an island (no neighbors)')
('WARNING: ', 6235, ' is an island (no neighbors)')
('WARNING: ', 16339, ' is an island (no neighbors)')
('WARNING: ', 47189, ' is an island (no neighbors)')
('WARNING: ', 13941, ' is an island (no neighbors)')
('WARNING: ', 3882, ' is an island (no neighbors)')
('WARNING: ', 6293, ' is an island (no neighbors)')
('WARNING: ', 27888, ' is an island (no neighbors)')
('WARNING: ', 53042, ' is an island (no neighbors)')
('WARNING: ', 1057, ' is an island (no neighbors)')
('WARNING: ', 3940, ' is an island (no neighbors)')
('WARNING: ', 13905, ' is an island (no neighbors)')
('WARNING: ', 23822, ' is an island (no neighbors)')
('WARNING: ', 5718, ' is an island (no neighbors)')
('WARNING: ', 56868, ' is an island (no neighbors)')
('WARNING: ', 32706, ' is an island (no neighbors)')
('WARNING: ', 14448, ' is an island (no neighbors)')
('WARNING: ', 55945, ' is an island (no neighbors)')
('WARNING: ', 62576, ' is an island (no neighbors)')
('WARNING: ', 44861, ' is an island (no neighbors)')
('WARNING: ', 21854, ' is an island (no neighbors)')
('WARNING: ', 23113, ' is an island (no neighbors)')
('WARNING: ', 2592, ' is an island (no neighbors)')
('WARNING: ', 49176, ' is an island (no neighbors)')
('WARNING: ', 32660, ' is an island (no neighbors)')
('WARNING: ', 9203, ' is an island (no neighbors)')
('WARNING: ', 53796, ' is an island (no neighbors)')
('WARNING: ', 64534, ' is an island (no neighbors)')
('WARNING: ', 29115, ' is an island (no neighbors)')
('WARNING: ', 2528, ' is an island (no neighbors)')
('WARNING: ', 24598, ' is an island (no neighbors)')
('WARNING: ', 63874, ' is an island (no neighbors)')
('WARNING: ', 32042, ' is an island (no neighbors)')
('WARNING: ', 52069, ' is an island (no neighbors)')
('WARNING: ', 23801, ' is an island (no neighbors)')
('WARNING: ', 64731, ' is an island (no neighbors)')
('WARNING: ', 86, ' is an island (no neighbors)')
('WARNING: ', 14261, ' is an island (no neighbors)')
('WARNING: ', 6860, ' is an island (no neighbors)')
('WARNING: ', 22952, ' is an island (no neighbors)')
('WARNING: ', 24468, ' is an island (no neighbors)')
('WARNING: ', 50437, ' is an island (no neighbors)')
('WARNING: ', 63320, ' is an island (no neighbors)')
('WARNING: ', 41634, ' is an island (no neighbors)')
('WARNING: ', 60750, ' is an island (no neighbors)')
('WARNING: ', 61119, ' is an island (no neighbors)')
('WARNING: ', 46676, ' is an island (no neighbors)')
('WARNING: ', 36421, ' is an island (no neighbors)')
('WARNING: ', 23825, ' is an island (no neighbors)')
('WARNING: ', 41172, ' is an island (no neighbors)')
('WARNING: ', 18303, ' is an island (no neighbors)')
('WARNING: ', 13912, ' is an island (no neighbors)')
('WARNING: ', 4172, ' is an island (no neighbors)')
('WARNING: ', 24961, ' is an island (no neighbors)')
('WARNING: ', 23214, ' is an island (no neighbors)')
('WARNING: ', 39387, ' is an island (no neighbors)')
('WARNING: ', 32648, ' is an island (no neighbors)')
('WARNING: ', 60358, ' is an island (no neighbors)')
('WARNING: ', 7390, ' is an island (no neighbors)')
('WARNING: ', 11307, ' is an island (no neighbors)')
('WARNING: ', 5428, ' is an island (no neighbors)')
('WARNING: ', 39900, ' is an island (no neighbors)')
('WARNING: ', 56365, ' is an island (no neighbors)')
('WARNING: ', 27436, ' is an island (no neighbors)')
('WARNING: ', 41391, ' is an island (no neighbors)')
('WARNING: ', 5672, ' is an island (no neighbors)')
('WARNING: ', 41177, ' is an island (no neighbors)')
('WARNING: ', 9594, ' is an island (no neighbors)')
('WARNING: ', 5732, ' is an island (no neighbors)')
('WARNING: ', 45602, ' is an island (no neighbors)')
('WARNING: ', 11247, ' is an island (no neighbors)')
('WARNING: ', 48575, ' is an island (no neighbors)')
('WARNING: ', 6793, ' is an island (no neighbors)')
('WARNING: ', 44374, ' is an island (no neighbors)')
('WARNING: ', 34127, ' is an island (no neighbors)')
('WARNING: ', 3222, ' is an island (no neighbors)')
('WARNING: ', 13757, ' is an island (no neighbors)')
('WARNING: ', 33136, ' is an island (no neighbors)')
('WARNING: ', 42374, ' is an island (no neighbors)')
('WARNING: ', 9279, ' is an island (no neighbors)')
('WARNING: ', 5598, ' is an island (no neighbors)')
('WARNING: ', 54538, ' is an island (no neighbors)')
('WARNING: ', 12353, ' is an island (no neighbors)')
('WARNING: ', 24114, ' is an island (no neighbors)')
('WARNING: ', 55732, ' is an island (no neighbors)')
('WARNING: ', 49963, ' is an island (no neighbors)')
('WARNING: ', 6280, ' is an island (no neighbors)')
('WARNING: ', 10950, ' is an island (no neighbors)')
('WARNING: ', 54892, ' is an island (no neighbors)')
('WARNING: ', 6724, ' is an island (no neighbors)')
('WARNING: ', 41647, ' is an island (no neighbors)')
('WARNING: ', 48736, ' is an island (no neighbors)')
('WARNING: ', 44469, ' is an island (no neighbors)')
('WARNING: ', 6756, ' is an island (no neighbors)')
('WARNING: ', 45762, ' is an island (no neighbors)')
('WARNING: ', 34138, ' is an island (no neighbors)')
('WARNING: ', 4475, ' is an island (no neighbors)')
('WARNING: ', 28846, ' is an island (no neighbors)')
('WARNING: ', 335, ' is an island (no neighbors)')
('WARNING: ', 800, ' is an island (no neighbors)')
('WARNING: ', 57483, ' is an island (no neighbors)')
('WARNING: ', 33140, ' is an island (no neighbors)')
('WARNING: ', 6839, ' is an island (no neighbors)')
('WARNING: ', 15296, ' is an island (no neighbors)')
('WARNING: ', 17523, ' is an island (no neighbors)')
('WARNING: ', 42265, ' is an island (no neighbors)')
('WARNING: ', 6247, ' is an island (no neighbors)')
('WARNING: ', 60766, ' is an island (no neighbors)')
('WARNING: ', 8972, ' is an island (no neighbors)')
('WARNING: ', 49665, ' is an island (no neighbors)')
('WARNING: ', 24115, ' is an island (no neighbors)')
('WARNING: ', 27966, ' is an island (no neighbors)')
('WARNING: ', 63882, ' is an island (no neighbors)')
('WARNING: ', 56765, ' is an island (no neighbors)')
('WARNING: ', 21843, ' is an island (no neighbors)')
('WARNING: ', 7106, ' is an island (no neighbors)')
('WARNING: ', 21879, ' is an island (no neighbors)')
('WARNING: ', 12139, ' is an island (no neighbors)')
('WARNING: ', 50898, ' is an island (no neighbors)')
('WARNING: ', 21349, ' is an island (no neighbors)')
('WARNING: ', 8550, ' is an island (no neighbors)')
('WARNING: ', 5572, ' is an island (no neighbors)')
('WARNING: ', 22125, ' is an island (no neighbors)')
('WARNING: ', 58079, ' is an island (no neighbors)')
('WARNING: ', 49959, ' is an island (no neighbors)')
('WARNING: ', 1043, ' is an island (no neighbors)')
('WARNING: ', 9207, ' is an island (no neighbors)')
('WARNING: ', 43149, ' is an island (no neighbors)')
('WARNING: ', 63105, ' is an island (no neighbors)')
('WARNING: ', 39118, ' is an island (no neighbors)')
('WARNING: ', 30824, ' is an island (no neighbors)')
('WARNING: ', 64539, ' is an island (no neighbors)')
('WARNING: ', 2579, ' is an island (no neighbors)')
('WARNING: ', 27239, ' is an island (no neighbors)')
('WARNING: ', 59942, ' is an island (no neighbors)')
('WARNING: ', 22133, ' is an island (no neighbors)')
('WARNING: ', 45630, ' is an island (no neighbors)')
('WARNING: ', 57335, ' is an island (no neighbors)')
('WARNING: ', 62564, ' is an island (no neighbors)')
('WARNING: ', 55604, ' is an island (no neighbors)')
('WARNING: ', 33720, ' is an island (no neighbors)')
('WARNING: ', 1059, ' is an island (no neighbors)')
('WARNING: ', 55897, ' is an island (no neighbors)')
('WARNING: ', 31082, ' is an island (no neighbors)')
('WARNING: ', 39903, ' is an island (no neighbors)')
('WARNING: ', 56425, ' is an island (no neighbors)')
('WARNING: ', 13726, ' is an island (no neighbors)')
('WARNING: ', 32623, ' is an island (no neighbors)')
('WARNING: ', 60784, ' is an island (no neighbors)')
('WARNING: ', 36371, ' is an island (no neighbors)')
('WARNING: ', 46618, ' is an island (no neighbors)')
('WARNING: ', 44166, ' is an island (no neighbors)')
('WARNING: ', 41729, ' is an island (no neighbors)')
('WARNING: ', 44611, ' is an island (no neighbors)')
('WARNING: ', 55040, ' is an island (no neighbors)')
('WARNING: ', 55607, ' is an island (no neighbors)')
('WARNING: ', 52280, ' is an island (no neighbors)')
('WARNING: ', 11730, ' is an island (no neighbors)')
('WARNING: ', 1033, ' is an island (no neighbors)')
('WARNING: ', 462, ' is an island (no neighbors)')
('WARNING: ', 19988, ' is an island (no neighbors)')
('WARNING: ', 15364, ' is an island (no neighbors)')
('WARNING: ', 6821, ' is an island (no neighbors)')
('WARNING: ', 4935, ' is an island (no neighbors)')
('WARNING: ', 44034, ' is an island (no neighbors)')
('WARNING: ', 26517, ' is an island (no neighbors)')
('WARNING: ', 21726, ' is an island (no neighbors)')
('WARNING: ', 22395, ' is an island (no neighbors)')
('WARNING: ', 16174, ' is an island (no neighbors)')
('WARNING: ', 62362, ' is an island (no neighbors)')
('WARNING: ', 61968, ' is an island (no neighbors)')
('WARNING: ', 10852, ' is an island (no neighbors)')
('WARNING: ', 19366, ' is an island (no neighbors)')
('WARNING: ', 9069, ' is an island (no neighbors)')
('WARNING: ', 34463, ' is an island (no neighbors)')
('WARNING: ', 19922, ' is an island (no neighbors)')
('WARNING: ', 12362, ' is an island (no neighbors)')
('WARNING: ', 35235, ' is an island (no neighbors)')
('WARNING: ', 50433, ' is an island (no neighbors)')
('WARNING: ', 29874, ' is an island (no neighbors)')
('WARNING: ', 12858, ' is an island (no neighbors)')
('WARNING: ', 828, ' is an island (no neighbors)')
('WARNING: ', 3732, ' is an island (no neighbors)')
('WARNING: ', 63785, ' is an island (no neighbors)')
('WARNING: ', 5954, ' is an island (no neighbors)')
('WARNING: ', 12081, ' is an island (no neighbors)')
('WARNING: ', 5927, ' is an island (no neighbors)')
('WARNING: ', 53380, ' is an island (no neighbors)')
('WARNING: ', 60323, ' is an island (no neighbors)')
('WARNING: ', 57256, ' is an island (no neighbors)')
('WARNING: ', 60783, ' is an island (no neighbors)')
('WARNING: ', 27008, ' is an island (no neighbors)')
('WARNING: ', 13904, ' is an island (no neighbors)')
('WARNING: ', 5537, ' is an island (no neighbors)')
('WARNING: ', 62376, ' is an island (no neighbors)')
('WARNING: ', 44175, ' is an island (no neighbors)')
('WARNING: ', 22108, ' is an island (no neighbors)')
('WARNING: ', 15443, ' is an island (no neighbors)')
('WARNING: ', 51061, ' is an island (no neighbors)')
('WARNING: ', 48913, ' is an island (no neighbors)')
('WARNING: ', 43346, ' is an island (no neighbors)')
('WARNING: ', 41602, ' is an island (no neighbors)')
('WARNING: ', 49836, ' is an island (no neighbors)')
('WARNING: ', 14286, ' is an island (no neighbors)')
('WARNING: ', 9660, ' is an island (no neighbors)')
('WARNING: ', 44511, ' is an island (no neighbors)')
('WARNING: ', 19313, ' is an island (no neighbors)')
('WARNING: ', 60357, ' is an island (no neighbors)')
('WARNING: ', 33647, ' is an island (no neighbors)')
('WARNING: ', 28547, ' is an island (no neighbors)')
('WARNING: ', 29112, ' is an island (no neighbors)')
('WARNING: ', 38801, ' is an island (no neighbors)')
('WARNING: ', 3179, ' is an island (no neighbors)')
('WARNING: ', 2059, ' is an island (no neighbors)')
('WARNING: ', 22379, ' is an island (no neighbors)')
('WARNING: ', 49189, ' is an island (no neighbors)')
('WARNING: ', 6804, ' is an island (no neighbors)')
('WARNING: ', 57297, ' is an island (no neighbors)')
('WARNING: ', 59919, ' is an island (no neighbors)')
('WARNING: ', 6199, ' is an island (no neighbors)')
('WARNING: ', 60370, ' is an island (no neighbors)')
('WARNING: ', 12131, ' is an island (no neighbors)')
('WARNING: ', 10369, ' is an island (no neighbors)')
('WARNING: ', 4948, ' is an island (no neighbors)')
('WARNING: ', 25526, ' is an island (no neighbors)')
('WARNING: ', 26546, ' is an island (no neighbors)')
('WARNING: ', 15388, ' is an island (no neighbors)')
('WARNING: ', 4391, ' is an island (no neighbors)')
('WARNING: ', 49249, ' is an island (no neighbors)')
('WARNING: ', 51003, ' is an island (no neighbors)')
('WARNING: ', 5148, ' is an island (no neighbors)')
('WARNING: ', 28283, ' is an island (no neighbors)')
('WARNING: ', 6196, ' is an island (no neighbors)')
('WARNING: ', 42997, ' is an island (no neighbors)')
('WARNING: ', 44601, ' is an island (no neighbors)')
('WARNING: ', 10574, ' is an island (no neighbors)')
('WARNING: ', 49255, ' is an island (no neighbors)')
('WARNING: ', 35210, ' is an island (no neighbors)')
('WARNING: ', 586, ' is an island (no neighbors)')
('WARNING: ', 25657, ' is an island (no neighbors)')
('WARNING: ', 57770, ' is an island (no neighbors)')
('WARNING: ', 11306, ' is an island (no neighbors)')
('WARNING: ', 43142, ' is an island (no neighbors)')
('WARNING: ', 6079, ' is an island (no neighbors)')
('WARNING: ', 8350, ' is an island (no neighbors)')
('WARNING: ', 58096, ' is an island (no neighbors)')
('WARNING: ', 39899, ' is an island (no neighbors)')
('WARNING: ', 9011, ' is an island (no neighbors)')
('WARNING: ', 22617, ' is an island (no neighbors)')
('WARNING: ', 42333, ' is an island (no neighbors)')
('WARNING: ', 7162, ' is an island (no neighbors)')
('WARNING: ', 20318, ' is an island (no neighbors)')
('WARNING: ', 11518, ' is an island (no neighbors)')
('WARNING: ', 28286, ' is an island (no neighbors)')
('WARNING: ', 18848, ' is an island (no neighbors)')
('WARNING: ', 57240, ' is an island (no neighbors)')
('WARNING: ', 60018, ' is an island (no neighbors)')
('WARNING: ', 19842, ' is an island (no neighbors)')
('WARNING: ', 54974, ' is an island (no neighbors)')
('WARNING: ', 6165, ' is an island (no neighbors)')
('WARNING: ', 1574, ' is an island (no neighbors)')
('WARNING: ', 62364, ' is an island (no neighbors)')
('WARNING: ', 6315, ' is an island (no neighbors)')
('WARNING: ', 11342, ' is an island (no neighbors)')
('WARNING: ', 49389, ' is an island (no neighbors)')
('WARNING: ', 62357, ' is an island (no neighbors)')
('WARNING: ', 22815, ' is an island (no neighbors)')
('WARNING: ', 9654, ' is an island (no neighbors)')
('WARNING: ', 52978, ' is an island (no neighbors)')
('WARNING: ', 6166, ' is an island (no neighbors)')
('WARNING: ', 23346, ' is an island (no neighbors)')
('WARNING: ', 56165, ' is an island (no neighbors)')
('WARNING: ', 5683, ' is an island (no neighbors)')
('WARNING: ', 2020, ' is an island (no neighbors)')
('WARNING: ', 20993, ' is an island (no neighbors)')
('WARNING: ', 6522, ' is an island (no neighbors)')
('WARNING: ', 56888, ' is an island (no neighbors)')
('WARNING: ', 31380, ' is an island (no neighbors)')
('WARNING: ', 27893, ' is an island (no neighbors)')
('WARNING: ', 49019, ' is an island (no neighbors)')
('WARNING: ', 15679, ' is an island (no neighbors)')
('WARNING: ', 36995, ' is an island (no neighbors)')
('WARNING: ', 21805, ' is an island (no neighbors)')
('WARNING: ', 2181, ' is an island (no neighbors)')
('WARNING: ', 38545, ' is an island (no neighbors)')
('WARNING: ', 5554, ' is an island (no neighbors)')
('WARNING: ', 53001, ' is an island (no neighbors)')
('WARNING: ', 9661, ' is an island (no neighbors)')
('WARNING: ', 19438, ' is an island (no neighbors)')
('WARNING: ', 1572, ' is an island (no neighbors)')
('WARNING: ', 15974, ' is an island (no neighbors)')
('WARNING: ', 22124, ' is an island (no neighbors)')
('WARNING: ', 1576, ' is an island (no neighbors)')
('WARNING: ', 14672, ' is an island (no neighbors)')
('WARNING: ', 58098, ' is an island (no neighbors)')
('WARNING: ', 23919, ' is an island (no neighbors)')
('WARNING: ', 14940, ' is an island (no neighbors)')
('WARNING: ', 1039, ' is an island (no neighbors)')
('WARNING: ', 27969, ' is an island (no neighbors)')
('WARNING: ', 20065, ' is an island (no neighbors)')
('WARNING: ', 27026, ' is an island (no neighbors)')
('WARNING: ', 32641, ' is an island (no neighbors)')
('WARNING: ', 18442, ' is an island (no neighbors)')
('WARNING: ', 31454, ' is an island (no neighbors)')
('WARNING: ', 39868, ' is an island (no neighbors)')
('WARNING: ', 53165, ' is an island (no neighbors)')
('WARNING: ', 14243, ' is an island (no neighbors)')
('WARNING: ', 20003, ' is an island (no neighbors)')
('WARNING: ', 32165, ' is an island (no neighbors)')
('WARNING: ', 55898, ' is an island (no neighbors)')
('WARNING: ', 34124, ' is an island (no neighbors)')
('WARNING: ', 21847, ' is an island (no neighbors)')
('WARNING: ', 35669, ' is an island (no neighbors)')
('WARNING: ', 41799, ' is an island (no neighbors)')
('WARNING: ', 59783, ' is an island (no neighbors)')
('WARNING: ', 26891, ' is an island (no neighbors)')
('WARNING: ', 38771, ' is an island (no neighbors)')
('WARNING: ', 12136, ' is an island (no neighbors)')
('WARNING: ', 11065, ' is an island (no neighbors)')
('WARNING: ', 60871, ' is an island (no neighbors)')
('WARNING: ', 33486, ' is an island (no neighbors)')
('WARNING: ', 21001, ' is an island (no neighbors)')
('WARNING: ', 28321, ' is an island (no neighbors)')
('WARNING: ', 23815, ' is an island (no neighbors)')
('WARNING: ', 40255, ' is an island (no neighbors)')
('WARNING: ', 12358, ' is an island (no neighbors)')
('WARNING: ', 47186, ' is an island (no neighbors)')
('WARNING: ', 54979, ' is an island (no neighbors)')
('WARNING: ', 48337, ' is an island (no neighbors)')
('WARNING: ', 11744, ' is an island (no neighbors)')
('WARNING: ', 28907, ' is an island (no neighbors)')
('WARNING: ', 53985, ' is an island (no neighbors)')
('WARNING: ', 39716, ' is an island (no neighbors)')
('WARNING: ', 17516, ' is an island (no neighbors)')
('WARNING: ', 26232, ' is an island (no neighbors)')
('WARNING: ', 32146, ' is an island (no neighbors)')
('WARNING: ', 5739, ' is an island (no neighbors)')
('WARNING: ', 27694, ' is an island (no neighbors)')
('WARNING: ', 9008, ' is an island (no neighbors)')
('WARNING: ', 64542, ' is an island (no neighbors)')
('WARNING: ', 29343, ' is an island (no neighbors)')
('WARNING: ', 6717, ' is an island (no neighbors)')
('WARNING: ', 80, ' is an island (no neighbors)')
('WARNING: ', 10708, ' is an island (no neighbors)')
('WARNING: ', 41601, ' is an island (no neighbors)')
('WARNING: ', 28282, ' is an island (no neighbors)')
('WARNING: ', 33646, ' is an island (no neighbors)')
('WARNING: ', 19433, ' is an island (no neighbors)')
('WARNING: ', 24582, ' is an island (no neighbors)')
('WARNING: ', 49586, ' is an island (no neighbors)')
('WARNING: ', 4939, ' is an island (no neighbors)')
('WARNING: ', 31060, ' is an island (no neighbors)')
('WARNING: ', 64261, ' is an island (no neighbors)')
('WARNING: ', 27171, ' is an island (no neighbors)')
('WARNING: ', 37461, ' is an island (no neighbors)')
('WARNING: ', 5522, ' is an island (no neighbors)')
('WARNING: ', 7687, ' is an island (no neighbors)')
('WARNING: ', 60187, ' is an island (no neighbors)')
('WARNING: ', 5426, ' is an island (no neighbors)')
('WARNING: ', 54090, ' is an island (no neighbors)')
('WARNING: ', 10351, ' is an island (no neighbors)')
('WARNING: ', 27831, ' is an island (no neighbors)')
('WARNING: ', 31523, ' is an island (no neighbors)')
('WARNING: ', 12913, ' is an island (no neighbors)')
('WARNING: ', 32100, ' is an island (no neighbors)')
('WARNING: ', 60412, ' is an island (no neighbors)')
('WARNING: ', 45593, ' is an island (no neighbors)')
('WARNING: ', 39658, ' is an island (no neighbors)')
('WARNING: ', 21543, ' is an island (no neighbors)')
('WARNING: ', 58111, ' is an island (no neighbors)')
('WARNING: ', 16338, ' is an island (no neighbors)')
('WARNING: ', 17394, ' is an island (no neighbors)')
('WARNING: ', 38332, ' is an island (no neighbors)')
('WARNING: ', 2721, ' is an island (no neighbors)')
('WARNING: ', 32645, ' is an island (no neighbors)')
('WARNING: ', 35836, ' is an island (no neighbors)')
('WARNING: ', 4484, ' is an island (no neighbors)')
('WARNING: ', 61325, ' is an island (no neighbors)')
('WARNING: ', 38770, ' is an island (no neighbors)')
('WARNING: ', 49023, ' is an island (no neighbors)')
('WARNING: ', 13935, ' is an island (no neighbors)')
('WARNING: ', 37373, ' is an island (no neighbors)')
('WARNING: ', 57975, ' is an island (no neighbors)')
('WARNING: ', 16330, ' is an island (no neighbors)')
('WARNING: ', 29867, ' is an island (no neighbors)')
('WARNING: ', 23175, ' is an island (no neighbors)')
('WARNING: ', 44033, ' is an island (no neighbors)')
('WARNING: ', 40188, ' is an island (no neighbors)')
('WARNING: ', 54905, ' is an island (no neighbors)')
('WARNING: ', 7696, ' is an island (no neighbors)')
('WARNING: ', 1851, ' is an island (no neighbors)')
('WARNING: ', 22112, ' is an island (no neighbors)')
('WARNING: ', 49991, ' is an island (no neighbors)')
('WARNING: ', 24705, ' is an island (no neighbors)')
('WARNING: ', 38597, ' is an island (no neighbors)')
('WARNING: ', 42334, ' is an island (no neighbors)')
('WARNING: ', 35306, ' is an island (no neighbors)')
('WARNING: ', 35261, ' is an island (no neighbors)')
('WARNING: ', 19314, ' is an island (no neighbors)')
('WARNING: ', 5784, ' is an island (no neighbors)')
('WARNING: ', 14277, ' is an island (no neighbors)')
('WARNING: ', 4828, ' is an island (no neighbors)')
('WARNING: ', 22048, ' is an island (no neighbors)')
('WARNING: ', 52782, ' is an island (no neighbors)')
('WARNING: ', 19869, ' is an island (no neighbors)')
('WARNING: ', 40010, ' is an island (no neighbors)')
('WARNING: ', 43264, ' is an island (no neighbors)')
('WARNING: ', 17367, ' is an island (no neighbors)')
('WARNING: ', 18282, ' is an island (no neighbors)')
('WARNING: ', 27950, ' is an island (no neighbors)')
('WARNING: ', 43029, ' is an island (no neighbors)')
('WARNING: ', 3180, ' is an island (no neighbors)')
('WARNING: ', 8565, ' is an island (no neighbors)')
('WARNING: ', 32426, ' is an island (no neighbors)')
('WARNING: ', 5348, ' is an island (no neighbors)')
('WARNING: ', 29338, ' is an island (no neighbors)')
('WARNING: ', 45938, ' is an island (no neighbors)')
('WARNING: ', 63145, ' is an island (no neighbors)')
('WARNING: ', 54176, ' is an island (no neighbors)')
('WARNING: ', 6188, ' is an island (no neighbors)')
('WARNING: ', 54839, ' is an island (no neighbors)')
('WARNING: ', 12912, ' is an island (no neighbors)')
('WARNING: ', 57308, ' is an island (no neighbors)')
('WARNING: ', 4816, ' is an island (no neighbors)')
('WARNING: ', 16329, ' is an island (no neighbors)')
('WARNING: ', 16914, ' is an island (no neighbors)')
('WARNING: ', 803, ' is an island (no neighbors)')
('WARNING: ', 60577, ' is an island (no neighbors)')
('WARNING: ', 49969, ' is an island (no neighbors)')
('WARNING: ', 56833, ' is an island (no neighbors)')
('WARNING: ', 5442, ' is an island (no neighbors)')
('WARNING: ', 48227, ' is an island (no neighbors)')
('WARNING: ', 37395, ' is an island (no neighbors)')
('WARNING: ', 2658, ' is an island (no neighbors)')
('WARNING: ', 48917, ' is an island (no neighbors)')
('WARNING: ', 56366, ' is an island (no neighbors)')
('WARNING: ', 55245, ' is an island (no neighbors)')
('WARNING: ', 9272, ' is an island (no neighbors)')
('WARNING: ', 21014, ' is an island (no neighbors)')
('WARNING: ', 42458, ' is an island (no neighbors)')
('WARNING: ', 29321, ' is an island (no neighbors)')
('WARNING: ', 63673, ' is an island (no neighbors)')
('WARNING: ', 23367, ' is an island (no neighbors)')
('WARNING: ', 27887, ' is an island (no neighbors)')
('WARNING: ', 5866, ' is an island (no neighbors)')
('WARNING: ', 3211, ' is an island (no neighbors)')
('WARNING: ', 22035, ' is an island (no neighbors)')
('WARNING: ', 23829, ' is an island (no neighbors)')
('WARNING: ', 51197, ' is an island (no neighbors)')
('WARNING: ', 49051, ' is an island (no neighbors)')
('WARNING: ', 33134, ' is an island (no neighbors)')
('WARNING: ', 36502, ' is an island (no neighbors)')

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/libpysal/weights/util.py:826: UserWarning: The weights matrix is not fully connected: 
 There are 1225 disconnected components.
 There are 571 islands with ids: 14865, 8557, 6769, 15460, 24589, 24581, 18310, 5382, 45376, 21821, 6288, 53823, 56887, 3308, 39216, 26231, 60581, 34683, 22120, 23900, 21844, 40189, 2791, 4448, 21588, 62377, 34004, 24983, 53791, 55983, 12892, 55989, 59932, 6234, 2527, 15318, 8492, 38604, 7107, 5517, 41014, 16140, 14624, 11310, 27248, 36391, 45846, 37112, 21808, 31079, 1736, 48743, 2541, 15969, 47256, 5540, 38807, 33648, 36209, 5476, 2182, 18307, 61370, 37151, 12900, 60893, 9670, 30717, 21018, 54196, 21838, 31522, 55603, 13010, 39485, 9093, 26791, 15302, 5547, 55728, 53157, 17599, 19321, 51062, 33652, 13835, 63592, 13459, 41015, 42118, 33636, 8375, 55941, 6268, 8925, 21841, 47813, 18441, 32163, 1805, 56839, 61971, 60492, 5528, 56855, 14898, 22130, 32159, 3452, 64506, 3431, 61338, 40568, 55725, 45990, 21597, 23945, 62597, 54184, 49169, 6235, 16339, 47189, 13941, 3882, 6293, 27888, 53042, 1057, 3940, 13905, 23822, 5718, 56868, 32706, 14448, 55945, 62576, 44861, 21854, 23113, 2592, 49176, 32660, 9203, 53796, 64534, 29115, 2528, 24598, 63874, 32042, 52069, 23801, 64731, 86, 14261, 6860, 22952, 24468, 50437, 63320, 41634, 60750, 61119, 46676, 36421, 23825, 41172, 18303, 13912, 4172, 24961, 23214, 39387, 32648, 60358, 7390, 11307, 5428, 39900, 56365, 27436, 41391, 5672, 41177, 9594, 5732, 45602, 11247, 48575, 6793, 44374, 34127, 3222, 13757, 33136, 42374, 9279, 5598, 54538, 12353, 24114, 55732, 49963, 6280, 10950, 54892, 6724, 41647, 48736, 44469, 6756, 45762, 34138, 4475, 28846, 335, 800, 57483, 33140, 6839, 15296, 17523, 42265, 6247, 60766, 8972, 49665, 24115, 27966, 63882, 56765, 21843, 7106, 21879, 12139, 50898, 21349, 8550, 5572, 22125, 58079, 49959, 1043, 9207, 43149, 63105, 39118, 30824, 64539, 2579, 27239, 59942, 22133, 45630, 57335, 62564, 55604, 33720, 1059, 55897, 31082, 39903, 56425, 13726, 32623, 60784, 36371, 46618, 44166, 41729, 44611, 55040, 55607, 52280, 11730, 1033, 462, 19988, 15364, 6821, 4935, 44034, 26517, 21726, 22395, 16174, 62362, 61968, 10852, 19366, 9069, 34463, 19922, 12362, 35235, 50433, 29874, 12858, 828, 3732, 63785, 5954, 12081, 5927, 53380, 60323, 57256, 60783, 27008, 13904, 5537, 62376, 44175, 22108, 15443, 51061, 48913, 43346, 41602, 49836, 14286, 9660, 44511, 19313, 60357, 33647, 28547, 29112, 38801, 3179, 2059, 22379, 49189, 6804, 57297, 59919, 6199, 60370, 12131, 10369, 4948, 25526, 26546, 15388, 4391, 49249, 51003, 5148, 28283, 6196, 42997, 44601, 10574, 49255, 35210, 586, 25657, 57770, 11306, 43142, 6079, 8350, 58096, 39899, 9011, 22617, 42333, 7162, 20318, 11518, 28286, 18848, 57240, 60018, 19842, 54974, 6165, 1574, 62364, 6315, 11342, 49389, 62357, 22815, 9654, 52978, 6166, 23346, 56165, 5683, 2020, 20993, 6522, 56888, 31380, 27893, 49019, 15679, 36995, 21805, 2181, 38545, 5554, 53001, 9661, 19438, 1572, 15974, 22124, 1576, 14672, 58098, 23919, 14940, 1039, 27969, 20065, 27026, 32641, 18442, 31454, 39868, 53165, 14243, 20003, 32165, 55898, 34124, 21847, 35669, 41799, 59783, 26891, 38771, 12136, 11065, 60871, 33486, 21001, 28321, 23815, 40255, 12358, 47186, 54979, 48337, 11744, 28907, 53985, 39716, 17516, 26232, 32146, 5739, 27694, 9008, 64542, 29343, 6717, 80, 10708, 41601, 28282, 33646, 19433, 24582, 49586, 4939, 31060, 64261, 27171, 37461, 5522, 7687, 60187, 5426, 54090, 10351, 27831, 31523, 12913, 32100, 60412, 45593, 39658, 21543, 58111, 16338, 17394, 38332, 2721, 32645, 35836, 4484, 61325, 38770, 49023, 13935, 37373, 57975, 16330, 29867, 23175, 44033, 40188, 54905, 7696, 1851, 22112, 49991, 24705, 38597, 42334, 35306, 35261, 19314, 5784, 14277, 4828, 22048, 52782, 19869, 40010, 43264, 17367, 18282, 27950, 43029, 3180, 8565, 32426, 5348, 29338, 45938, 63145, 54176, 6188, 54839, 12912, 57308, 4816, 16329, 16914, 803, 60577, 49969, 56833, 5442, 48227, 37395, 2658, 48917, 56366, 55245, 9272, 21014, 42458, 29321, 63673, 23367, 27887, 5866, 3211, 22035, 23829, 51197, 49051, 33134, 36502.
  w = W(neighbors, weights, ids, **kwargs)
/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/libpysal/weights/distance.py:844: UserWarning: The weights matrix is not fully connected: 
 There are 1225 disconnected components.
 There are 571 islands with ids: 14865, 8557, 6769, 15460, 24589, 24581, 18310, 5382, 45376, 21821, 6288, 53823, 56887, 3308, 39216, 26231, 60581, 34683, 22120, 23900, 21844, 40189, 2791, 4448, 21588, 62377, 34004, 24983, 53791, 55983, 12892, 55989, 59932, 6234, 2527, 15318, 8492, 38604, 7107, 5517, 41014, 16140, 14624, 11310, 27248, 36391, 45846, 37112, 21808, 31079, 1736, 48743, 2541, 15969, 47256, 5540, 38807, 33648, 36209, 5476, 2182, 18307, 61370, 37151, 12900, 60893, 9670, 30717, 21018, 54196, 21838, 31522, 55603, 13010, 39485, 9093, 26791, 15302, 5547, 55728, 53157, 17599, 19321, 51062, 33652, 13835, 63592, 13459, 41015, 42118, 33636, 8375, 55941, 6268, 8925, 21841, 47813, 18441, 32163, 1805, 56839, 61971, 60492, 5528, 56855, 14898, 22130, 32159, 3452, 64506, 3431, 61338, 40568, 55725, 45990, 21597, 23945, 62597, 54184, 49169, 6235, 16339, 47189, 13941, 3882, 6293, 27888, 53042, 1057, 3940, 13905, 23822, 5718, 56868, 32706, 14448, 55945, 62576, 44861, 21854, 23113, 2592, 49176, 32660, 9203, 53796, 64534, 29115, 2528, 24598, 63874, 32042, 52069, 23801, 64731, 86, 14261, 6860, 22952, 24468, 50437, 63320, 41634, 60750, 61119, 46676, 36421, 23825, 41172, 18303, 13912, 4172, 24961, 23214, 39387, 32648, 60358, 7390, 11307, 5428, 39900, 56365, 27436, 41391, 5672, 41177, 9594, 5732, 45602, 11247, 48575, 6793, 44374, 34127, 3222, 13757, 33136, 42374, 9279, 5598, 54538, 12353, 24114, 55732, 49963, 6280, 10950, 54892, 6724, 41647, 48736, 44469, 6756, 45762, 34138, 4475, 28846, 335, 800, 57483, 33140, 6839, 15296, 17523, 42265, 6247, 60766, 8972, 49665, 24115, 27966, 63882, 56765, 21843, 7106, 21879, 12139, 50898, 21349, 8550, 5572, 22125, 58079, 49959, 1043, 9207, 43149, 63105, 39118, 30824, 64539, 2579, 27239, 59942, 22133, 45630, 57335, 62564, 55604, 33720, 1059, 55897, 31082, 39903, 56425, 13726, 32623, 60784, 36371, 46618, 44166, 41729, 44611, 55040, 55607, 52280, 11730, 1033, 462, 19988, 15364, 6821, 4935, 44034, 26517, 21726, 22395, 16174, 62362, 61968, 10852, 19366, 9069, 34463, 19922, 12362, 35235, 50433, 29874, 12858, 828, 3732, 63785, 5954, 12081, 5927, 53380, 60323, 57256, 60783, 27008, 13904, 5537, 62376, 44175, 22108, 15443, 51061, 48913, 43346, 41602, 49836, 14286, 9660, 44511, 19313, 60357, 33647, 28547, 29112, 38801, 3179, 2059, 22379, 49189, 6804, 57297, 59919, 6199, 60370, 12131, 10369, 4948, 25526, 26546, 15388, 4391, 49249, 51003, 5148, 28283, 6196, 42997, 44601, 10574, 49255, 35210, 586, 25657, 57770, 11306, 43142, 6079, 8350, 58096, 39899, 9011, 22617, 42333, 7162, 20318, 11518, 28286, 18848, 57240, 60018, 19842, 54974, 6165, 1574, 62364, 6315, 11342, 49389, 62357, 22815, 9654, 52978, 6166, 23346, 56165, 5683, 2020, 20993, 6522, 56888, 31380, 27893, 49019, 15679, 36995, 21805, 2181, 38545, 5554, 53001, 9661, 19438, 1572, 15974, 22124, 1576, 14672, 58098, 23919, 14940, 1039, 27969, 20065, 27026, 32641, 18442, 31454, 39868, 53165, 14243, 20003, 32165, 55898, 34124, 21847, 35669, 41799, 59783, 26891, 38771, 12136, 11065, 60871, 33486, 21001, 28321, 23815, 40255, 12358, 47186, 54979, 48337, 11744, 28907, 53985, 39716, 17516, 26232, 32146, 5739, 27694, 9008, 64542, 29343, 6717, 80, 10708, 41601, 28282, 33646, 19433, 24582, 49586, 4939, 31060, 64261, 27171, 37461, 5522, 7687, 60187, 5426, 54090, 10351, 27831, 31523, 12913, 32100, 60412, 45593, 39658, 21543, 58111, 16338, 17394, 38332, 2721, 32645, 35836, 4484, 61325, 38770, 49023, 13935, 37373, 57975, 16330, 29867, 23175, 44033, 40188, 54905, 7696, 1851, 22112, 49991, 24705, 38597, 42334, 35306, 35261, 19314, 5784, 14277, 4828, 22048, 52782, 19869, 40010, 43264, 17367, 18282, 27950, 43029, 3180, 8565, 32426, 5348, 29338, 45938, 63145, 54176, 6188, 54839, 12912, 57308, 4816, 16329, 16914, 803, 60577, 49969, 56833, 5442, 48227, 37395, 2658, 48917, 56366, 55245, 9272, 21014, 42458, 29321, 63673, 23367, 27887, 5866, 3211, 22035, 23829, 51197, 49051, 33134, 36502.
  W.__init__(

# Check if the weights matrix is proper
print(f"Weights matrix properties:")
print(f"Number of observations: {knn_w.n}")
print(f"Number of neighbors per observation: {k}")

Weights matrix properties:
Number of observations: 10000
Number of neighbors per observation: 8

Now let’s calculate Global Moran’s I and interpret

y = glas_sample['glas_z']
moran = esda.Moran(y, knn_w)

print(f"Global Moran's I: {moran.I:.3f}")
print(f"p-value: {moran.p_sim:.3f}")
print(f"z-score: {moran.z_sim:.3f}")

Global Moran's I: 0.991
p-value: 0.001
z-score: 183.492

if moran.p_sim < 0.05:
    if moran.I > 0:
        print("There is statistically significant positive spatial autocorrelation in elevation values.")
        print("This means that similar elevation values tend to cluster together spatially.")
    else:
        print("There is statistically significant negative spatial autocorrelation in elevation values.")
        print("This means that dissimilar elevation values tend to be near each other (checkerboard pattern).")
else:
    print("There is no statistically significant spatial autocorrelation in elevation values.")
    print("The spatial distribution of elevation values appears to be random.")

There is statistically significant positive spatial autocorrelation in elevation values.
This means that similar elevation values tend to cluster together spatially.

f,ax=plt.subplots(figsize=(10,4))
sns.kdeplot(moran.sim,shade=True,ax=ax)
ax.axvline(moran.I,color='red',linestyle='--')
ax.vlines(moran.sim,ymin=0,ymax=3,color='black')
ax.set_xlabel('Moran\'s I')
ax.set_ylabel('Frequency')
ax.set_title('Moran\'s I Simulation Distribution')

/tmp/ipykernel_3561439/895129676.py:2: FutureWarning: 

`shade` is now deprecated in favor of `fill`; setting `fill=True`.
This will become an error in seaborn v0.14.0; please update your code.

  sns.kdeplot(moran.sim,shade=True,ax=ax)

Text(0.5, 1.0, "Moran's I Simulation Distribution")

f,ax=plt.subplots(figsize=(10,10))
moran_scatterplot(moran, p=0.05, ax=ax, aspect_equal=True)
ax.text(0.16, 0.95, "LH", transform=ax.transAxes, fontsize=25)
ax.text(0.65, 0.95, "HH", transform=ax.transAxes, fontsize=25)
ax.text(0.16, 0.05, "LL", transform=ax.transAxes, fontsize=25)
ax.text(0.65, 0.05, "HL", transform=ax.transAxes, fontsize=25)

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/splot/_viz_esda_mpl.py:131: UserWarning: `p` is only used for plotting `esda.moran.Moran_Local`
or `Moran_Local_BV` objects
  warnings.warn(

Text(0.65, 0.05, 'HL')

Spatial autocorrelation - Local Moran’s I

We’ll use the same sample and weights matrix as above

lisa = esda.Moran_Local(y, knn_w)

f,ax=plt.subplots(figsize=(10,10))
moran_scatterplot(lisa, p=0.05, ax=ax)
ax.text(1, 2, "HH", fontsize=25)
ax.text(1, -2, "HL", fontsize=25)
ax.text(-1, 2, "LH", fontsize=25)
ax.text(-1, -2, "LL", fontsize=25)

Text(-1, -2, 'LL')

# Create a classification of significant Local Moran's I results
sig = 0.05
sig_lisa = lisa.p_sim < sig
hotspot = sig_lisa * (lisa.q == 1)  # High-High clusters
coldspot = sig_lisa * (lisa.q == 3)  # Low-Low clusters
doughnut = sig_lisa * (lisa.q == 2)  # Low-High outliers
diamond = sig_lisa * (lisa.q == 4)  # High-Low outliers

# Create a new GeoDataFrame with LISA results
glas_lisa = glas_sample.copy()
glas_lisa['lisa_p'] = lisa.p_sim
glas_lisa['lisa_q'] = lisa.q
glas_lisa['significant'] = sig_lisa
glas_lisa['cluster_type'] = 'Not Significant'
glas_lisa.loc[hotspot, 'cluster_type'] = 'High-High'
glas_lisa.loc[coldspot, 'cluster_type'] = 'Low-Low'
glas_lisa.loc[doughnut, 'cluster_type'] = 'Low-High'
glas_lisa.loc[diamond, 'cluster_type'] = 'High-Low'
glas_lisa

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry	lisa_p	lisa_q	significant	cluster_type
15542	2004.429831	731738.318132	40.634266	-105.955516	3187.86	3189.93	13.01	12	POINT (446409.665 76585.637)	0.001	1	True	High-High
5692	2003.766250	731495.681086	43.334905	-120.342822	1323.71	1324.80	1.15	31	POINT (-730844.923 402143.028)	0.075	3	False	Not Significant
34538	2005.890662	732272.091524	39.276701	-106.181434	3889.66	3894.25	11.44	12	POINT (436063.724 -75975.063)	0.001	1	True	High-High
14865	2004.409427	731730.850395	45.459046	-121.662681	890.63	872.95	0.72	31	POINT (-809285.618 649025.642)	0.001	3	True	Low-Low
52951	2007.793761	732966.722938	35.047811	-119.473550	801.42	791.08	8.73	31	POINT (-745626.774 -525756.897)	0.001	3	True	Low-Low
...	...	...	...	...	...	...	...	...	...	...	...	...	...
9167	2003.853784	731527.631050	36.662171	-118.390415	3587.19	3592.78	4.88	31	POINT (-633135.494 -354993.718)	0.001	1	True	High-High
35630	2006.163054	732371.514708	38.139578	-119.725524	2536.30	2532.79	8.94	31	POINT (-735984.476 -180594.846)	0.001	1	True	High-High
42622	2006.474926	732485.347975	43.149502	-113.372769	1465.09	1463.40	0.24	31	POINT (-170174.25 346393.909)	0.010	3	True	Low-Low
40746	2006.440673	732472.845738	39.985041	-113.675834	1282.22	1283.12	0.00	31	POINT (-204187.075 -6472.453)	0.001	3	True	Low-Low
34032	2005.882575	732269.139828	47.566385	-123.303839	1592.93	1591.70	11.31	31	POINT (-904984.716 897177.191)	0.400	3	False	Not Significant

10000 rows × 13 columns

f, ax = plt.subplots(figsize=(15, 10))

# Plot only significant clusters with appropriate colors
colors = {'High-High': 'red', 'Low-Low': 'blue', 'Low-High': 'lightblue', 'High-Low': 'moccasin', 'Not Significant': 'gray'}
for cluster_type, color in colors.items():
    cluster_points = glas_lisa[glas_lisa['cluster_type'] == cluster_type]
    if len(cluster_points) > 0:
        cluster_points.plot(ax=ax, color=color, label=cluster_type, markersize=20, edgecolor='black')

ax.set_title('LISA Clusters for GLAS Elevation')
ax.legend(fontsize=12)

ctx.add_basemap(ax, crs=glas_gdf.crs, source=ctx.providers.CartoDB.Voyager)

# Calculate percentage of each cluster type
total_sig = sum(glas_lisa['significant'])
for cluster_type in ['High-High', 'Low-Low', 'Low-High', 'High-Low']:
    count = sum(glas_lisa['cluster_type'] == cluster_type)
    percent = count / total_sig * 100 if total_sig > 0 else 0
    print(f"{cluster_type}: {count} points ({percent:.1f}% of significant clusters)")

# Interpret results
print("\nInterpretation of LISA clusters:")
print("- High-High clusters (hotspots): Areas with high elevation surrounded by high elevation")
print("- Low-Low clusters (coldspots): Areas with low elevation surrounded by low elevation")
print("- Low-High outliers: Low elevation areas surrounded by high elevation (valleys)")
print("- High-Low outliers: High elevation areas surrounded by low elevation (peaks)")

High-High: 2551 points (37.6% of significant clusters)
Low-Low: 4035 points (59.4% of significant clusters)
Low-High: 0 points (0.0% of significant clusters)
High-Low: 203 points (3.0% of significant clusters)

Interpretation of LISA clusters:
- High-High clusters (hotspots): Areas with high elevation surrounded by high elevation
- Low-Low clusters (coldspots): Areas with low elevation surrounded by low elevation
- Low-High outliers: Low elevation areas surrounded by high elevation (valleys)
- High-Low outliers: High elevation areas surrounded by low elevation (peaks)

Spatial autocorrelation: Variogram analysis

Now let’s focus on a smaller region around Wyoming

#sierras_gdf = gpd.GeoDataFrame(geometry=[Polygon([(-120, 35), (-115, 35), (-115, 41), (-120, 41)])], crs='EPSG:4326') sierras_gdf = sierras_gdf.to_crs(WUS_aea_crs)
states_url = 'http://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_040_00_500k.json'
states_gdf = gpd.read_file(states_url)
wa_state_gdf = states_gdf.loc[states_gdf['NAME'] == 'Wyoming']
wa_state_gdf = wa_state_gdf.to_crs(WUS_aea_crs)
wa_glas_gdf = glas_gdf.clip(wa_state_gdf.geometry)
wa_glas_gdf

	decyear	ordinal	lat	lon	glas_z	dem_z	dem_z_std	lulc	geometry
62340	2008.957958	733392.612512	41.246035	-109.477271	1933.91	1932.70	2.55	31	POINT (149271.215 133209.437)
17403	2004.466642	731751.790861	41.254643	-109.496467	1951.60	1950.47	8.84	31	POINT (147653.616 134139.503)
17404	2004.466642	731751.790861	41.257727	-109.497060	1932.47	1935.36	8.08	31	POINT (147597.42 134483.139)
55604	2007.841519	732984.154510	41.086350	-109.433132	2122.50	2126.97	9.19	31	POINT (153312.722 115440.406)
14873	2004.410758	731731.337459	42.017654	-109.413979	1989.46	1990.01	0.10	31	POINT (152755.811 219532.112)
...	...	...	...	...	...	...	...	...	...
8927	2003.845505	731524.609145	44.977203	-108.039953	1580.09	1576.41	6.44	31	POINT (253923.854 552915.201)
36995	2006.198541	732384.467417	44.977144	-108.033952	1658.68	1661.88	11.17	31	POINT (254395.712 552925.655)
8926	2003.845505	731524.609145	44.978749	-108.039639	1666.83	1662.22	12.32	31	POINT (253942.319 553088.305)
23819	2005.176556	732011.442943	44.832693	-108.699409	1380.74	1380.29	4.05	31	POINT (202552.643 535133.893)
39862	2006.425643	732467.359674	44.833433	-108.732082	1357.49	1357.32	0.33	31	POINT (199976.024 535142.955)

5708 rows × 9 columns

f,ax=plt.subplots(figsize=(10,8))

wa_glas_gdf.plot(ax=ax,column='glas_z',cmap='inferno',legend=True,markersize=1,legend_kwds={'label': "Elevation (m)"})
ctx.add_basemap(ax,crs=wa_glas_gdf.crs.to_string())

ax.set_aspect('equal')

ax.set_title('GLAS Elevation')

f.tight_layout()

wa_sample_size = min(10000, len(wa_glas_gdf))  # Adjust based on computational resources
wa_glas_sample = wa_glas_gdf.sample(wa_sample_size, random_state=42)

Try adjusting the parameters of the variogram calculation! Try maxlag distances 1km, 5km, 50km, and ‘auto’

coords = np.vstack((wa_glas_sample.geometry.x, wa_glas_sample.geometry.y)).T
values = wa_glas_sample['glas_z'].values

# Create a variogram object
print("Calculating empirical variogram...")
V = skg.Variogram(
    coords, 
    values,
    model='spherical',  # Starting with a spherical model
    n_lags=30,          # Number of distance classes
    maxlag=5*1000,       # Maximum lag distance
    use_nugget=True     # Include a nugget effect
)

# Fit the variogram model
print("Fitting variogram model...")
V.fit()
V.plot()

Calculating empirical variogram...
Fitting variogram model...

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/skgstat/plotting/variogram_plot.py:123: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  fig.show()

# Print variogram parameters
print("Variogram Model Parameters:")
print(f"Model type: {V.describe()['model']}")
print(f"Nugget: {V.describe()['nugget']:.4f}")
print(f"Sill: {V.describe()['sill']:.4f}")
print(f"Effective range: {V.describe()['effective_range']:.4f} meters")

Variogram Model Parameters:
Model type: spherical
Nugget: 0.0000
Sill: 13055.0435
Effective range: 2367.0805 meters

# Variogram model interpretation
print("\nVariogram Interpretation:")
print(f"Nugget/Sill Ratio: {V.describe()['nugget'] / V.describe()['sill'] * 100:.1f}%")
if V.describe()['nugget'] / V.describe()['sill'] < 0.25:
    print("Low Nugget/Sill ratio indicates strong spatial dependence")
elif V.describe()['nugget'] / V.describe()['sill'] < 0.75:
    print("Moderate Nugget/Sill ratio indicates moderate spatial dependence")
else:
    print("High Nugget/Sill ratio indicates weak spatial dependence")
print(f"Range indicates that spatial dependence extends to approximately {V.describe()['effective_range']:.4f} meters")

Variogram Interpretation:
Nugget/Sill Ratio: 0.0%
Low Nugget/Sill ratio indicates strong spatial dependence
Range indicates that spatial dependence extends to approximately 2367.0805 meters

Interpretation time! Why does the range change with maxlag distance? What does this tell us about our data? For your run with the largest maxlag distance, why does the variance go back down at large lag distances??

We can try a directional variogram too!

Vnorth = skg.DirectionalVariogram(coords, values, azimuth=90, tolerance=90, maxlag=5*1000, n_lags=30)
Veast = skg.DirectionalVariogram(coords, values, azimuth=0, tolerance=90, maxlag=5*1000, n_lags=30)

f, ax = plt.subplots()
vnplot = Vnorth.plot(axes=ax)
veplot = Veast.plot(axes=ax)

ax.plot(Vnorth.bins, Vnorth.experimental, '.--r', label='North-South')
ax.plot(Veast.bins, Veast.experimental, '.--b', label='East-West')
ax.legend()

f.tight_layout()

/home/eric/miniconda3/envs/uwgda2025/lib/python3.12/site-packages/skgstat/plotting/variogram_plot.py:123: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  fig.show()

Now let’s try some spatial interpolation

First, let’s try a couple of deterministic approaches

dem_da = easysnowdata.topography.get_copernicus_dem(bbox_input=wa_state_gdf.to_crs('EPSG:4326'),resolution=90).load()
dem_da = dem_da.rio.reproject(WUS_aea_crs,resolution=(3000,3000),resampling=rasterio.enums.Resampling.average)
dem_da

<xarray.DataArray 'data' (y: 157, x: 196)> Size: 123kB
array([[      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       [      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       [      nan,       nan,       nan, ...,       nan,       nan,
              nan],
       ...,
       [2249.9893, 2225.0647, 2231.8096, ...,       nan,       nan,
              nan],
       [2214.296 , 2222.7998, 2245.6206, ...,       nan,       nan,
              nan],
       [2167.871 , 2249.182 , 2320.737 , ...,       nan,       nan,
              nan]], dtype=float32)
Coordinates:
  * x            (x) float64 2kB 1.82e+04 2.12e+04 ... 6.002e+05 6.032e+05
  * y            (y) float64 1kB 5.729e+05 5.699e+05 ... 1.079e+05 1.049e+05
    time         datetime64[ns] 8B 2021-04-22
    spatial_ref  int64 8B 0

xarray.DataArray

'data'

• y: 157
• x: 196

• nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nan

  array([[      nan,       nan,       nan, ...,       nan,       nan,
                nan],
         [      nan,       nan,       nan, ...,       nan,       nan,
                nan],
         [      nan,       nan,       nan, ...,       nan,       nan,
                nan],
         ...,
         [2249.9893, 2225.0647, 2231.8096, ...,       nan,       nan,
                nan],
         [2214.296 , 2222.7998, 2245.6206, ...,       nan,       nan,
                nan],
         [2167.871 , 2249.182 , 2320.737 , ...,       nan,       nan,
                nan]], dtype=float32)

• Coordinates: (4)
  • x
    (x)
    float64
    1.82e+04 2.12e+04 ... 6.032e+05

    axis :

    X

    long_name :

    x coordinate of projection

    standard_name :

    projection_x_coordinate

    units :

    metre

    array([ 18199.524504,  21199.524504,  24199.524504,  27199.524504,
            30199.524504,  33199.524504,  36199.524504,  39199.524504,
            42199.524504,  45199.524504,  48199.524504,  51199.524504,
            54199.524504,  57199.524504,  60199.524504,  63199.524504,
            66199.524504,  69199.524504,  72199.524504,  75199.524504,
            78199.524504,  81199.524504,  84199.524504,  87199.524504,
            90199.524504,  93199.524504,  96199.524504,  99199.524504,
           102199.524504, 105199.524504, 108199.524504, 111199.524504,
           114199.524504, 117199.524504, 120199.524504, 123199.524504,
           126199.524504, 129199.524504, 132199.524504, 135199.524504,
           138199.524504, 141199.524504, 144199.524504, 147199.524504,
           150199.524504, 153199.524504, 156199.524504, 159199.524504,
           162199.524504, 165199.524504, 168199.524504, 171199.524504,
           174199.524504, 177199.524504, 180199.524504, 183199.524504,
           186199.524504, 189199.524504, 192199.524504, 195199.524504,
           198199.524504, 201199.524504, 204199.524504, 207199.524504,
           210199.524504, 213199.524504, 216199.524504, 219199.524504,
           222199.524504, 225199.524504, 228199.524504, 231199.524504,
           234199.524504, 237199.524504, 240199.524504, 243199.524504,
           246199.524504, 249199.524504, 252199.524504, 255199.524504,
           258199.524504, 261199.524504, 264199.524504, 267199.524504,
           270199.524504, 273199.524504, 276199.524504, 279199.524504,
           282199.524504, 285199.524504, 288199.524504, 291199.524504,
           294199.524504, 297199.524504, 300199.524504, 303199.524504,
           306199.524504, 309199.524504, 312199.524504, 315199.524504,
           318199.524504, 321199.524504, 324199.524504, 327199.524504,
           330199.524504, 333199.524504, 336199.524504, 339199.524504,
           342199.524504, 345199.524504, 348199.524504, 351199.524504,
           354199.524504, 357199.524504, 360199.524504, 363199.524504,
           366199.524504, 369199.524504, 372199.524504, 375199.524504,
           378199.524504, 381199.524504, 384199.524504, 387199.524504,
           390199.524504, 393199.524504, 396199.524504, 399199.524504,
           402199.524504, 405199.524504, 408199.524504, 411199.524504,
           414199.524504, 417199.524504, 420199.524504, 423199.524504,
           426199.524504, 429199.524504, 432199.524504, 435199.524504,
           438199.524504, 441199.524504, 444199.524504, 447199.524504,
           450199.524504, 453199.524504, 456199.524504, 459199.524504,
           462199.524504, 465199.524504, 468199.524504, 471199.524504,
           474199.524504, 477199.524504, 480199.524504, 483199.524504,
           486199.524504, 489199.524504, 492199.524504, 495199.524504,
           498199.524504, 501199.524504, 504199.524504, 507199.524504,
           510199.524504, 513199.524504, 516199.524504, 519199.524504,
           522199.524504, 525199.524504, 528199.524504, 531199.524504,
           534199.524504, 537199.524504, 540199.524504, 543199.524504,
           546199.524504, 549199.524504, 552199.524504, 555199.524504,
           558199.524504, 561199.524504, 564199.524504, 567199.524504,
           570199.524504, 573199.524504, 576199.524504, 579199.524504,
           582199.524504, 585199.524504, 588199.524504, 591199.524504,
           594199.524504, 597199.524504, 600199.524504, 603199.524504])

  • y
    (y)
    float64
    5.729e+05 5.699e+05 ... 1.049e+05

    axis :

    Y

    long_name :

    y coordinate of projection

    standard_name :

    projection_y_coordinate

    units :

    metre

    array([572916.818115, 569916.818115, 566916.818115, 563916.818115,
           560916.818115, 557916.818115, 554916.818115, 551916.818115,
           548916.818115, 545916.818115, 542916.818115, 539916.818115,
           536916.818115, 533916.818115, 530916.818115, 527916.818115,
           524916.818115, 521916.818115, 518916.818115, 515916.818115,
           512916.818115, 509916.818115, 506916.818115, 503916.818115,
           500916.818115, 497916.818115, 494916.818115, 491916.818115,
           488916.818115, 485916.818115, 482916.818115, 479916.818115,
           476916.818115, 473916.818115, 470916.818115, 467916.818115,
           464916.818115, 461916.818115, 458916.818115, 455916.818115,
           452916.818115, 449916.818115, 446916.818115, 443916.818115,
           440916.818115, 437916.818115, 434916.818115, 431916.818115,
           428916.818115, 425916.818115, 422916.818115, 419916.818115,
           416916.818115, 413916.818115, 410916.818115, 407916.818115,
           404916.818115, 401916.818115, 398916.818115, 395916.818115,
           392916.818115, 389916.818115, 386916.818115, 383916.818115,
           380916.818115, 377916.818115, 374916.818115, 371916.818115,
           368916.818115, 365916.818115, 362916.818115, 359916.818115,
           356916.818115, 353916.818115, 350916.818115, 347916.818115,
           344916.818115, 341916.818115, 338916.818115, 335916.818115,
           332916.818115, 329916.818115, 326916.818115, 323916.818115,
           320916.818115, 317916.818115, 314916.818115, 311916.818115,
           308916.818115, 305916.818115, 302916.818115, 299916.818115,
           296916.818115, 293916.818115, 290916.818115, 287916.818115,
           284916.818115, 281916.818115, 278916.818115, 275916.818115,
           272916.818115, 269916.818115, 266916.818115, 263916.818115,
           260916.818115, 257916.818115, 254916.818115, 251916.818115,
           248916.818115, 245916.818115, 242916.818115, 239916.818115,
           236916.818115, 233916.818115, 230916.818115, 227916.818115,
           224916.818115, 221916.818115, 218916.818115, 215916.818115,
           212916.818115, 209916.818115, 206916.818115, 203916.818115,
           200916.818115, 197916.818115, 194916.818115, 191916.818115,
           188916.818115, 185916.818115, 182916.818115, 179916.818115,
           176916.818115, 173916.818115, 170916.818115, 167916.818115,
           164916.818115, 161916.818115, 158916.818115, 155916.818115,
           152916.818115, 149916.818115, 146916.818115, 143916.818115,
           140916.818115, 137916.818115, 134916.818115, 131916.818115,
           128916.818115, 125916.818115, 122916.818115, 119916.818115,
           116916.818115, 113916.818115, 110916.818115, 107916.818115,
           104916.818115])

  • time
    ()
    datetime64[ns]
    2021-04-22

    array('2021-04-22T00:00:00.000000000', dtype='datetime64[ns]')

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",40.0675002],PARAMETER["longitude_of_center",-111.2695312],PARAMETER["standard_parallel_1",33.5538453],PARAMETER["standard_parallel_2",46.581155],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    unknown

    horizontal_datum_name :

    World Geodetic System 1984

    projected_crs_name :

    unknown

    grid_mapping_name :

    albers_conical_equal_area

    standard_parallel :

    (33.5538453, 46.581155)

    latitude_of_projection_origin :

    40.0675002

    longitude_of_central_meridian :

    -111.2695312

    false_easting :

    0.0

    false_northing :

    0.0

    spatial_ref :

    PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",40.0675002],PARAMETER["longitude_of_center",-111.2695312],PARAMETER["standard_parallel_1",33.5538453],PARAMETER["standard_parallel_2",46.581155],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

    GeoTransform :

    16699.524503911936 3000.0 0.0 574416.8181151337 0.0 -3000.0

    array(0)

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([18199.524503911936, 21199.524503911936, 24199.524503911936,
           27199.524503911936, 30199.524503911936,  33199.52450391193,
            36199.52450391193,  39199.52450391193,  42199.52450391193,
            45199.52450391193,
           ...
             576199.524503912,   579199.524503912,   582199.524503912,
             585199.524503912,   588199.524503912,   591199.524503912,
             594199.524503912,   597199.524503912,   600199.524503912,
             603199.524503912],
          dtype='float64', name='x', length=196))

  • y
    PandasIndex

    PandasIndex(Index([ 572916.8181151337,  569916.8181151337,  566916.8181151337,
            563916.8181151337,  560916.8181151337,  557916.8181151337,
            554916.8181151337,  551916.8181151337,  548916.8181151337,
            545916.8181151337,
           ...
           131916.81811513368, 128916.81811513368, 125916.81811513368,
           122916.81811513368, 119916.81811513368, 116916.81811513368,
           113916.81811513368, 110916.81811513368, 107916.81811513368,
           104916.81811513368],
          dtype='float64', name='y', length=157))

• Attributes: (0)

x = wa_glas_gdf['geometry'].x
y = wa_glas_gdf['geometry'].y
z = wa_glas_gdf['glas_z']
bounds = wa_glas_gdf.total_bounds
dx,dy = (3000,3000)

mpl_extent = (bounds[0],bounds[2],bounds[1],bounds[3])

xi = np.arange(np.floor(bounds[0]), np.ceil(bounds[2]),dx)
yi = np.arange(np.ceil(bounds[3]),np.floor(bounds[1]),-dy)
xx, yy = np.meshgrid(xi, yi)

interp_nearest = scipy.interpolate.griddata((x,y), z, (xx, yy), method='nearest')
interp_linear = scipy.interpolate.griddata((x,y), z, (xx, yy), method='linear')
interp_cubic = scipy.interpolate.griddata((x,y), z, (xx, yy), method='cubic')
interp_func_rbf = scipy.interpolate.Rbf(x, y, z, function='linear')
interp_rbf = interp_func_rbf(xx, yy)

f,axs=plt.subplots(2,3,figsize=(20,8))
wa_glas_gdf.plot(ax=axs[0,0],column='glas_z',cmap='inferno',vmin=0,vmax=4000, markersize=2)#,legend=True,markersize=1,legend_kwds={'label': "Elevation (m)"})
axs[0,0].set_title('Original Data')
axs[0,0].set_aspect('equal')

axs[0,1].imshow(interp_nearest,extent=(bounds[0],bounds[2],bounds[1],bounds[3]),cmap='inferno',vmin=0,vmax=4000)
axs[0,1].set_title('Nearest Neighbor Interpolation')
axs[0,1].set_aspect('equal')  

axs[0,2].imshow(interp_linear,extent=(bounds[0],bounds[2],bounds[1],bounds[3]),cmap='inferno',vmin=0,vmax=4000)
axs[0,2].set_title('Linear Interpolation')
axs[0,2].set_aspect('equal')

axs[1,0].imshow(interp_cubic,extent=(bounds[0],bounds[2],bounds[1],bounds[3]),cmap='inferno',vmin=0,vmax=4000)
axs[1,0].set_title('Cubic Interpolation')
axs[1,0].set_aspect('equal')

axs[1,1].imshow(interp_rbf,extent=(bounds[0],bounds[2],bounds[1],bounds[3]),cmap='inferno',vmin=0,vmax=4000)
axs[1,1].set_title('RBF Interpolation')
axs[1,1].set_aspect('equal')

dem_da.plot.imshow(ax=axs[1,2],cmap='inferno',vmin=0,vmax=4000)
axs[1,2].set_title('DEM Data')
axs[1,2].set_aspect('equal')

for ax in axs.flat:
    wa_state_gdf.boundary.plot(ax=ax,color='blue',linewidth=2)

And now we’ll try a geostatistical method for interpolation, Ordinary Kriging

ok = skg.OrdinaryKriging(V, min_points=5, max_points=15, mode='exact') # try V or DV!

# build the target grid
coords = np.stack([wa_glas_gdf.geometry.x.values, wa_glas_gdf.geometry.y.values], axis=1)
col = 'glas_z'
vals = wa_glas_gdf[col].values
field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)
s2 = ok.sigma.reshape(xx.shape)

Warning: for 28521 locations, not enough neighbors were found within the range.

What do you think this warning means???

f, axs = plt.subplots(1, 5, figsize=(20, 6), sharey=True)

# Plot GLAS Data
wa_glas_gdf.plot(ax=axs[0], column='glas_z', cmap='inferno', vmin=0, vmax=4000, markersize=2)
axs[0].set_title('GLAS Data')
axs[0].set_aspect('equal')

# Plot DEM Data
dem_da.plot.imshow(ax=axs[1], cmap='inferno', vmin=0, vmax=4000, add_colorbar=False)
axs[1].set_title('DEM Data')
axs[1].set_aspect('equal')

# Plot RBF Interpolation
im_rbf = axs[2].imshow(interp_rbf, extent=(bounds[0], bounds[2], bounds[1], bounds[3]), cmap='inferno', vmin=0, vmax=4000)
axs[2].set_title('RBF Interpolation')
axs[2].set_aspect('equal')

# Plot Kriging Interpolation
im_kriging = axs[3].imshow(field, vmin=0, vmax=4000, cmap='inferno', extent=mpl_extent)
axs[3].set_title('Kriging Interpolation')
axs[3].set_aspect('equal')
axs[3].plot(coords[:, 0], coords[:, 1], 'k.', markersize=1)


# Plot Kriging Variance
im_variance = axs[4].imshow(s2, cmap='YlGn_r', extent=mpl_extent)
axs[4].set_title('Kriging Variance')
axs[4].set_aspect('equal')
axs[4].plot(coords[:, 0], coords[:, 1], 'k.', markersize=1)

# Add state boundary to all plots
for ax in axs.flat:
    wa_state_gdf.boundary.plot(ax=ax, color='blue', linewidth=2)
    ax.axis('off')

# Add colorbars
cbar_ax1 = f.add_axes([0.89, 0.15, 0.02, 0.7])  # Position for inferno colorbar
cbar_ax2 = f.add_axes([0.96, 0.15, 0.02, 0.7])  # Position for YlGn_r colorbar

f.colorbar(im_rbf, cax=cbar_ax1, label='Elevation (m)')
f.colorbar(im_variance, cax=cbar_ax2, label='Variance')

f.tight_layout(rect=[0, 0, 0.9, 1])  # Adjust layout to make room for colorbars

/tmp/ipykernel_3561439/3288337770.py:43: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  f.tight_layout(rect=[0, 0, 0.9, 1])  # Adjust layout to make room for colorbars

Where’s the interpolation?? Why might Kriging perform very poorly here?

Think about the range we calculated in the variogram relative to the spacing of the points! If we wanted to successfully use Kriging, we would either need to increase our spatial sampling density, or we would need the geophysical variable we are analyzing to be more correlated at longer lag distances! It also doesn’t help that our variogram’s stationarity assumption probably isn’t true.

Let’s see what the variogram and Kriging would have looked like if the elevation point values were smoothed! (This simulates positive autocorrelation at longer lag distances)

from scipy.spatial import cKDTree 

centroids = np.array([wa_glas_sample.geometry.x, wa_glas_sample.geometry.y]).T
kdtree = cKDTree(centroids)

def spatial_smoothing_gaussian(values, centroids, k=10, bandwidth=0.1):
    smoothed_values = np.zeros_like(values)
    for i, centroid in enumerate(centroids):
        distances, indices = kdtree.query(centroid, k=k)
        weights = np.exp(-distances**2 / (2 * bandwidth**2))
        smoothed_values[i] = np.sum(values[indices] * weights) / np.sum(weights)
    return smoothed_values

wa_glas_sample_smoothed = wa_glas_sample.copy()
wa_glas_sample_smoothed['glas_z'] = spatial_smoothing_gaussian(wa_glas_sample['glas_z'].values, centroids, k=1000, bandwidth=5000000)

f,axs=plt.subplots(1,2,figsize=(10,5),layout='constrained')
wa_glas_gdf.plot(ax=axs[0],column='glas_z',cmap='inferno', markersize=2, vmin=0, vmax=4000)
axs[0].set_title('Original Data')
wa_glas_sample_smoothed.plot(ax=axs[1], column='glas_z', cmap='inferno', legend=True, markersize=2,vmin=0, vmax=4000)
axs[1].set_title('Spatially Smoothed Data')

for ax in axs:
    wa_state_gdf.boundary.plot(ax=ax, color='blue', linewidth=2)
    ax.axis('off')
    ax.set_aspect('equal')

coords = np.vstack((wa_glas_sample_smoothed.geometry.x, wa_glas_sample_smoothed.geometry.y)).T
values = wa_glas_sample_smoothed['glas_z'].values

# Create a variogram object
print("Calculating empirical variogram...")
V = skg.Variogram(
    coords, 
    values,
    model='spherical',  # Starting with a spherical model
    n_lags=30,          # Number of distance classes
    maxlag='auto',#100*1000,       # Maximum lag distance
    use_nugget=True     # Include a nugget effect
)

# Fit the variogram model
print("Fitting variogram model...")
V.fit()
V.plot(show=False)

Calculating empirical variogram...
Fitting variogram model...

The Kernel crashed while executing code in the current cell or a previous cell. 

Please review the code in the cell(s) to identify a possible cause of the failure. 

Click <a href='https://aka.ms/vscodeJupyterKernelCrash'>here</a> for more info. 

View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details.

# Print variogram parameters
print("Variogram Model Parameters:")
print(f"Model type: {V.describe()['model']}")
print(f"Nugget: {V.describe()['nugget']:.4f}")
print(f"Sill: {V.describe()['sill']:.4f}")
print(f"Effective range: {V.describe()['effective_range']:.4f} meters")

Variogram Model Parameters:
Model type: spherical
Nugget: 0.0000
Sill: 774792.5096
Effective range: 628858.9018 meters

ok = skg.OrdinaryKriging(V, min_points=5, max_points=15, mode='exact')
field = ok.transform(xx.flatten(), yy.flatten()).reshape(xx.shape)
s2 = ok.sigma.reshape(xx.shape)

f,axs=plt.subplots(1,3,figsize=(20,5))

wa_glas_sample_smoothed.plot(ax=axs[0], column='glas_z', cmap='inferno', legend=True, markersize=20, edgecolor='k',vmin=0, vmax=4000, legend_kwds={'label': "Elevation (m)"})
axs[0].set_title('Smoothed GLAS Data')

im_kriging = axs[1].imshow(field, vmin=0, vmax=4000, cmap='inferno', extent=mpl_extent)
axs[1].set_title('Kriging Interpolation')
wa_glas_sample_smoothed.plot(ax=axs[1], column='glas_z', cmap='inferno', legend=True, markersize=20, edgecolor='k',vmin=0, vmax=4000, legend_kwds={'label': "Elevation (m)"})

im_variance = axs[2].imshow(s2, cmap='YlGn_r', extent=mpl_extent)
axs[2].set_title('Kriging Variance')
wa_glas_sample_smoothed.plot(ax=axs[2], color='black',markersize=5)

for ax in axs.flat:
    wa_state_gdf.boundary.plot(ax=ax, color='blue', linewidth=2)
    ax.axis('off')
    ax.set_aspect('equal')

f.tight_layout()

Looks better! At least we have predicitons now. However, the stationarity assumption is still violated because the elevation mean changes throughout the raster. When you have an underlying trend, other types of Kriging are usually more appropriate, like Universal Kriging, because they essentially assess the autocorrelation of the de-trended data. Check out the different types here.