<jupyter_start><jupyter_text>Poisson Kriging - Area to Area Kriging Table of Contents:1. Load areal and point data,2. Load semivariogram (regularized),3. Remove 10% of areal dataset,4. Predict values at unknown locations,5. Analyse Forecast Bias and Root Mean Squared Error of prediction. Level: Advanced Changelog| Date | Change description | Author ||------|--------------------|--------|| 2021-05-28 | Updated paths for input/output data | @szymon-datalions || 2021-05-11 | Refactored TheoreticalSemivariogram class | @szymon-datalions || 2021-03-31 | Update related to the change of semivariogram weighting. Updated cancer rates data. | @szymon-datalions | IntroductionTo start this tutorial it is required to understand concepts in the **Ordinary and Simple Kriging** and **Semivariogram Regularization** tutorials.Poisson Kriging technique is used to model spatial count data. We are analyzing special case where data is counted over areas. Those areas may have irregular shapes and sizes beacuse they represent administrative regions. In this tutorial we try to predict rates of Breast Cancer in Pennsylvania counties. Along with areal data we use U.S. Census 2010 data for population blocks.> Breast cancer rates data is stored in the shapefile in folder `sample_data/areal_data/cancer_data.shp`.> Population blocks data is stored in the shapefile in folder `sample_data/population_data/cancer_population_base.shp`In tutorial **Poisson Kriging - Centroid based approach** we've made assumption that all areas have the same size and shape. This wasn't true. Now we use Area-to-Area Kriging technique and predict and evaluate missing values.Area-to-Area and Area-to-Point Poisson Kriging are slower than simplified Kriging with Centroids, but gives more reliable results because they are tuned to areal size and shape.This tutorial covers following steps:1. Read and explore data,2. Load semivariogram model,3. Prepare training and test data,4. Predict values of unknown locations and calculate forecast bias and root mean squared error,5. Analyse error metrics.  1) Read and explore data<jupyter_code>import numpy as np
import pandas as pd
import geopandas as gpd

from pyinterpolate.io_ops import prepare_areal_shapefile, get_points_within_area
from pyinterpolate.semivariance import TheoreticalSemivariogram
from pyinterpolate.kriging import ArealKriging
areal_data_file = '../sample_data/areal_data/cancer_data.shp'
point_data_file = '../sample_data/population_data/cancer_population_base.shp'

areal_id_column_name = 'FIPS'
areal_val_column_name = 'rate'
points_val_column_name = 'POP10'


areal_data = prepare_areal_shapefile(areal_data_file,
                                     id_column_name=areal_id_column_name,
                                     value_column_name=areal_val_column_name)
point_data = get_points_within_area(areal_data_file, point_data_file, areal_id_col_name=areal_id_column_name,
                                   points_val_col_name=points_val_column_name)
# Lets take a look into a map of areal counts

gdf = gpd.read_file(areal_data_file)
gdf.plot(column='rate', cmap='Spectral_r', legend=True)<jupyter_output><empty_output><jupyter_text>Clarification: It is good idea to look into the spatial patterns in dataset and to visually check if our data do not have any NaN values. We use geopandas GeoDataFrame plot function for it with color map which is diverging color classes well.  2) Load semivariogram model In this step we load regularized semivariogram from the tutorial **Semivariogram Regularization (Intermediate)**. You can always perform semivariogram regulariztion along with the Poisson Kriging but it is a very long process and it is more convenient to separate those two steps.<jupyter_code>semivariogram = TheoreticalSemivariogram()  # Create TheoreticalSemivariogram object
semivariogram.import_model('output/regularized_model.csv')  # Load regularized semivariogram<jupyter_output><empty_output><jupyter_text>3) Prepare training and test dataWe simply remove 10% of random ID's from areal dataset (the same for points) and create four arrays: two training arrays with areal and point geometry and values and two test arrays with the same categories of data.<jupyter_code># Remove 10% of rows (values) to test our model

def create_test_areal_set(areal_dataset, points_dataset, frac=0.1):
    """
    Function creates set with random areas for tests
    
    INPUT:
    :param areal_dataset: (numpy array) array with rows of areal records,
    :param points_dataset: (numpy array) array with rows of points records,
    :param frac: (float) number of elements in a test set
    
    OUTPUT:
    return: [training_set_areal(numpy array),
             test_set_areal (numpy array),
             training_set_points (numpy array),
             test_set_points (numpy array)]
    """
    
    all_ids = areal_dataset[:, 0]
    training_set_size = int(len(all_ids) * (1 - frac))
    training_ids = np.random.choice(all_ids,
                                    size=training_set_size,
                                    replace=False)
    training_areas = np.array(
            [a for a in areal_dataset if a[0] in training_ids])
    test_areas = np.array(
            [a for a in areal_dataset if a[0] not in training_ids])
    training_points = np.array(
            [pt for pt in points_dataset if pt[0] in training_ids])
    test_points = np.array(
            [pt for pt in points_dataset if pt[0] not in training_ids])
    

    output = [training_areas, test_areas, training_points, test_points]
    return output

ar_train, ar_test, pt_train, pt_test = create_test_areal_set(areal_data, point_data)<jupyter_output><empty_output><jupyter_text>4) Predict values at unknown locations and calculate forecast bias and root mean squared errorDo you know scikit-learn's fit-predict-transform estimation procedure? With Kriging it is more complicated because we are fitting semivariogram model to data and we pass this fitted model into kriging estimator. That's why package with *TheoreticalSemivariogram()* class is named **fit_semivariance**.You may start work with predictions with fitted semivariogram model. Area-to-Area Poisson Kriging takes three arguments during initialization:=> semivariogram model (fitted semivariogram model),=> known areas (training set),=> known points (training set).After initilization you may perform prediction. This process requires you to control two hyperparameters:=> **number of observations** (the most important parameter - how many neighbors are affecting your area of analysis),=> minimum **search radius** (minimum search radius for analysis from your area of interest, if there are less areas than number of observations then the new, next closest neighbors are included in the analysis),Why **radius** parameter is less important than **number of observations**? It is a tradeoff between mathemtical accuracy and reality. Usually your areas of interest have different shapes and sizes and with very large county and small search radius from the centroid you may not find any neighbour!**The main difference between Area-to-Area and Centroid Poisson Kriging is that for the former you pass only variable with smaller blocks which are representing fine scale process withing area with missing value. Centroid PK is "feed" with missing area geometry too.**Algorithm in the cell below iteratively picks one area from the test set and perform prediction. Then forecast bias which is difference between real and predicted value is calculated. Forecast Bias is useful to understand if our algorithm predictions are too low or too high (under- or overestimation). The next parameter is Root Mean Squared Error. This kind of measure tells us more about outliers and very large differences between predictions and real values. We will see it in the last part of tutorial.Your work with Poisson Kriging (or Kriging) will be usually the same:- prepare training and test data,- use training data to prepare semivariogram model,- test different hyperparameters with a test set to find optimal number of neighbors which are affecting your area of analysis,- forecast into unknown areas OR perform data smoothing.<jupyter_code>number_of_obs = 8
radius = 60000  # 60 km

kriging_model = ArealKriging(semivariogram_model=semivariogram,
    known_areas=ar_train,
    known_areas_points=pt_train,
    kriging_type='ata')

kriged_predictions = []
for unknown_area in ar_test:
    unknown_points = pt_test[pt_test[:, 0] == unknown_area[0]][0]

    # Predict
    try:
        predicted = kriging_model.predict(unknown_points,
                                          number_of_obs,
                                          radius)
        predicted_rates = int(predicted[0])
        fb = (unknown_area[-1] - predicted_rates)
        rmse = np.sqrt((fb)**2)
        
        kriged_predictions.append([unknown_area[0], predicted_rates, unknown_area[-1], fb, rmse])
    except ValueError:
        fb = np.nan
        rmse = np.nan
        kriged_predictions.append([unknown_area[0], np.nan, unknown_area[-1], fb, rmse])  # All NaN's
kriged_predictions = np.array(kriged_predictions)
pred_df = pd.DataFrame(data=kriged_predictions[:, 1:],
                       index=kriged_predictions[:, 0],
                       columns=['predicted', 'real', 'forecast bias', 'rmse'])  # Store results in DataFrame<jupyter_output><empty_output><jupyter_text>5) Analyze Forecast Bias and Root Mean Squared Error of predictionThe last step is analysis of errors. We plot two histograms: of forecast bias and of root mean squared error then we calculate base statistics of a dataset.<jupyter_code># Show histograms of errors

pred_df['forecast bias'].plot.hist(bins=10)
pred_df['rmse'].plot.hist(bins=10)
pred_df.describe()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06654437990420871, 'fraction_numerical': 0.008458167736675837, 'mean_word_length': 4.138219895287958, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 28, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '256546', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2767, 'n_tokens_neox': 2623, 'n_words': 1074, 'path': 'tutorials/Poisson Kriging - Area to Area (Advanced).ipynb', 'repo_name': 'scottgallacher-3/pyinterpolate', 'size': '9813'}