공간 통계로 위치 데이터의 패턴과 관계를 해석하는 방법

위치가 달라지면 분석 결과도 달라진다

도시·환경·보건·네트워크 계획처럼 위치가 결과에 영향을 주는 데이터는 일반적인 통계 모델만으로 해석하기 어렵다. 가까운 관측치가 함께 높거나 낮을 수 있고, 측정하지 않은 지역의 값을 추정해야 하며, 설명변수의 영향도 지역마다 달라질 수 있다.

이때 Moran’s I는 공간적 군집과 분산을 검정하고, Kriging은 연속적인 공간 표면을 보간하며, Geographically Weighted Regression(GWR)은 관측 지점별 회귀계수의 차이를 추정한다. 세 기법은 서로 대체재라기보다 분석 목적에 따라 이어지는 도구다.

패턴 탐지, 보간, 관계 추정의 역할

Moran’s I는 연속 변수의 전역 또는 국소 공간 자기상관을 검정하는 지표다. 공간 가중행렬 W에 인접 관계나 거리 기반 관계를 담아 무작위 공간 배치와 비교하고, 관측값이 군집하는지 흩어지는지 판단한다.

Kriging은 세미베리오그램으로 공간 상관 구조를 추정한 뒤 최적 선형 불편 예측(BLUP)으로 연속 표면을 만든다. 예측값뿐 아니라 예측분산도 제공하므로, 결과의 불확실성을 함께 다룰 수 있다.

GWR은 각 관측 지점에서 국소 회귀를 수행한다. 커널과 대역폭을 통해 인접 관측치에 가중치를 주고, 변수 영향이 공간적으로 일정하지 않은 상황을 계수 변화로 드러낸다.

Moran’s I는 가중행렬 설계에서 시작한다

Moran’s I의 결과는 W를 어떻게 만들었는지에 직접 영향을 받는다. Queen·Rook 인접, 거리 임계값, K-NN 방식과 행표준화 여부는 모두 통계량과 유의성 해석에 관여한다.

전역 Moran’s I는 전체 영역의 패턴을 확인하는 데 쓰고, Local Moran’s I(LISA)는 Hot/Cold spot과 경계(Boundary)를 찾는 데 사용한다. 유의성은 퍼뮤테이션 테스트로 p-값을 추정한다.

Kriging은 변이함수와 검증 품질에 좌우된다

Kriging에서는 경험적 세미베리오그램을 만든 뒤 구형, 지수, 가우시안 모델을 적합한다. 이 과정에서 Nugget, Sill, Range 파라미터를 추정하며, 등방성 또는 이방성 여부와 방향성도 검토한다.

분석 목적에 따라 Ordinary Kriging 또는 Universal Kriging을 고른 뒤 공분산 행렬 해법으로 예측과 분산을 산출한다. LOOCV를 통해 RMSE와 MAE를 확인하고, 그 결과를 샘플링 설계 개선에 반영한다.

GWR은 지역별 계수 차이를 읽는 모델이다

GWR은 고정 또는 적응형 대역폭과 Gaussian·Bisquare 커널을 선택해 국소 가중 구조를 만든다. 대역폭은 AICc 또는 크로스밸리데이션으로 최적화한다.

적합 뒤에는 지역 R2와 잔차의 공간 자기상관을 확인해야 한다. 지역 다중공선성은 Local VIF로 점검하며, 계수 지도는 정책 타게팅이나 가설 검증에 활용할 수 있다. 다만 지역별 계수를 곧바로 인과관계로 해석해서는 안 된다.

좌표계와 관측 품질을 먼저 정리한다

거리 기반 분석은 투영 좌표(EPSG: UTM 등)를 사용하고 단위를 일관되게 맞추는 편이 낫다. CRS가 맞지 않으면 거리 왜곡이나 수렴 실패가 발생할 수 있다.

이상치와 누락값을 처리하고, 표본이 특정 지역에 지나치게 몰리지 않도록 분포를 확인한다. 측정 오차와 Nugget을 분리해 해석하는 일도 필요하다.

분야별 적용 모습

도시·부동산에서는 Moran’s I로 임대료나 가격 클러스터를 찾고, GWR로 신규 입지의 영향을 분석하며, Kriging으로 미측정 지역의 가격을 보간할 수 있다.

환경·농업에서는 Kriging으로 대기오염 또는 토양성분 농도면을 추정하고, LISA로 Hotspot 알림을 만들며, GWR 계수로 비료·관수 처방 구역을 나눌 수 있다. 보건·역학에서는 질병 발생 클러스터와 환경·사회경제 요인의 지역별 영향을 함께 검토한다.

네트워크·인프라 영역에서는 Kriging으로 커버리지 갭을 보간하고, GWR로 수요와 공급이 불균형한 지역을 진단하며, Moran’s I로 서비스 수준의 공간 패턴을 살핀다.

기대할 수 있는 분석 효과

Kriging 기반 면 추정에서는 RMSE 10~~30% 감소와 불확실성 지도 제공에 따른 샘플 재배치 효율 15~~25% 향상을 기대할 수 있다. GWR 도입 시에는 설명력(R2) 0.10~0.30p 개선과 잔차 자기상관 유의성 감소(p<0.05) 빈도 증대가 제시된다.

수치 외에도 공간 이질성을 시각화해 정책 타게팅을 정교하게 만들 수 있다. 잔차 SAC나 이방성처럼 가정 위반의 신호를 일찍 발견해 데이터 기반 의사결정의 신뢰성을 높이는 효과도 있다.

분석 흐름을 목적에 맞게 나눈다

입력에는 지점 좌표와 속성, 좌표계(CRS), 누락·이상치 상태가 포함된다. EDA와 CRS 정합 뒤 W를 구성하고 Global/Local Moran’s I를 수행한 다음, 분석 목적에 따라 보간 또는 관계 추정으로 분기한다.

연속면 예측이 필요하면 변이함수를 적합하고 Kriging 예측·분산 산출과 CV 검증을 수행한다. 인과 또는 설명이 목적이라면 대역폭 탐색 후 GWR을 적합하고, 잔차 SAC·지역 R2·VIF를 진단한다. 최종 산출물은 LISA 클러스터 지도, Kriging 예측면·분산면, GWR 계수 지도와 지역 성능 지표가 된다.

⤢✕

Python으로 핵심 분석을 실행하는 예시

전제조건은 다음과 같다.

Python 3.10+, geopandas, libpysal, esda, splot, pykrige, mgwr, shapely, numpy, scikit-learn
좌표계: 거리 연산 전 투영 좌표(EPSG: UTM/5179 등)로 변환 권장

설치:

pip install geopandas libpysal esda splot pykrige mgwr shapely numpy scikit-learn

예시 데이터는 value 필드를 가진 points.geojson이며 WGS84 좌표라고 가정한다.

Moran’s I의 전역·국소 검정은 다음과 같이 수행한다.

import geopandas as gpd
from libpysal.weights import KNN
from esda.moran import Moran, Moran_Local
import numpy as np

gdf = gpd.read_file("points.geojson").to_crs(3857)  # 투영 좌표 변환
gdf = gdf.dropna(subset=["value"])
y = gdf["value"].values

# 거리 스케일에 둔감한 KNN 가중행렬(예: k=8)
w = KNN.from_array(np.c_[gdf.geometry.x, gdf.geometry.y], k=8)
w.transform = "R"  # 행표준화

mi = Moran(y, w, permutations=999)
print("Global Moran's I:", mi.I, "p-value:", mi.p_norm, "p-perm:", mi.p_sim)

lisa = Moran_Local(y, w, permutations=999)
gdf["lisa_q"] = lisa.q    # 사분면(Hot/Cold 등)
gdf["lisa_p"] = lisa.p_sim
# 유의한 클러스터만 필터: gdf[gdf["lisa_p"] < 0.05]

Ordinary Kriging으로 예측면과 분산면을 만드는 코드는 다음과 같다.

import numpy as np
from pykrige.ok import OrdinaryKriging

# 포인트 좌표/값
x = gdf.geometry.x.values
ycoord = gdf.geometry.y.values
z = gdf["value"].values

# 그리드 정의
gridx = np.linspace(x.min(), x.max(), 100)
gridy = np.linspace(ycoord.min(), ycoord.max(), 100)

# 단순 모델(실무는 변이함수 적합/교차검증 권장)
OK = OrdinaryKriging(
    x, ycoord, z,
    variogram_model="spherical",  # "exponential","gaussian" 등
    enable_plotting=False,
    coordinates_type="euclidean",
)
z_pred, z_var = OK.execute("grid", gridx, gridy)

# z_pred, z_var를 래스터/등치선으로 시각화 가능
print(z_pred.shape, z_var.shape)

GWR에서는 설명변수와 좌표를 준비한 뒤 대역폭을 찾고 지역별 계수를 저장한다.

import numpy as np
from mgwr.gwr import GWR
from mgwr.sel_bw import Sel_BW
from sklearn.preprocessing import StandardScaler

# 설명변수 X, 종속변수 y 준비(예: x1, x2)
# gdf에 x1, x2 컬럼 존재 가정
cols = ["x1", "x2"]
gdf = gdf.dropna(subset=cols+["value"])
coords = np.c_[gdf.geometry.x.values, gdf.geometry.y.values]
y = gdf["value"].values.reshape(-1,1)
X = gdf[cols].values

# 스케일링 권장(계수 비교/수치 안정성)
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 대역폭 선택(AICc 최소)
bw = Sel_BW(coords, y, X_scaled, spherical=False).search(bw_min=30)
print("Optimal bandwidth:", bw)

model = GWR(coords, y, X_scaled, bw=bw, fixed=False, kernel='bisquare')
res = model.fit()
gdf["gwr_R2_local"] = res.localR2
for i, c in enumerate(cols):
    gdf[f"beta_{c}"] = res.params[:, i+1]  # intercept 제외

# 잔차 자기상관 재검정 권장

모델 선택 시 비교할 기준

항목	Moran’s I	Kriging	GWR
목적	공간 자기상관 탐지/검정	연속면 예측 및 불확실성 정량화	계수의 공간적 이질성 추정
입력 요건	값 + 공간가중행렬 W	좌표 + 값(연속)	좌표 + 종속 + 설명변수
모델/가정	랜덤성 대비 상관 구조	약한 정상성, 변이함수 모델	국소 정상성, 커널·대역폭
성능/확장성	O(n^2)~W희소성에 의존	O(n^3) 해법, 수천~수만 제약	수천 단위 권장, BW에 민감
일관성/안정성	퍼뮤테이션 검정으로 견고	모델 부적합 시 편향/분산 증가	다중공선성/희소구간 불안정
운영 편의	구성 단순, 해석 용이	모델링·CV 요구, 운영 난이도 중	하이퍼파라미터·진단 부담 큼

운영에서 놓치기 쉬운 진단 항목

거리 기반 모델에는 투영 좌표와 통일된 단위가 필요하다. 이상치·누락값을 처리하고, 표본이 공간적으로 균형을 이루는지도 확인한다. 가변성을 해석할 때는 Nugget과 Sill을 구분한다.

예측이 우선이면 Kriging, 설명과 정책 타게팅이 우선이면 GWR, 초기 패턴 스크리닝에는 Moran’s I가 적합하다. 이방성 또는 경계 효과가 있다면 방향성 변이함수나 배리오그램 클라우드를 다시 검토한다.

LOOCV와 외부 검증을 병행하고, 불확실성 지도를 의사결정 리스크 관리에 사용한다. GWR은 지역 다중공선성과 잔차 SAC를 반드시 진단해야 한다. Kriging은 높은 예측 정확도와 모델링 복잡도·계산비용 증가 사이의 균형이 필요하며, GWR은 높은 설명력·지역 해석력 대신 과적합·연산 비용·해석 난이도를 감수해야 한다.

CRS 불일치가 있으면 to_crs(투영) 뒤에 거리 기반 기법을 수행한다. PyKrige에서 수렴 실패나 음의 분산이 발생하면 변이함수 모델 교체, Nugget 설정, 이상치 제거를 검토한다. GWR의 지역 행렬이 특이가 되면 대역폭 증가, 변수 축소·정규화, ridge 옵션 고려(MGWR 등)가 필요하다. 과적합을 막기 위해 LOOCV/AICc 기반 하이퍼파라미터 선택과 외부 검증 세트를 확보한다.