Moran’s I와 Kriging으로 공간 패턴을 진단하는 GIS 분석

위치 관계가 분석 결과를 바꾸는 지점

공간 통계는 관측값이 놓인 위치와 값 사이의 의존성을 수리적으로 다룬다. 값의 분포가 우연한지, 인접 지역끼리 비슷한 값을 보이는지 진단하고, 관측되지 않은 영역은 연속면으로 추정한다.

이 과정에서 Moran’s I는 공간 자기상관을 확인하는 도구이고, Kriging은 관측 지점 사이를 보간하는 방법이다. GIS 데이터 분석에서는 좌표계, 거리 단위, 공간 가중치 설계가 두 방법 모두의 전제가 된다.

Moran’s I로 전역과 국지 패턴 읽기

Moran’s I는 관측값이 공간적으로 유사하거나 상이하게 분포하는 정도를 [-1, 1] 구간에서 나타내는 전역(Globlal) 통계량이다. 양(+) 값은 유사한 값이 이웃해 군집화되는 경향을, 음(-) 값은 서로 다른 값이 인접하는 분산 경향을 뜻한다.

유의성은 무작위 재배열(permutation)로 p-값을 산출해 검정한다. 지역(Local) 확장인 LISA(Local Moran’s I)는 핫/콜드스팟과 고-저, 저-고 이웃 조합을 찾아내는 데 사용한다.

거리별 변이를 모델에 담는 방법

반변량함수(semivariogram)는 거리에 따라 달라지는 분산, 즉 공분산의 역을 경험적으로 추정한다. 구형·지수·가우시안 등의 모형을 적합하며, 너깃(nugget), 범위(range), 천이(sill)가 주요 모수다.

Kriging은 이 공분산 구조를 바탕으로 하는 선형 최소분산 불편추정(BLUP)이다. Ordinary, Universal, Simple Kriging으로 구분되며, 국지 이웃(neighborhood)을 활용해 예측값과 변동성(분산)을 함께 추정한다.

GIS 데이터와 공간 가중치가 만드는 분석 조건

벡터와 래스터 데이터는 좌표참조계(CRS)를 일치시켜야 하며, 거리 기반 분석에서는 거리 단위를 다룰 수 있는 투영 좌표가 필요하다. EPSG:5179 같은 투영 좌표 사용이 권장된다.

공간 가중치 행렬 W는 Moran’s I, 공간 회귀, 공간 검정의 핵심 입력이다. 인접 관계는 Queen/Rook 접촉 방식으로, 거리 관계는 거리 기반 또는 knn 방식으로 정의할 수 있다.

분석 파이프라인에서 점검할 설정

데이터 전처리에서는 CRS를 통일하고 이상치·결측치와 경계효과(border effect)를 확인한다. 샘플 밀도가 불균형하면 가중 샘플링 또는 다중 해상도 분석을 적용할 수 있다. 로그 변환이나 Box-Cox 변환은 정규성과 분산 안정화에 활용한다.

이웃은 knn(k=6~12) 또는 거리 임계값으로 설정하며, 희소(sparse) 구조와 행 표준화(row-standardized)를 권장한다. k나 거리 컷오프를 바꿨을 때 Moran’s I가 얼마나 안정적인지 민감도 분석으로 살핀다.

경험 반변량함수를 계산한 뒤 구형·지수·가우시안 모형을 적합한다. 잔차와 Cressie-Hawkins 등의 견고 추정도 적용 대상이다. 방향성, 즉 비등방성이 드러나면 방향별 변이도와 회전 모수를 설정한다.

Kriging에서는 탐색 반경, 최소·최대 이웃 수, 블록 크기(블록 Kriging)를 정한다. LOOCV 또는 K-fold 교차검증으로 RMSE·MAE와 예측분산의 보정 계수를 점검한다. 잔차의 Moran’s I는 공간 구조가 모형에 충분히 반영됐는지 확인하는 데 쓰이며, 결과에는 예측분산 기반의 불확실성 맵을 포함한다.

데이터 수집부터 검정, 보간, 평가까지를 자동화하고 환경·버전을 고정하면 재현성과 감사추적성을 확보할 수 있다.

공간 분포를 해석하고 운영에 연결하는 사례

도시 환경 품질 지도화에서는 미세먼지와 열섬 지표의 핫스팟을 Moran’s I로 진단하고, Kriging으로 고해상도 분포를 예측한다. 대기질 센서의 간헐적 결측을 보완하고 취약지역을 선제적으로 관리하는 데 연결할 수 있다.

수자원·토양 관리에서는 지하수 질산염과 토양 유기물 함량의 공간 분포를 추정한다. 시료 채취 지점을 최적화하고 가변율(Variable-rate) 처방을 만드는 데 활용된다.

공중보건·안전 분야에서는 병원성 발생률의 공간 자기상관을 검정하고 위험도 연속면을 산출한다. 그 결과는 자원 배치 우선순위와 모니터링 구간을 정하는 근거가 된다.

IDW 대비 교차검증 RMSE는 1030% 감소할 수 있으며, 예측분산 기반으로 샘플링을 재설계하면 추가 515% 개선할 수 있다. 변이도 구조가 안정적일 때는 샘플 수를 20~40% 절감하면서 동등 RMSE를 달성할 수 있다. 핫/콜드스팟과 불확실성을 함께 제시하면 의사결정 신뢰도를 높이고, 표준화된 파이프라인은 재현성과 감사추적성에도 기여한다.

관측 데이터에서 결과 맵까지

입력은 속성값을 포함한 포인트 관측, 경계 폴리곤, CRS 정보다. 전처리 후 W를 구성하고 Moran’s I 전역·지역 검정을 수행한다. 이어 반변량함수를 추정·적합한 뒤 Kriging 보간과 교차검증을 진행한다.

비정상성(추세)이 발견되면 Universal Kriging을 사용한다. 모형 적합이 실패하면 견고 변이도나 대체 모형을 검토하고, 과적합이 나타나면 이웃 수와 반경을 규제한다. 최종 산출물은 클러스터 맵, 예측면, 예측분산, RMSE·MAE 평가 지표, 잔차 Moran’s I다.

⤢✕

PySAL과 PyKrige로 구현하는 분석 흐름

분석 환경은 Python 3.10+, geopandas 0.14+, libpysal 4.10+, esda 2.5+, pykrige 1.7+, numpy, scikit-learn(옵션), shapely 2.0+를 전제로 한다. 입력 데이터는 value 컬럼을 가진 관측 포인트 GeoJSON이며 좌표계는 임의일 수 있다. 분석 시에는 EPSG:5179(Korea 2000 / Unified CS)로 투영한다.

# pip install geopandas libpysal esda pykrige scikit-learn numpy shapely
import numpy as np
import geopandas as gpd
from libpysal.weights import KNN
from esda.moran import Moran, Moran_Local
from pykrige.ok import OrdinaryKriging
from sklearn.model_selection import KFold
from sklearn.metrics import mean_squared_error

# 1) 데이터 로드 및 투영
gdf = gpd.read_file("samples.geojson").to_crs(5179)  # meter 단위
gdf = gdf.dropna(subset=["value"]).copy()
gdf["x"] = gdf.geometry.x
gdf["y"] = gdf.geometry.y

# 2) 공간 가중치(W) 구성 및 Moran's I
w = KNN.from_dataframe(gdf, k=8)
w.transform = "R"  # 행 표준화
y = gdf["value"].values
mi = Moran(y, w, permutations=999)
print(f"Global Moran's I: {mi.I:.4f}, p-value: {mi.p_sim:.4f}")

li = Moran_Local(y, w, permutations=999)
gdf["lisa_p"] = li.p_sim
gdf["lisa_q"] = li.q  # 1: HH, 2: LH, 3: LL, 4: HL

# 3) Kriging 설정 및 그리드 예측
x, ycoord, z = gdf["x"].values, gdf["y"].values, gdf["value"].values
OK = OrdinaryKriging(
    x, ycoord, z,
    variogram_model="spherical",  # 'exponential', 'gaussian' 가능
    enable_plotting=False, verbose=False
)

# 그리드 정의(해상도는 데이터 밀도에 맞춰 조정)
nx, ny = 200, 200
xmin, ymin, xmax, ymax = gdf.total_bounds
gridx = np.linspace(xmin, xmax, nx)
gridy = np.linspace(ymin, ymax, ny)
zgrid, vgrid = OK.execute("grid", gridx, gridy)  # zgrid: 예측, vgrid: 분산

# 4) 간단 K-fold 교차검증
kf = KFold(n_splits=5, shuffle=True, random_state=42)
preds, trues = [], []
for train_idx, test_idx in kf.split(gdf):
    OK_cv = OrdinaryKriging(
        x[train_idx], ycoord[train_idx], z[train_idx],
        variogram_model="spherical", enable_plotting=False, verbose=False
    )
    z_pred, _ = OK_cv.execute("points", x[test_idx], ycoord[test_idx])
    preds.extend(np.asarray(z_pred).ravel())
    trues.extend(z[test_idx])

rmse = mean_squared_error(trues, preds, squared=False)
print(f"Kriging 5-fold RMSE: {rmse:.3f}")

# 5) 결과 활용(예: 넘파이 저장; 래스터 저장은 rasterio 사용)
np.save("kriging_pred.npy", np.asarray(zgrid))
np.save("kriging_var.npy", np.asarray(vgrid))
gdf.to_file("lisa_result.geojson", driver="GeoJSON")

데이터에 뚜렷한 추세가 있으면 추세항 1차/2차를 둔 Universal Kriging으로 전환한다. 비등방성이 감지되면 좌표 회전과 방향별 모수 설정을 적용한다. 성능 문제가 생기면 탐색 반경과 이웃 수를 제한하고, 그리드 해상도를 낮추거나 타일 처리한다.

Moran’s I와 Kriging이 맡는 역할

운영 관점	Moran’s I	Kriging
성능(연산 복잡도)	O(N·k), 희소 W로 빠름	O(N^3) 경향, 국지 이웃/타일링으로 완화
확장성(대용량)	수십~수백만 포인트까지 가능(스트리밍/샘플링)	수만~수십만 포인트에서 국지 Kriging 권장
일관성(추정 안정성)	가중치 설계 민감도 중간	변이도 모형·샘플 밀도에 고민감
안정성(모형 민감도)	permutation 검정으로 견고	모형 부정합 시 과소/과대 분산 가능
운영 편의	설정 단순(W만 정의)	모수 다수(변이도, 이웃, 반경)·검증 필수

Moran’s I는 공간 자기상관의 유무와 패턴을 정량적으로 진단하고, LISA는 국지 클러스터를 파악한다. 변이도 구조가 확인되면 Kriging으로 연속면과 불확실성을 함께 산출할 수 있다. 전처리→W→검정→변이도→Kriging→검증→리포트 순서의 표준 파이프라인을 두고, 환경 고정과 민감도 분석을 지속하는 방식이 운영에 적합하다.