# UniversalKriging, gauged precipiation interpolation with dem as a covariate

**URL:** <https://discourse.julialang.org/t/universalkriging-gauged-precipiation-interpolation-with-dem-as-a-covariate/130519>\
**Category:** General Usage\
**Tags:** question, kriging\
**Created:** [July 6, 2025, 8:37am UTC](https://discourse.julialang.org/t/universalkriging-gauged-precipiation-interpolation-with-dem-as-a-covariate/130519 "2025-07-06T08:37:49Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![kongdd](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kongdd/32/17058_2.png) [@kongdd](https://discourse.julialang.org/u/kongdd)\
**Post date:** [July 6, 2025, 8:37am UTC](https://discourse.julialang.org/t/universalkriging-gauged-precipiation-interpolation-with-dem-as-a-covariate/130519/1 "2025-07-06T08:37:49Z")

</div>

The following script is the UniversalKriging script gived by grok 3.5, which is trying to interpolate  
gauged precipiation with dem as a covariate. Any idea how to achieve this in Julia?

```python
import numpy as np
from pykrige.uk import UniversalKriging

# Example data: x, y (coordinates), z (precipitation), h (elevation)
x = np.array([100, 200, 300, 400, 500]) # Station x coordinates
y = np.array([50, 150, 250, 350, 450]) # Station y coordinates
z = np.array([20, 25, 30, 35, 40]) # Precipitation (mm)
h = np.array([500, 600, 700, 800, 900]) # Elevation (m)

# Target grid points
gridx = np.arange(100, 500, 10)
gridy = np.arange(50, 450, 10)
grid_h = np.random.rand(len(gridx), len(gridy)) * 1000 # Assumed elevation grid

# Universal Kriging, specifying elevation as external drift
UK = UniversalKriging(
    x, y, z,
    variogram_model='spherical', # Variogram model
    drift_terms=['external_drift'], # Use elevation as external drift
    external_drift=h, # Station elevation
    external_drift_x=x, external_drift_y=y # Coordinates corresponding to elevation
)

# Perform interpolation
z_pred, variance = UK.execute('grid', gridx, gridy, external_drift=grid_h)

# Visualization
import matplotlib.pyplot as plt
plt.contourf(gridx, gridy, z_pred, cmap='viridis')
plt.colorbar(label='Precipitation (mm)')
plt.scatter(x, y, c='red', marker='o', label='Stations')
plt.title('Precipitation Interpolation (Considering Elevation)')
plt.legend()
plt.show()

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [July 6, 2025, 6:54pm UTC](https://discourse.julialang.org/t/universalkriging-gauged-precipiation-interpolation-with-dem-as-a-covariate/130519/2 "2025-07-06T18:54:16Z")

</div>

@kongdd I recommend reading the GeoStats.jl documentation and examples related to Kriging. If something is not clear, please reach out.