Regresión Ponderada Geografica GWR#
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from mgwr.gwr import GWR, MGWR
from mgwr.sel_bw import Sel_BW
from matplotlib_scalebar.scalebar import ScaleBar
from mpl_toolkits.axes_grid1 import make_axes_locatable
from tabulate import tabulate
np.float = float
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 4
2 import geopandas as gpd
3 import matplotlib.pyplot as plt
----> 4 from sklearn.preprocessing import StandardScaler
5 from mgwr.gwr import GWR, MGWR
6 from mgwr.sel_bw import Sel_BW
ModuleNotFoundError: No module named 'sklearn'
gdf = gpd.read_file("https://github.com/algarciach/AnalisisGeoespacial/raw/main/Covid19_model/Data/covid19_municipios_antioquia.gpkg")
gdf.info()
<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 125 entries, 0 to 124
Data columns (total 15 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 codigo_municipio 125 non-null object
1 nombre_municipio 125 non-null object
2 codigo_subregion 125 non-null object
3 nombre_subregion 125 non-null object
4 area_municipio 125 non-null float64
5 altitud 125 non-null float64
6 temperatura 125 non-null float64
7 humedad_relativa 125 non-null float64
8 poblacion 125 non-null int64
9 urbanizacion 125 non-null float64
10 densidad 125 non-null float64
11 muertes_covid19 125 non-null int64
12 recuperados_covid19 125 non-null int64
13 cfr 125 non-null float64
14 geometry 125 non-null geometry
dtypes: float64(7), geometry(1), int64(3), object(4)
memory usage: 14.8+ KB
gdf['y_log'] = np.log(gdf['cfr'] + 1)
g_y = gdf['y_log'].values.reshape(-1, 1)
# Define coordinates as points
u = gdf.centroid.x
v = gdf.centroid.y
g_coords = list(zip(u, v))
/var/folders/nq/kj34bm5140bgwbktdx1z94cm0000gn/T/ipykernel_999/1898457282.py:2: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
u = gdf.centroid.x
/var/folders/nq/kj34bm5140bgwbktdx1z94cm0000gn/T/ipykernel_999/1898457282.py:3: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
v = gdf.centroid.y
variables = ["altitud", "temperatura", "humedad_relativa", "urbanizacion", "densidad"]
scaler = StandardScaler()
X_scaled = scaler.fit_transform(gdf[variables])
g_X_num = np.column_stack((np.ones(125), X_scaled))
g_coords = np.array(g_coords)
gdf = gpd.GeoDataFrame(
geometry=gpd.points_from_xy(g_coords[:, 0], g_coords[:, 1]),
crs="EPSG:4326"
)
gdf = gdf.to_crs("EPSG:32618")
g_coords_utm = np.array([[p.x, p.y] for p in gdf.geometry])
print("g_coords_utm shape:", g_coords_utm.shape)
g_coords_utm shape: (125, 2)
model = GWR(g_coords_utm, g_y, g_X_num, bw=12500, fixed=True, kernel='gaussian').fit()
print(f"UTM BW=12,500 m, R²: {model.R2}, AIC: {model.aic}")
model.summary()
UTM BW=12,500 m, R²: 0.8811742465106677, AIC: 11.443115021614148
===========================================================================
Model type Gaussian
Number of observations: 125
Number of covariates: 6
Global Regression Results
---------------------------------------------------------------------------
Residual sum of squares: 12.537
Log-likelihood: -33.639
AIC: 79.277
AICc: 82.234
BIC: -562.033
R2: 0.197
Adj. R2: 0.164
Variable Est. SE t(Est/SE) p-value
------------------------------- ---------- ---------- ---------- ----------
X0 1.324 0.029 45.606 0.000
X1 -0.249 0.210 -1.188 0.235
X2 -0.124 0.211 -0.587 0.557
X3 -0.050 0.033 -1.531 0.126
X4 -0.009 0.032 -0.282 0.778
X5 -0.045 0.033 -1.354 0.176
Geographically Weighted Regression (GWR) Results
---------------------------------------------------------------------------
Spatial kernel: Fixed gaussian
Bandwidth used: 12500.000
Diagnostic information
---------------------------------------------------------------------------
Residual sum of squares: 1.856
Effective number of parameters (trace(S)): 90.475
Degree of freedom (n - trace(S)): 34.525
Sigma estimate: 0.232
Log-likelihood: 85.754
AIC: 11.443
AICc: 531.611
BIC: 270.164
R2: 0.881
Adjusted R2: 0.560
Adj. alpha (95%): 0.003
Adj. critical t value (95%): 2.995
Summary Statistics For GWR Parameter Estimates
---------------------------------------------------------------------------
Variable Mean STD Min Median Max
-------------------- ---------- ---------- ---------- ---------- ----------
X0 2.309 9.755 -40.998 1.201 60.584
X1 -0.162 7.189 -22.123 -0.173 42.306
X2 -0.094 6.447 -20.837 -0.174 35.768
X3 -0.193 0.651 -2.593 -0.118 1.027
X4 -0.020 0.505 -2.482 -0.016 1.955
X5 5.345 44.285 -187.444 -0.032 282.946
===========================================================================
Calculate MWGR#
selector = Sel_BW(g_coords_utm, g_y, g_X_num, multi=True, fixed=True, kernel='gaussian', spherical=False)
raw_bw = selector.search(criterion='AICc')
print("Bandas sin restricción:", raw_bw)
# Forzar límites de 10 km (10,000 m) a 150 km (150,000 m)
bw_min = 10000
bw_max = 150000
bw_clipped = np.clip(raw_bw, bw_min, bw_max)
print("Bandas ajustadas al rango [10km, 150km]:", bw_clipped)
Bandas sin restricción: [751688.53 751657.91 147042.79 751688.53 751688.53 751367.84]
Bandas ajustadas al rango [10km, 150km]: [150000. 150000. 147042.79 150000. 150000. 150000. ]
model = MGWR(
g_coords_utm,
g_y,
g_X_num,
selector=selector,
fixed=True
).fit()
model.summary()
===========================================================================
Model type Gaussian
Number of observations: 125
Number of covariates: 6
Global Regression Results
---------------------------------------------------------------------------
Residual sum of squares: 12.537
Log-likelihood: -33.639
AIC: 79.277
AICc: 82.234
BIC: -562.033
R2: 0.197
Adj. R2: 0.164
Variable Est. SE t(Est/SE) p-value
------------------------------- ---------- ---------- ---------- ----------
X0 1.324 0.029 45.606 0.000
X1 -0.249 0.210 -1.188 0.235
X2 -0.124 0.211 -0.587 0.557
X3 -0.050 0.033 -1.531 0.126
X4 -0.009 0.032 -0.282 0.778
X5 -0.045 0.033 -1.354 0.176
Multi-Scale Geographically Weighted Regression (MGWR) Results
---------------------------------------------------------------------------
Spatial kernel: Fixed bisquare
Criterion for optimal bandwidth: AICc
Score of Change (SOC) type: Smoothing f
Termination criterion for MGWR: 1e-05
MGWR bandwidths
---------------------------------------------------------------------------
Variable Bandwidth ENP_j Adj t-val(95%) Adj alpha(95%)
X0 751688.530 1.021 1.988 0.049
X1 751657.910 1.018 1.987 0.049
X2 147042.790 4.785 2.600 0.010
X3 751688.530 1.029 1.992 0.049
X4 751688.530 1.047 1.999 0.048
X5 751367.840 1.001 1.980 0.050
Diagnostic information
---------------------------------------------------------------------------
Residual sum of squares: 12.426
Effective number of parameters (trace(S)): 9.900
Degree of freedom (n - trace(S)): 115.100
Sigma estimate: 0.329
Log-likelihood: -33.086
AIC: 87.973
AICc: 90.267
BIC: 118.802
R2 0.204
Adjusted R2 0.135
Summary Statistics For MGWR Parameter Estimates
---------------------------------------------------------------------------
Variable Mean STD Min Median Max
-------------------- ---------- ---------- ---------- ---------- ----------
X0 1.318 0.000 1.318 1.318 1.318
X1 -0.058 0.000 -0.059 -0.058 -0.057
X2 0.064 0.012 0.046 0.062 0.107
X3 -0.041 0.000 -0.041 -0.041 -0.040
X4 -0.012 0.000 -0.013 -0.012 -0.012
X5 -0.050 0.000 -0.050 -0.050 -0.050
===========================================================================
local_params = model.params
param_names = ["altitud", "temperatura", "humedad_relativa", "urbanizacion", "densidad"]
fig, axes = plt.subplots(1, 5, figsize=(18, 6))
data_clean = gdf.to_crs(epsg=4326)
for i, param_name in enumerate(param_names):
ax = axes[i]
data_clean[f'coef_{param_name}'] = local_params[:, i+1]
data_clean.plot(column=f'coef_{param_name}', ax=ax, cmap=plt.cm.jet, legend=True)
scalebar = ScaleBar(111319.49079327357, "m", location='lower right', scale_loc='top', length_fraction=0.25, font_properties={"size": 10})
ax.add_artist(scalebar)
plt.tight_layout()
plt.show()
/usr/local/Caskroom/miniforge/base/envs/geo/lib/python3.11/site-packages/matplotlib_scalebar/scalebar.py:457: UserWarning: Drawing scalebar on axes with unequal aspect ratio; either call ax.set_aspect(1) or suppress the warning with rotation='horizontal-only'.
warnings.warn(
/usr/local/Caskroom/miniforge/base/envs/geo/lib/python3.11/site-packages/matplotlib_scalebar/scalebar.py:457: UserWarning: Drawing scalebar on axes with unequal aspect ratio; either call ax.set_aspect(1) or suppress the warning with rotation='horizontal-only'.
warnings.warn(
/usr/local/Caskroom/miniforge/base/envs/geo/lib/python3.11/site-packages/matplotlib_scalebar/scalebar.py:457: UserWarning: Drawing scalebar on axes with unequal aspect ratio; either call ax.set_aspect(1) or suppress the warning with rotation='horizontal-only'.
warnings.warn(
/usr/local/Caskroom/miniforge/base/envs/geo/lib/python3.11/site-packages/matplotlib_scalebar/scalebar.py:457: UserWarning: Drawing scalebar on axes with unequal aspect ratio; either call ax.set_aspect(1) or suppress the warning with rotation='horizontal-only'.
warnings.warn(
fig, ax = plt.subplots(figsize=(10, 10))
data_clean['fitted'] = model.predy
data_clean.plot(column='fitted', ax=ax, cmap=plt.cm.jet, legend=True)
scalebar = ScaleBar(111319.49079327357, "m", location='lower right', scale_loc='top', length_fraction=0.25, font_properties={"size": 10})
ax.add_artist(scalebar)
plt.tight_layout()
plt.show()
/home/C0d3l7/MyCodes/Python/Geospatial/.geo/lib/python3.10/site-packages/matplotlib_scalebar/scalebar.py:457: UserWarning: Drawing scalebar on axes with unequal aspect ratio; either call ax.set_aspect(1) or suppress the warning with rotation='horizontal-only'.
warnings.warn(
