import numpy as np, pandas as pd, matplotlib.pyplot as plt, geopandas as gpd, contextily as cx
import pointproc as P
rng = np.random.default_rng(1)
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; GREY="#cccccc"
pts = pd.read_csv("bei_points.csv").to_numpy(float); gcov = pd.read_csv("bei_grid.csv")
xr=(0,1000); yr=(0,500)
gr = P.make_grid(pts, gcov["x"].to_numpy(), gcov["y"].to_numpy(), gcov[["elev","grad"]].to_numpy(), xr, yr, 25.0)
nx,ny=gr["nx"],gr["ny"]
def gridmap(ax, vals, title, cmap="viridis", pts_overlay=False):
    im=ax.imshow(vals.reshape(nx,ny).T, origin="lower", extent=[*xr,*yr], cmap=cmap, aspect="equal")
    if pts_overlay: ax.scatter(pts[:,0],pts[:,1],s=1,color="k",alpha=.3)
    ax.set_title(title,fontsize=10); ax.set_xticks([]); ax.set_yticks([]); plt.colorbar(im,ax=ax,shrink=0.5)
print(f"{len(pts)} trees on a {xr[1]}x{yr[1]} m plot; overall intensity {len(pts)/(xr[1]*yr[1])*1e4:.1f} per hectare")
print("A point process models the pattern of locations itself -- random, clustered, or regular?")

3604 trees on a 1000x500 m plot; overall intensity 72.1 per hectare
A point process models the pattern of locations itself -- random, clustered, or regular?

# locator: where the BCI plot is (a point on an OpenStreetMap basemap of the Panama Canal region)
loc = gpd.GeoDataFrame(geometry=gpd.points_from_xy([-79.8461],[9.1543]), crs=4326).to_crs(3857)
fig,ax=plt.subplots(figsize=(5.6,5.2))
cxc,cyc=float(loc.geometry.x.iloc[0]),float(loc.geometry.y.iloc[0]); dd=22000
ax.set_xlim(cxc-dd,cxc+dd); ax.set_ylim(cyc-dd,cyc+dd)
loc.plot(ax=ax, color=RED, marker="*", markersize=400, edgecolor="k", zorder=5)
try:
    cx.add_basemap(ax, source=cx.providers.OpenStreetMap.Mapnik, crs=3857, attribution_size=5)
except Exception as e:
    print("basemap tiles unavailable:", str(e)[:80])
ax.set_title("Barro Colorado Island 50-ha plot (star) — Gatún Lake, Panama Canal"); ax.set_axis_off()
plt.tight_layout(); plt.show()
print("The plot sits on an island in the Panama Canal. The tree coordinates in the maps below are LOCAL plot metres,")
print("not lat/long, so those maps use the 1000 x 500 m plot frame rather than this geographic basemap.")

The plot sits on an island in the Panama Canal. The tree coordinates in the maps below are LOCAL plot metres,
not lat/long, so those maps use the 1000 x 500 m plot frame rather than this geographic basemap.

fig,ax=plt.subplots(1,2,figsize=(9.8,2.7))
ax[0].scatter(pts[:,0],pts[:,1],s=2,color=GREEN); ax[0].set_aspect("equal"); ax[0].set_title(f"{len(pts)} tree locations"); ax[0].set_xticks([]);ax[0].set_yticks([])
gridmap(ax[1], gr["cov"][:,0], "Elevation (m), with trees", cmap="terrain", pts_overlay=True)
plt.tight_layout(); plt.show()
print("The pattern looks clumped and terrain-linked -- but 'looks clustered' needs a test against complete spatial")
print("randomness, and the terrain link needs a model. Both follow.")

The pattern looks clumped and terrain-linked -- but 'looks clustered' needs a test against complete spatial
randomness, and the terrain link needs a model. Both follow.

radii=np.linspace(2,100,25)
L=P.ripley_L(pts,xr,yr,radii); lo,hi=P.csr_envelope(len(pts),xr,yr,radii,rng,nsim=39)
fig,ax=plt.subplots(figsize=(5.6,3.2))
ax.fill_between(radii,lo,hi,color=GREY,alpha=.7,label="CSR envelope"); ax.plot(radii,L,color=RED,lw=2,label="observed L(r)-r")
ax.axhline(0,color="k",lw=.8); ax.set_xlabel("distance r (m)"); ax.set_ylabel("L(r) - r")
ax.set_title("Ripley's L: trees are strongly clustered (curve far above CSR)"); ax.legend(frameon=False,fontsize=8)
plt.tight_layout(); plt.show()
print(f"The observed L sits far ABOVE the CSR envelope at every scale (L(30)-30 = {L[np.argmin(abs(radii-30))]:.0f} m): the pattern")
print("is strongly clustered, not random. Clustering can come from covariates (terrain) OR from unmeasured spatial")
print("structure -- the next two models separate the two.")

The observed L sits far ABOVE the CSR envelope at every scale (L(30)-30 = 17 m): the pattern
is strongly clustered, not random. Clustering can come from covariates (terrain) OR from unmeasured spatial
structure -- the next two models separate the two.

elev=(gr["cov"][:,0]-gr["cov"][:,0].mean())/gr["cov"][:,0].std()
grad=(gr["cov"][:,1]-gr["cov"][:,1].mean())/gr["cov"][:,1].std()
X=np.column_stack([np.ones(len(elev)),elev,grad])
b,se=P.inhom_poisson(gr["counts"],X,gr["area"])
print(f"inhomogeneous Poisson intensity ~ elevation + gradient:")
print(f"  elevation {b[1]:+.2f} (se {se[1]:.2f})   gradient {b[2]:+.2f} (se {se[2]:.2f})   [both 'significant' -- but see below]")
fitted=np.exp(X@b)*gr["area"]
fig,ax=plt.subplots(1,2,figsize=(9.8,2.7))
gridmap(ax[0], gr["counts"], "Observed tree counts per cell", cmap="viridis")
gridmap(ax[1], fitted, "Fitted intensity (covariates only)", cmap="viridis")
plt.tight_layout(); plt.show()
print("Trees are denser at higher elevation and on steeper slopes, and the fitted surface captures the broad gradient.")
print("But it is SMOOTH -- it cannot reproduce the fine clumping, and (crucially) its standard errors assume the points")
print("are independent given the covariates, which the Ripley test already told us is false. Hence the Cox process.")

inhomogeneous Poisson intensity ~ elevation + gradient:
  elevation +0.16 (se 0.02)   gradient +0.33 (se 0.02)   [both 'significant' -- but see below]

Trees are denser at higher elevation and on steeper slopes, and the fitted surface captures the broad gradient.
But it is SMOOTH -- it cannot reproduce the fine clumping, and (crucially) its standard errors assume the points
are independent given the covariates, which the Ripley test already told us is false. Hence the Cox process.

res=P.lgcp_gibbs(gr["counts"],X,gr["area"],gr["W"].astype(float),rng,draws=2000,burn=2000)
be=res["beta"].mean(0); bl,bh=np.percentile(res["beta"],[2.5,97.5],axis=0); sd=res["spatial_sd"].mean()
print(f"LGCP:  elevation {be[1]:+.2f} [{bl[1]:+.2f},{bh[1]:+.2f}]   gradient {be[2]:+.2f} [{bl[2]:+.2f},{bh[2]:+.2f}]")
print(f"  latent spatial-field SD = {sd:.2f}  -> large residual clustering NOT explained by terrain")
print(f"  covariate CrIs are far wider than the naive Poisson SEs (~{se[1]:.2f}): the LGCP corrects the overconfidence.")
phi=res["phi"].mean(0); tot=np.exp(X@be+phi)*gr["area"]
fig,ax=plt.subplots(1,2,figsize=(9.8,2.7))
gridmap(ax[0], phi, "Latent spatial field (residual clustering)", cmap="RdBu_r")
gridmap(ax[1], tot, "LGCP total intensity (covariates + field)", cmap="viridis")
plt.tight_layout(); plt.show()
print("The latent field (left) shows strong clumps the covariates missed; the total-intensity surface (right) now")
print("reproduces the real clustering. The LGCP's message: terrain explains part of where the trees are, but most of")
print("the pattern is spatial clustering (SD %.2f) beyond it -- and ignoring it makes a plain Poisson overconfident."%sd)

LGCP:  elevation +0.41 [+0.16,+0.66]   gradient +0.47 [+0.34,+0.59]
  latent spatial-field SD = 1.28  -> large residual clustering NOT explained by terrain
  covariate CrIs are far wider than the naive Poisson SEs (~0.02): the LGCP corrects the overconfidence.

The latent field (left) shows strong clumps the covariates missed; the total-intensity surface (right) now
reproduces the real clustering. The LGCP's message: terrain explains part of where the trees are, but most of
the pattern is spatial clustering (SD 1.28) beyond it -- and ignoring it makes a plain Poisson overconfident.

import pymc as pm
adj=gr["W"].astype(int); off=np.log(gr["area"]); y=gr["counts"].astype(int)
with pm.Model() as mod:
    b0=pm.Normal("b0",0,5); bc=pm.Normal("bc",0,5,shape=2); tau=pm.HalfNormal("tau",2)
    phi_=pm.ICAR("phi",W=adj)
    eta=off+b0+bc[0]*elev+bc[1]*grad+tau*phi_
    pm.Poisson("y",mu=pm.math.exp(eta),observed=y)
    idata=pm.sample(3000,tune=4000,chains=4,target_accept=0.99,random_seed=3,progressbar=False)
bc_pm=idata.posterior["bc"].mean(("chain","draw")).values; tau_pm=float(idata.posterior["tau"].mean())
import arviz as az
def diag(names):
    r = max(float(v.max()) for v in az.rhat(idata, var_names=names).data_vars.values())
    e = min(float(v.min()) for v in az.ess(idata, var_names=names).data_vars.values())
    return r, e
rs, es = diag(["b0","bc","tau"]); rf, ef = diag(["phi"])
print(f"PyMC convergence -- reported quantities (b0, coefficients, field scale): max r-hat {rs:.3f}, min ESS {es:.0f}")
print(f"                   the {idata.posterior['phi'].shape[-1]} latent cells individually:            max r-hat {rf:.3f}, min ESS {ef:.0f}")
worst = max((("b0","bc","tau")), key=lambda v: float(az.rhat(idata, var_names=[v])[v].max()))
print(f"   the slowest-mixing parameter is {worst!r} -- the intercept, the covariate coefficients and the field")
print("   scale all trade off against the level of the latent field, which is what limits mixing here rather")
print("   than the 800 cells themselves. Both sets are close enough to 1.01 for the comparison below to stand.")
print(f"elevation: from-scratch {be[1]:+.2f}   PyMC {bc_pm[0]:+.2f}")
print(f"gradient : from-scratch {be[2]:+.2f}   PyMC {bc_pm[1]:+.2f}")
print(f"spatial-field SD: from-scratch {sd:.2f}   PyMC {tau_pm:.2f}")
print("PyMC's ICAR-LGCP agrees: a strong latent field on top of a weaker covariate signal. Point-process intensity")
print("estimation and areal disease-mapping are the same Poisson-plus-spatial-field model, one on points, one on regions.")

g++ not available, if using conda: `conda install gxx`

Initializing NUTS using jitter+adapt_diag...

Multiprocess sampling (4 chains in 4 jobs)

NUTS: [b0, bc, tau, phi]

Sampling 4 chains for 4_000 tune and 3_000 draw iterations (16_000 + 12_000 draws total) took 40 seconds.

The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

PyMC convergence -- reported quantities (b0, coefficients, field scale): max r-hat 1.017, min ESS 200
                   the 800 latent cells individually:            max r-hat 1.015, min ESS 288
   the slowest-mixing parameter is 'bc' -- the intercept, the covariate coefficients and the field
   scale all trade off against the level of the latent field, which is what limits mixing here rather
   than the 800 cells themselves. Both sets are close enough to 1.01 for the comparison below to stand.
elevation: from-scratch +0.41   PyMC +0.47
gradient : from-scratch +0.47   PyMC +0.51
spatial-field SD: from-scratch 1.28   PyMC 1.29
PyMC's ICAR-LGCP agrees: a strong latent field on top of a weaker covariate signal. Point-process intensity
estimation and areal disease-mapping are the same Poisson-plus-spatial-field model, one on points, one on regions.

# R's kppm reports sigma2 = 1.58 and a clustering scale of 48 m (see pointproc_R.ipynb);
# its trend coefficients are in raw units, so multiply by each covariate's SD to compare.
R_SIGMA2, R_SCALE, R_ELEV_RAW, R_GRAD_RAW = 1.58, 48.0, 0.021, 5.846
sd_e = gr['cov'][:,0].std(); sd_g = gr['cov'][:,1].std()
print(f"{'':22s}{'field SD':>10s}{'elevation':>12s}{'gradient':>11s}")
print(f"{'from scratch':22s}{sd:10.2f}{be[1]:12.2f}{be[2]:11.2f}")
print(f"{'PyMC':22s}{tau_pm:10.2f}{bc_pm[0]:12.2f}{bc_pm[1]:11.2f}")
print(f"{'R (kppm, converted)':22s}{np.sqrt(R_SIGMA2):10.2f}{R_ELEV_RAW*sd_e:12.2f}{R_GRAD_RAW*sd_g:11.2f}")
print()
print('The three agree on the SIZE of the clustering to within a couple of percent, and disagree on the')
print('covariate coefficients by up to a factor of two. That is the honest summary of this model: how much')
print('clustering there is beyond terrain is well identified; how to split the signal between terrain and')
print('the latent field is not, because a smooth covariate and a smooth field are largely interchangeable.')
print()
print('Note also that R and Python disagree about what accounting for clustering even DOES. Here the')
print(f'coefficients MOVED -- elevation {b[1]:+.2f} under the plain Poisson against {be[1]:+.2f} under the LGCP -- while')
print('kppm returns coefficients identical to ppm and changes only their standard errors. Neither is wrong:')
print('kppm fixes the trend at the Poisson estimate and corrects the variance afterwards, whereas estimating')
print('beta and the field jointly lets the field re-allocate signal. Whether clustering changes your estimate')
print('or only your uncertainty is a property of the estimator, not of the trees.')
print()
print('One more thing worth flagging: the coefficient went UP. In the areal BYM models a latent field')
print('typically ABSORBS covariate signal and shrinks the coefficient. Terrain is itself spatially smooth, so')
print('field and covariate compete -- but the competition here resolved toward a STRONGER terrain association')
print('once the fine clumping had somewhere else to go. The direction is not predictable in advance.')

                        field SD   elevation   gradient
from scratch                1.28        0.41       0.47
PyMC                        1.29        0.47       0.51
R (kppm, converted)         1.26        0.17       0.33

The three agree on the SIZE of the clustering to within a couple of percent, and disagree on the
covariate coefficients by up to a factor of two. That is the honest summary of this model: how much
clustering there is beyond terrain is well identified; how to split the signal between terrain and
the latent field is not, because a smooth covariate and a smooth field are largely interchangeable.

Note also that R and Python disagree about what accounting for clustering even DOES. Here the
coefficients MOVED -- elevation +0.16 under the plain Poisson against +0.41 under the LGCP -- while
kppm returns coefficients identical to ppm and changes only their standard errors. Neither is wrong:
kppm fixes the trend at the Poisson estimate and corrects the variance afterwards, whereas estimating
beta and the field jointly lets the field re-allocate signal. Whether clustering changes your estimate
or only your uncertainty is a property of the estimator, not of the trees.

One more thing worth flagging: the coefficient went UP. In the areal BYM models a latent field
typically ABSORBS covariate signal and shrinks the coefficient. Terrain is itself spatially smooth, so
field and covariate compete -- but the competition here resolved toward a STRONGER terrain association
once the fine clumping had somewhere else to go. The direction is not predictable in advance.

gr2 = P.make_grid(pts, gcov['x'].to_numpy(), gcov['y'].to_numpy(),
                  gcov[['elev','grad']].to_numpy(), xr, yr, 12.5)
e2 = (gr2['cov'][:,0]-gr2['cov'][:,0].mean())/gr2['cov'][:,0].std()
g2 = (gr2['cov'][:,1]-gr2['cov'][:,1].mean())/gr2['cov'][:,1].std()
X2 = np.column_stack([np.ones(len(e2)), e2, g2])
res2 = P.lgcp_gibbs(gr2['counts'], X2, gr2['area'], gr2['W'].astype(float),
                    np.random.default_rng(5), draws=1000, burn=1000)
be2 = res2['beta'].mean(0); sd2 = res2['spatial_sd'].mean()
print(f"{'cell size':>12s}{'cells':>8s}{'field SD':>11s}{'elevation':>12s}{'gradient':>11s}")
print(f"{'25.0 m':>12s}{len(gr['counts']):8d}{sd:11.2f}{be[1]:12.2f}{be[2]:11.2f}")
print(f"{'12.5 m':>12s}{len(gr2['counts']):8d}{sd2:11.2f}{be2[1]:12.2f}{be2[2]:11.2f}")
print()
print(f'Quadrupling the number of cells leaves the clustering essentially where it was -- field SD {sd:.2f} -> {sd2:.2f}')
print(f'and gradient {be[2]:+.2f} -> {be2[2]:+.2f} -- while elevation moves {be[1]:+.2f} -> {be2[1]:+.2f}. The ICAR field is defined ON')
print('the grid, so exact invariance was never expected; what the check establishes is that the substantive')
print('conclusion does not rest on the discretisation. And it reproduces the same split a third time:')
print('the amount of clustering is stable across engines, priors and now resolution, while the share')
print('attributed to terrain rather than to the field is the quantity that keeps moving.')

   cell size   cells   field SD   elevation   gradient
      25.0 m     800       1.28        0.41       0.47
      12.5 m    3200       1.32        0.60       0.46

Quadrupling the number of cells leaves the clustering essentially where it was -- field SD 1.28 -> 1.32
and gradient +0.47 -> +0.46 -- while elevation moves +0.41 -> +0.60. The ICAR field is defined ON
the grid, so exact invariance was never expected; what the check establishes is that the substantive
conclusion does not rest on the discretisation. And it reproduces the same split a third time:
the amount of clustering is stable across engines, priors and now resolution, while the share
attributed to terrain rather than to the field is the quantity that keeps moving.

Spatial Point Processes — CSR, Intensity, and Cox Processes¶

Modelling where events happen¶

1. Where this is, the pattern, and a covariate¶

2. Is it random? — Ripley's L against CSR¶

3. Inhomogeneous Poisson — intensity from covariates¶

4. Log-Gaussian Cox process — clustering beyond covariates¶

5. Cross-check in PyMC¶

6. What is identified, and what is not¶

Is the grid fine enough?¶

7. Summary¶