<jupyter_start><jupyter_text>Data reclassificationReclassifying data based on specific criteria is a common task when doing GIS analysis. The purpose of this lesson is to see how we can reclassify values based on some criteria which can be whatever, such as:```1. if travel time to my work is less than 30 minutes    AND    2. the rent of the apartment is less than 1000 € per month    ------------------------------------------------------    IF TRUE: ==> I go to view it and try to rent the apartment    IF NOT TRUE: ==> I continue looking for something else```In this tutorial, we will:1. Use classification schemes from the PySAL [mapclassify library](https://pysal.org/mapclassify/) to classify travel times into multiple classes.2. Create a custom classifier to classify travel times and distances in order to find out good locations to buy an apartment with these conditions:   - good public transport accessibility to city center   - bit further away from city center where the prices are presumably lower Input dataWe will use [Travel Time Matrix data from Helsinki](https://blogs.helsinki.fi/accessibility/helsinki-region-travel-time-matrix/) that contains travel time and distance information for routes between all 250 m x 250 m grid cell centroids (n = 13231) in the Capital Region of Helsinki by walking, cycling, public transportation and car.In this tutorial, we will use the geojson file generated in the previous section:`"data/TravelTimes_to_5975375_RailwayStation_Helsinki.geojson"`Alternatively, you can re-download [L4 data](https://github.com/AutoGIS/data/raw/master/L4_data.zip) and use `"data/Travel_times_to_5975375_RailwayStation.shp"` as input file in here. Common classifiers Classification schemes for thematic maps[PySAL](http://pysal.readthedocs.io/en/latest) -module is an extensive Python library for spatial analysis. It also includes all of the most common data classifiers that are used commonly e.g. when visualizing data. Available map classifiers in [pysal's mapclassify -module](https://pysal.readthedocs.io/en/v1.11.0/library/esda/mapclassify.html): - Box_Plot - Equal_Interval - Fisher_Jenks - Fisher_Jenks_Sampled - HeadTail_Breaks - Jenks_Caspall - Jenks_Caspall_Forced - Jenks_Caspall_Sampled - Max_P_Classifier - Maximum_Breaks - Natural_Breaks - Quantiles - Percentiles - Std_Mean - User_Defined- First, we need to read our Travel Time data from Helsinki:<jupyter_code>import geopandas as gpd

fp = "data/TravelTimes_to_5975375_RailwayStation_Helsinki.geojson"

# Read the GeoJSON file similarly as Shapefile
acc = gpd.read_file(fp)

# Let's see what we have
print(acc.head(2))<jupyter_output>car_m_d  car_m_t  car_r_d  car_r_t  from_id  pt_m_d  pt_m_t  pt_m_tt  \
0    29476       41    29483       46  5876274   29990      76       95   
1    29456       41    29462       46  5876275   29866      74       95   

   pt_r_d  pt_r_t  pt_r_tt    to_id  walk_d  walk_t    GML_ID   NAMEFIN  \
0   24984      77       99  5975375   25532     365  27517366  Helsinki   
1   24860      75       93  5975375   25408     363  27517366  Helsinki   

       NAMESWE NATCODE                                           geometry  
0  Helsingfors     091  POLYGON ((402250.000 6685750.000, 402024.224 6...  
1  Helsingfors     091  POLYGON ((402367.890 6685750.000, 402250.000 6...<jupyter_text>As we can see, there are plenty of different variables (see [from here the description](http://blogs.helsinki.fi/accessibility/helsinki-region-travel-time-matrix-2015) for all attributes) but what we are interested in are columns called `pt_r_tt` which is telling the time in minutes that it takes to reach city center from different parts of the city, and `walk_d` that tells the network distance by roads to reach city center from different parts of the city (almost equal to Euclidian distance).**The NoData values are presented with value -1**. - Thus we need to remove the No Data values first.<jupyter_code># Include only data that is above or equal to 0
acc = acc.loc[acc['pt_r_tt'] >=0]<jupyter_output><empty_output><jupyter_text>- Let's plot the data and see how it looks like- `cmap` parameter defines the color map. Read more about [choosing colormaps in matplotlib](https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html)- `scheme` option scales the colors according to a classification scheme (requires `mapclassify` module to be installed):<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

# Plot using 9 classes and classify the values using "Natural Breaks" classification
acc.plot(column="pt_r_tt", scheme="Natural_Breaks", k=9, cmap="RdYlBu", linewidth=0, legend=True)

# Use tight layout
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>As we can see from this map, the travel times are lower in the south where the city center is located but there are some areas of "good" accessibility also in some other areas (where the color is red).- Let's also make a plot about walking distances:<jupyter_code># Plot walking distance
acc.plot(column="walk_d", scheme="Natural_Breaks", k=9, cmap="RdYlBu", linewidth=0, legend=True)

# Use tight layour
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Okay, from here we can see that the walking distances (along road network) reminds more or less Euclidian distances.  Applying classifiers to dataAs mentioned, the `scheme` option defines the classification scheme using `pysal/mapclassify`. Let's have a closer look at how these classifiers work.<jupyter_code>import mapclassify<jupyter_output><empty_output><jupyter_text>- Natural Breaks<jupyter_code>mapclassify.NaturalBreaks(y=acc['pt_r_tt'], k=9)<jupyter_output><empty_output><jupyter_text>- Quantiles (default is 5 classes):<jupyter_code>mapclassify.Quantiles(y=acc['pt_r_tt'])<jupyter_output><empty_output><jupyter_text>- It's possible to extract the threshold values into an array:<jupyter_code>classifier = mapclassify.NaturalBreaks(y=acc['pt_r_tt'], k=9)
classifier.bins<jupyter_output><empty_output><jupyter_text>- Let's apply one of the `Pysal` classifiers into our data and classify the travel times by public transport into 9 classes- The classifier needs to be initialized first with `make()` function that takes the number of desired classes as input parameter<jupyter_code># Create a Natural Breaks classifier
classifier = mapclassify.NaturalBreaks.make(k=9)<jupyter_output><empty_output><jupyter_text>- Now we can apply that classifier into our data by using `apply` -function<jupyter_code># Classify the data
classifications = acc[['pt_r_tt']].apply(classifier)

# Let's see what we have
classifications.head()
type(classifications)<jupyter_output><empty_output><jupyter_text>Okay, so now we have a DataFrame where our input column was classified into 9 different classes (numbers 1-9) based on [Natural Breaks classification](http://wiki-1-1930356585.us-east-1.elb.amazonaws.com/wiki/index.php/Jenks_Natural_Breaks_Classification).- We can also add the classification values directly into a new column in our dataframe:<jupyter_code># Rename the column so that we know that it was classified with natural breaks
acc['nb_pt_r_tt'] = acc[['pt_r_tt']].apply(classifier)

# Check the original values and classification
acc[['pt_r_tt', 'nb_pt_r_tt']].head()<jupyter_output><empty_output><jupyter_text>Great, now we have those values in our accessibility GeoDataFrame. Let's visualize the results and see how they look.<jupyter_code># Plot
acc.plot(column="nb_pt_r_tt", linewidth=0, legend=True)

# Use tight layout
plt.tight_layout()
# plot them side-by-side
%matplotlib inline
import matplotlib.pyplot as plt

#basic config
fig, (ax1,ax2) = plt.subplots(nrows=1, ncols=2, figsize=(20, 16))
#ax1, ax2 = axes

#1st plot with 9 classes
ax1 = acc.plot(ax=ax1, column="nb_pt_r_tt", linewidth=0, legend=True)
ax1.set_title('Natural Breaks')

#2nd plot with original data
ax2 = acc.plot(ax=ax2, column="pt_r_tt", linewidth=0, legend=True)
ax2.set_title('Original')

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>And here we go, now we have a map where we have used one of the common classifiers to classify our data into 9 classes.  Plotting a histogramA histogram is a graphic representation of the distribution of the data. When classifying the data, it's always good to consider how the data is distributed, and how the classification shceme divides values into different ranges. - plot the histogram using [pandas.DataFrame.plot.hist](https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.plot.hist.html)- Number of histogram bins (groups of data) can be controlled using the parameter `bins`:<jupyter_code># Histogram for public transport rush hour travel time
acc['pt_r_tt'].plot.hist(bins=50)<jupyter_output><empty_output><jupyter_text>Let's also add threshold values on thop of the histogram as vertical lines.- Natural Breaks:<jupyter_code># Define classifier
classifier = mapclassify.NaturalBreaks(y=acc['pt_r_tt'], k=9)

# Plot histogram for public transport rush hour travel time
acc['pt_r_tt'].plot.hist(bins=50)

# Add vertical lines for class breaks
for value in classifier.bins:
    plt.axvline(value, color='k', linestyle='dashed', linewidth=1)<jupyter_output><empty_output><jupyter_text>- Quantiles:<jupyter_code># Define classifier
classifier = mapclassify.Quantiles(y=acc['pt_r_tt'])

# Plot histogram for public transport rush hour travel time
acc['pt_r_tt'].plot.hist(bins=50)

for value in classifier.bins:
    plt.axvline(value, color='k', linestyle='dashed', linewidth=1)<jupyter_output><empty_output><jupyter_text>**Task**Select another column from the data (for example, travel times by car: `car_r_t`). Do the following visualizations using one of the classification schemes available from [pysal/mapclassify](https://github.com/pysal/mapclassify):    - histogram with vertical lines showing the classification bins, `Equal_Interval`, `Fisher_Jenks`,     `Fisher_Jenks_Sampled`, `HeadTail_Breaks`, `Jenks_Caspall`, `Jenks_Caspall_Forced`, `Jenks_Caspall_Sampled`, `Max_P_Classifier`, `Maximum_Breaks`, `Natural_Breaks`, `Quantiles`, `Percentiles`, `Std_Mean`.- thematic map using the classification scheme<jupyter_code>import random

def classify_and_plot(gdf, classification = mapclassify.NaturalBreaks):
    """
    A function that accepts a geodataframe; chooses a column and generates histogram and thematic map plots based on a 
    classification choice
    """
    
    list_of_column_choices = ['car_m_t', 'car_r_t', 'from_id', 'pt_m_tt', 'pt_r_tt', 'walk_t']
    
    random_column_choice = random.choice(list_of_column_choices)
    
    #add classification to gpd
    classifier = classification(y=acc[random_column_choice], k=7)
    classifier = classification.make(k=7)
    
    # add a new 'temp' column to store classified
    gdf['temp'] = gdf[[random_column_choice]].apply(classifier) 
    
    #basic plot config
    fig, (ax1,ax2, ax3) = plt.subplots(nrows=3, ncols=1, figsize=(15, 11))
    #ax1, ax2 = axes

    #1st plot with random classsification
    ax2 = gdf.plot(ax=ax2, column='temp', linewidth=0, legend=True)
    
    #set title
    printTitle = ''
    if random_column_choice == 'car_m_t':
        printTitle = 'Travel time in minutes from origin to destination \nby private car in midday traffic'
    if random_column_choice == 'car_r_t':
        printTitle = 'Travel time in minutes from origin to destination \nby private car in rush hour traffic'
    if random_column_choice == 'pt_m_tt':
        printTitle = 'Travel time in minutes from origin to destination \nby public transportation in midday traffic'
    if random_column_choice == 'pt_r_tt':
        printTitle = 'Travel time in minutes from origin to destination \nby public transportation in rush hour traffic'
    if random_column_choice == 'walk_t':
        printTitle = 'Travel time in minutes from origin to destination \nby walking'
    ax2.set_title(printTitle, pad=20)

    #2nd plot with original data
    ax1 = gdf.plot(ax=ax1, column="pt_r_tt", linewidth=0, legend=True)
    ax1.set_title('Original')

    # Plot histogram for public transport rush hour travel time
    ax3 = gdf[random_column_choice].plot.hist(ax=ax3, bins=50)
    
    classifier2 = classification(y=acc[random_column_choice])

    # Add vertical lines for class breaks
    for value in classifier2.bins:
        plt.axvline(value, color='k', linestyle='dashed', linewidth=1)
        
    fig.tight_layout()
classify_and_plot(acc, classification = mapclassify.EqualInterval)<jupyter_output><empty_output><jupyter_text>Creating a custom classifier**Multicriteria data classification**Let's create a function where we classify the geometries into two classes based on a given `threshold` -parameter. If the area of a polygon is lower than the threshold value (average size of the lake), the output column will get a value 0, if it is larger, it will get a value 1. This kind of classification is often called a [binary classification](https://en.wikipedia.org/wiki/Binary_classification).First we need to create a function for our classification task. This function takes a single row of the GeoDataFrame as input, plus few other parameters that we can use.It also possible to do classifiers with multiple criteria easily in Pandas/Geopandas by extending the example that we started earlier. Now we will modify our binaryClassifier function a bit so that it classifies the data based on two columns.- Let's call it `custom_classifier` that does the binary classification based on two treshold values:<jupyter_code>def custom_classifier(row, src_col1, src_col2, threshold1, threshold2, output_col):
   # 1. If the value in src_col1 is LOWER than the threshold1 value
   # 2. AND the value in src_col2 is HIGHER than the threshold2 value, give value 1, otherwise give 0
   if row[src_col1] < threshold1 and row[src_col2] > threshold2:
       # Update the output column with value 0
       row[output_col] = 1
   # If area of input geometry is higher than the threshold value update with value 1
   else:
       row[output_col] = 0

   # Return the updated row
   return row<jupyter_output><empty_output><jupyter_text>Now we have defined the function, and we can start using it.- Let's do our classification based on two criteria and find out grid cells where the **travel time is lower or equal to 20 minutes** but they are further away **than 4 km (4000 meters) from city center**.- Let's create an empty column for our classification results called `"suitable_area"`.<jupyter_code># Create column for the classification results
acc["suitable_area"] = None

# Use the function
acc = acc.apply(custom_classifier, src_col1='pt_r_tt', 
                src_col2='walk_d', threshold1=20, threshold2=4000, 
                output_col="suitable_area", axis=1)

# See the first rows
acc.head(2)<jupyter_output><empty_output><jupyter_text>Okey we have new values in `suitable_area` -column.- How many Polygons are suitable for us? Let's find out by using a Pandas function called `value_counts()` that return the count of different values in our column.<jupyter_code># Get value counts
acc['suitable_area'].value_counts()<jupyter_output><empty_output><jupyter_text>Okay, so there seems to be nine suitable locations for us where we can try to find an appartment to buy.- Let's see where they are located:<jupyter_code># Plot
acc.plot(column="suitable_area", linewidth=0);

# Use tight layour
plt.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08602846054333764, 'fraction_numerical': 0.02483829236739974, 'mean_word_length': 4.771183277342292, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 84, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '63185', 'license': 'MIT', 'n_tokens_mistral': 5015, 'n_tokens_neox': 4603, 'n_words': 1667, 'path': 'source/notebooks/L4/reclassify.ipynb', 'repo_name': 'AdrianKriger/Automating-GIS-Processess', 'size': '15460'}