Conditioning with External Drift Kriging using GSTools

[LorenzEh]

Hi, I want to impute missing data on abortion rates of US-Counties. I´m comparing different Methods, one of them is Kriging. I used Pykrige for External Dirft Kriging:

def Kriging_with_SGS(df, variable, kriged_var, variogram, lags, drifts=None):
    interpolation_points = df[~df[variable].isnull()]
    interpolation_points.reset_index(drop=True, inplace=True)

    known_x = interpolation_points.geometry.centroid.x.values
    known_y = interpolation_points.geometry.centroid.y.values
    values = interpolation_points[variable].values
    external_drifts = []

    if drifts is not None:
        
        for drift in drifts:
            if drift in interpolation_points.columns:
                external_drift = interpolation_points[drift].values
                external_drifts.append(external_drift)
            else:
                print(f"The drift variable '{drift}' is missing.")

    if drifts:
        uk_all = UniversalKriging(
            known_x,
            known_y,
            values,
            variogram_model=variogram,
            drift_terms=["point_log"],
            point_drift=external_drifts if external_drifts else None,
            enable_plotting=True,
            nlags=lags,
            verbose=True
        )
        
    else:
        uk_all = UniversalKriging(
            known_x,
            known_y,
            values,
            variogram_model=variogram,
            enable_plotting=True,
            nlags=lags,
            verbose=True
        )

    df[kriged_var] = uk_all.execute(
        'points',
        df.geometry.centroid.x.values,
        df.geometry.centroid.y.values
    )

With this code I managed to get the kriging results for the Counties which have missing values in the target variable. But I have the feeling, that I´m loosing to much variance in the far field. Therefore I want to include conditioning in my method:

def cond_krig(df, variable, iterations, drift):
    # create a df with only those values which exist:
    interpolation_points = df[~df[variable].isnull()]
    interpolation_points.reset_index(drop=True, inplace=True)

    known_x = interpolation_points.geometry.centroid.x.values
    known_y = interpolation_points.geometry.centroid.y.values
    
    cond_val = interpolation_points[variable].values
    
    # define the grid with all points that I have
    all_x = df.geometry.centroid.x.values 
    all_y = df.geometry.centroid.y.values

    
    ext_drift = interpolation_points[drift].values 
    grid_drift = df[drift].values
    
    model = gs.Matern(dim=2, 
                      var=0.09, 
                      len_scale=11284.23, 
                      nugget=0, 
                      nu=0.20)  

    krig = gs.krige.ExtDrift(model = model, 
                             cond_pos = (known_x, known_y), 
                             cond_val = cond_val, 
                             ext_drift = ext_drift)
    
    krig.unstructured([all_x, all_y], ext_drift=grid_drift)
    
    cond_srf = gs.CondSRF(krig)  
    
    cond_srf.set_pos([all_x, all_y], 
                     "unstructured")

    seed = gs.random.MasterRNG(20170519)
    
    for i in range(iterations):
        cond_srf(seed=seed(), store=[f"f{i}", False, False])
        print("___")
        print(i)
        print("___")
        
    fields = [cond_srf[f"f{i}"] for i in range(iterations)]
    
    filtered_df["ensemble mean"] = np.mean(fields, axis=0)
    filtered_df["kriged field"] = cond_srf.krig.field 

I´m getting: Value Error:Krige: wrong number of ext. drift values. Thus I have two questions:

1. What am I doing wrong in the second code, I provided? I know that I can make some fine-tuning with latlon and other parameters but there seems to be a substential error in my code.
2. Is it possible to use the kriging model from Pykrige for conditioning? Something along these lines:

krig = uk_all.execute(
    'points',
    df.geometry.centroid.x.values,
    df.geometry.centroid.y.values
)[0]

    
cond_srf = gs.CondSRF(krig)  
cond_srf.set_pos([df.geometry.centroid.x.values, 
                  df.geometry.centroid.y.values], 
                 "unstructured")

I´m happy about any answer or advice. Many thanks!

[LSchueler]

Hi LorenzEh,

sorry for the late reply. As I can't reproduce your code without the data, my first guess would be that the shape of the ext_drift is wrong. Can you give me the shapes of your input data, i.e. ext_drift, grid_drift, and known_x and known_y?

And in which line is the error raised?

[LorenzEh]

Thanks for your answer. The external drift values for the whole grid (grid_drift) don´t contain nan values and differ from the condition values. I´m planning to use "population of women age 15-44" and "distance to abortion clinic" as external dirfts. Both of them highly correlate with abortion rates, and mapping them reveals similar spatial patterns (eg. counties with high abortion rates have high populations of women between 15 and 44).

I replicated your code with my data and it works perfectly fine. But I´m still struggeling to combine external Kriging with conditioning. I´m getting ValueError: Krige: wrong number of ext. drift values. when trying to loop:

import numpy as np
import matplotlib.pyplot as pt
import gstools as gs

X = 100_000
Y = 100_000
N = 100
N_ALL = 1000
VAL_MAX = 20

rng = np.random.default_rng()

known_x = rng.uniform(0, X, N)
known_y = rng.uniform(0, Y, N)

all_x = rng.uniform(0, X, N_ALL)
all_y = rng.uniform(0, Y, N_ALL)

cond_val = rng.uniform(0., VAL_MAX, size=N)

# some external information/ function, differing from your `cond_val`s
ext_drift = np.linspace(0., VAL_MAX, N)**0.5
ext_drift_all = np.linspace(0., VAL_MAX, N_ALL)**0.5

model = gs.Matern(
    dim=2,
    var=0.09,
    len_scale=11284.23,
    nugget=0.,
    nu=0.2,
)

krige = gs.krige.ExtDrift(
    model=model,
    cond_pos=(known_x, known_y),
    cond_val=cond_val,
    ext_drift=ext_drift,
)
krige((all_x, all_y), ext_drift=ext_drift_all)

cond_srf = gs.CondSRF(krige)  

cond_srf.set_pos((all_x, all_y), 
                 "unstructured")

seed = gs.random.MasterRNG(20170519)
fields = {} 

# here is the problem: 
for i in range(5):
    cond_srf(seed=seed(), store=[f"f{i}"])
    fields[f"field_{i}"] = cond_srf[f"f{i}"]
    print("___")
    print(i)
    print("___")
    
fields = [cond_srf[f"f{i}"] for i in range(5)]

mean = np.mean(fields, axis=0)
extra_variance = np.var(fields, axis=0) 
krig_and_variance = cond_srf.krige.field + extra_variance

I did a similar thing using Ordinary Kriging, without getting an error:

def geokrig(df, variable, iterations, model): # das funktioniert und schaut gut aus! chunk_size as a __call__ function  unbedingt auschecken!
    
    intpoints = df[~df[variable].isnull()]
    intpoints.reset_index(drop=True, inplace=True)
    lat = intpoints["geometry"].geometry.centroid.x
    lon = intpoints["geometry"].geometry.centroid.y  
    int_data = intpoints[variable].values

    normalizer = gs.normalizer.Manly(int_data) 

    bin_center, vario = gs.vario_estimate(
                        (lat, lon), 
                        int_data, 
                        latlon=True, 
                        mesh_type="unstructured",
                        geo_scale=gs.DEGREE_SCALE,
                        normalizer=normalizer)
                        
    print("estimated bin number:", len(bin_center))
    print("maximal bin distance:", max(bin_center))
    
    # fit a Matern model to the variogram
    model = model(dim = 2, latlon=True, geo_scale=gs.DEGREE_SCALE)
    model.fit_variogram(bin_center, vario)  
    
    # plotting:
    ax = model.plot("vario_yadrenko", x_max=max(bin_center))
    ax.scatter(bin_center, vario)
    print(model)
    
    # kriging:
    krig = gs.krige.Ordinary(model = model, 
                             cond_pos = (lat, lon), 
                             cond_val = int_data)
                             # normalizer = normalizer) 

    g_lat = df["geometry"].geometry.centroid.x
    g_lon = df["geometry"].geometry.centroid.y
    
    krig.set_pos((g_lat, g_lon), mesh_type="unstructured") 
    krig(return_var=True, store="krig") 
    krig(only_mean=True, store="mean_field")
    # krig(mesh_type = "unstrctured") 
    
    cond_srf = gs.CondSRF(krige = krig, 
                          latlon = True, 
                          normalizer = normalizer) 
    
    cond_srf.set_pos(pos = (g_lat, g_lon), 
                     mesh_type = "unstructured")
    
    #cond_srf.transform("lognormal")
    seed = gs.random.MasterRNG(20170519)
    fields = {} 
    
    for i in range(iterations):
        cond_srf(seed=seed(), store=[f"f{i}", False, False])
        fields[f"field_{i}"] = cond_srf[f"f{i}"]
        print("_____________")
        print("Interation:", i)
        
    fields = [cond_srf[f"f{i}"] for i in range(iterations)]
    
          
    df["cond_krig"] = krig["krig"] 
    extra_variance = np.var(fields, axis=0) 
    df["cond_krig_variance"] = cond_srf.krige.field + extra_variance
    df["ensemble_mean"] = np.mean(fields, axis=0) 
    # df["kriged_values_check"] = cond_srf.krige.field 
    # filtered_df["conditional ensamble mean"] = cond_srf.mean()
    
    fields_df = pd.DataFrame(fields)
    fields_df = fields_df.T.reset_index(drop=True)
    
    return fields_df

Many Thanks! P.S. including cost of contraceptions seems like a great idea, so thanks for that input as well! ;^)

[LSchueler]

Ah, okay, I somehow assumed that the drift argument in your functions only indexes a single column, sorry for that. And I'm sorry for assuming that you didn't really understand external drift kriging.

I have the strong feeling that you found a bug in GSTools and that we never actually checked if the conditioning also works with external drift kriging. But I have to dig into that a little bit further. I'll come back to you soon.

[MuellerSeb]

Hey there,

jumping in here to point out, that all key-word arguments in CondSRF.__call__ are forwarded to the underlying kriging class.
Since you are using external drift kriging, you need to specify the external drift when calling the conditioned field class:

cond_srf((all_x, all_y), ext_drift=ext_drift_all)

That should work and explains the error message (which could be better in this case of course...).

Since we have the set_pos method, it would be a nice feature to also be able to set the external drift there, since it only depends on the given positions.

Hope that helps.

Cheers, Sebastian

[LSchueler]

Nice, thanks! I can confirm that if you change the line in the for loop cond_srf(seed=seed(), store=[f"f{i}"]) to cond_srf(seed=seed(), store=[f"f{i}"], ext_drift=ext_drift_all) (or grid_drift in your original code), it works.

[LorenzEh]

Yes, it works now! Many thanks for your answers & the great tool you provide!

Best whishes,
Lorenz