Skip to content

Commit

Permalink
Update environment
Browse files Browse the repository at this point in the history
  • Loading branch information
juliohm committed Oct 21, 2024
1 parent fbaeaf0 commit 9c83ec6
Show file tree
Hide file tree
Showing 8 changed files with 304 additions and 254 deletions.
2 changes: 1 addition & 1 deletion 03-geoio.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ geotable = GeoIO.load("data/earth.tif")
```

If we save the geotable to a `.geojson` file on disk, and then load it back,
we observe that the `CartesianGrid` gets replaced by a `GeometrySet`:
we observe that the grid gets replaced by a geometry set:

```{julia}
fname = tempname() * ".geojson"
Expand Down
2 changes: 0 additions & 2 deletions 05-transforms.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -362,8 +362,6 @@ are the `Continuous` and the `Categorical` scientific traits. To convert (or coe
scientific traits of columns in a geotable, we can use the `Coerce` transform:

```{julia}
using GeoStats.DataScienceTraits: Continuous
st = georef((a=[1,2,2,2,3,3], b=[1,2,3,4,5,6])) |> Coerce("b" => Continuous)
```

Expand Down
27 changes: 15 additions & 12 deletions 10-correlation.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -180,12 +180,15 @@ n = 0.1
m = GaussianVariogram
g = m(range=r, sill=s, nugget=n)
xs = range(1.0e-6, 30.0, length=100)
ys = g.(xs)
fig = Mke.Figure()
Mke.Axis(fig[1,1], xlabel="h", ylabel="γ(h)", limits = (nothing, nothing, 0, 1.1))
Mke.plot!(g, maxlag=25, label="model")
Mke.vlines!(r, color="black", linestyle=:dash, label="range")
Mke.hlines!(s, color="teal", linestyle=:dash, label="sill")
Mke.hlines!(n, color="teal", linestyle=:dot, label="nugget")
ax = Mke.Axis(fig[1,1], xlabel="h", ylabel="γ(h)", limits = (nothing, nothing, 0, 1.1))
Mke.lines!(ax, xs, ys, color="slategray", label="model")
Mke.vlines!(ax, r, color="black", linestyle=:dash, label="range")
Mke.hlines!(ax, s, color="teal", linestyle=:dash, label="sill")
Mke.hlines!(ax, n, color="teal", linestyle=:dot, label="nugget")
Mke.axislegend("Elements", position = :rc)
fig
```
Expand Down Expand Up @@ -314,9 +317,9 @@ geospatial correlation. The most widely used are the `GaussianVariogram`, the
fig = Mke.Figure()
Mke.Axis(fig[1,1])
Mke.plot!(γ1, maxlag=2, vcolor = "teal", label = "Gaussian")
Mke.plot!(γ2, maxlag=2, vcolor = "slategray3", label = "Spherical")
Mke.plot!(γ3, maxlag=2, vcolor = "brown", label = "Exponential")
varioplot!(γ1, maxlag=2.0, color = "teal", label = "Gaussian")
varioplot!(γ2, maxlag=2.0, color = "slategray3", label = "Spherical")
varioplot!(γ3, maxlag=2.0, color = "brown", label = "Exponential")
Mke.axislegend("Model", position = :rb)
fig
```
Expand Down Expand Up @@ -406,7 +409,7 @@ g = EmpiricalVariogram(img, :Z, maxlag = 50.0)
```

```{julia}
Mke.plot(g)
varioplot(g)
```

::: {.callout-note}
Expand All @@ -424,8 +427,8 @@ specific directions:
gₕ = DirectionalVariogram((1.0, 0.0), img, :Z, maxlag = 50.0)
gᵥ = DirectionalVariogram((0.0, 1.0), img, :Z, maxlag = 50.0)
Mke.plot(gₕ, hshow = false, vcolor = "maroon")
Mke.plot!(gᵥ, hshow = false, vcolor = "slategray")
varioplot(gₕ, showhist = false, color = "maroon")
varioplot!(gᵥ, showhist = false, color = "slategray")
Mke.current_figure()
```

Expand All @@ -447,7 +450,7 @@ solid curve over the heatmap:
```{julia}
fig = Mke.Figure()
ax = Mke.PolarAxis(fig[1,1])
Mke.plot!(ax, gₚ)
varioplot!(ax, gₚ)
fig
```

Expand Down
6 changes: 3 additions & 3 deletions 11-interpolation.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -220,7 +220,7 @@ g = EmpiricalVariogram(samples, "Z", maxlag = 100.0)
```

```{julia}
Mke.plot(g)
varioplot(g)
```

A better alternative in this case is to use the robust `:cressie`
Expand All @@ -231,7 +231,7 @@ g = EmpiricalVariogram(samples, "Z", maxlag = 100.0, estimator = :cressie)
```

```{julia}
Mke.plot(g)
varioplot(g)
```

After estimating the empirical variogram, the next step consists of fitting
Expand All @@ -242,7 +242,7 @@ a theoretical model. The behavior near the origin resembles a `SphericalVariogra
```

```{julia}
Mke.plot(γ, maxlag = 100.0)
varioplot(γ, maxlag = 100.0)
```

Now that we extracted the geospatial correlation from the samples, we can
Expand Down
22 changes: 9 additions & 13 deletions 12-mining.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -155,18 +155,11 @@ elements, we visualize the `values` of the drill hole samples with the
`pairplot`:

```{julia}
dtable |> Select("Cu", "Au", "Ag", "S") |> DropUnits() |> values |> pairplot
dtable |> Select("Cu", "Au", "Ag", "S") |> values |> pairplot
```

We can observe that the distribution is very skewed.

::: {.callout-note}

The `DropUnits` transform can be useful to drop units from the columns of a
table before calling functions that do not support units yet (e.g., `pairplot`).

:::

### Domain of interpolation

Before we can interpolate these variables, we need to define our domain of interpolation.
Expand Down Expand Up @@ -274,7 +267,7 @@ In order to remove compositional data constraints, we will perform the centered
transform (`CLR`) from the [CoDa.jl](https://github.com/JuliaEarth/CoDa.jl) module:

```{julia}
grades = dtable |> Select("Cu", "Au", "Ag", "S") |> DropUnits()
grades = dtable |> Select("Cu", "Au", "Ag", "S")
grades |> CLR() |> values |> pairplot
```
Expand Down Expand Up @@ -306,12 +299,15 @@ samples

#### Geospatial correlation

Let's fit a theoretical variogram for all four (independent) variables:
Let's fit a theoretical variogram for all four (independent) variables
up to a given maximum lag:

```{julia}
maxlag = 300.0u"m"
ns = setdiff(names(samples), ["geometry"])
gs = [EmpiricalVariogram(samples, n, estimator = :cressie) for n in ns]
gs = [EmpiricalVariogram(samples, n, maxlag = maxlag) for n in ns]
γs = [GeoStatsFunctions.fit(Variogram, g, h -> 1 / h^2) for g in gs]
```
Expand All @@ -327,8 +323,8 @@ model. The constant `100` was chosen based on visual inspection of the

```{julia}
function gammaplot(n, g, γ)
Mke.plot(g, axis = (; title = n))
Mke.plot!(γ, maxlag = 300, vcolor = "teal")
varioplot(g, axis = (; title = n))
varioplot!(γ, maxlag = maxlag, color = "teal")
Mke.current_figure()
end
Expand Down
Loading

0 comments on commit 9c83ec6

Please sign in to comment.