I want to learn to use Mathematica to deduce abstract theoretical knowledge about spatial statistics. I found a demo file about spatial statistics in the 2015 Wolfram Technical Conference materials
https://library.wolfram.com/infocenter/Conferences/9273

However, because of the lack of an important part, it is currently unable to run. I hope someone can help rebuild this demo or add missing parts. My ordinary kriging code. I think there should be better practices in Mathematica, such as symbolic operations. Can anyone provide examples?

Remove["Global`*"];
(*Download air quality monitoring data*)
points = (ToExpression[#] + 0.) &@
       Cases[#, {x1_, x2_, x3_ /; NumberQ@ToExpression@x3}] &@
     Values@KeyTake[#, {"Longitude", "Latitude", "PM2.5"}] &@
   Import["https://opendata.epa.gov.tw/api/v1/AQI?%24skip=0&%24top=\
1000&%24format=json", {"Json", "Data"}];
(*define function*)
sphereModel = 
  Piecewise[{{c0, h == 0}, {c0 + c ((3 h)/(2 a) - h^3/(2 a^3)), 
     0 < h <= a}, {c0 + c, h > a}}];
model[x_] := sphereModel /. Flatten@{fit, h -> x};
(*semivariogram modeling*)
{semivariogramCloud, xMax, 
   yMax} = {#, Max@#[[All, 1]], Max@#[[All, 2]] + 1} &@
   ParallelMap[{EuclideanDistance[#[[1, 1 ;; 2]], #[[2, 1 ;; 2]]], 
      Variance[{#[[1, 3]], #[[2, 3]]}]} &, 
    DeleteDuplicates@Map[Sort, Permutations[points, {2}]]];
semivariogram = 
  Select[Length@# != 0 &]@
   Map[Map[Mean, Transpose[#]] &, 
    Map[#[[1]] &, 
     BinLists[semivariogramCloud, {0, xMax, 0.01}, {0, yMax, yMax}]]];
fit = FindFit[
   semivariogram, {sphereModel, {c0 >= 0, c >= 0, 
     0.005 <= a <= xMax}}, {c0, c, a}, h, 
   PerformanceGoal -> "Quality"];
(*Ordinary Kriging  & Cross Validation*)
crossValidation = Table[TakeDrop[points, {i}], {i, 1, Length@points}];
Do[
  {validationPoint, observePoint} = {#[[1, 1]], #[[2]]} &@
    crossValidation[[i]];
  observePoint = 
   Select[EuclideanDistance[validationPoint[[1 ;; 2]], #[[1 ;; 2]]] <=
        fit[[3, 2]] &]@observePoint;
  If[Length@observePoint != 0,
   n = Length@observePoint;
   kk = Table[1, n + 1, n + 1];         
   kk[[n + 1, n + 1]] = 0;        k = Table[1, n + 1];
   kk[[;; n, ;; n]] = DistanceMatrix@observePoint[[All, {1, 2}]];
   kk[[;; n, ;; n]] = Table[model@kk[[i, j]], {i, 1, n}, {j, 1, n}];
   k[[;; n]] = 
    Table[EuclideanDistance[validationPoint[[1 ;; 2]], 
      observePoint[[i, {1, 2}]]], {i, 1, n}];
   k[[;; n]] = Table[model@k[[i]], {i, 1, n}];
   w = Dot[Inverse@kk, k];
   prediction = Dot[w[[;; n]], observePoint[[All, 3]]];
   standErr = Sqrt@Dot[Transpose@w, k]; 
   Which[i == 1, 
    results = {{prediction, validationPoint[[3]], standErr}}, True, 
    results = 
     AppendTo[results, {prediction, validationPoint[[3]], standErr}]];
   ];
  , {i, 1, Length@crossValidation}];
(*performance*)
lm = LinearModelFit[results[[All, {1, 2}]], x, x];
lm["RSquared"]