<jupyter_start><jupyter_text>Spatial join[Spatial join](http://wiki.gis.com/wiki/index.php/Spatial_Join) isyet another classic GIS problem. Getting attributes from one layer andtransferring them into another layer based on their spatial relationshipis something you most likely need to do on a regular basis.In the previous section we learned how to perform **a Point in Polygon query**.We could now apply those techniques and create our own function to perform **a spatial join** between two layers based on theirspatial relationship. We could, for example, join the attributes of a polygon layer into a point layer where each point would get theattributes of a polygon that ``contains`` the point.Luckily, [spatial join is already implemented in Geopandas](http://geopandas.org/mergingdata.htmlspatial-joins), thus we do not need to create it ourselves. There are three possible types ofjoin that can be applied in spatial join that are determined with ``op`` -parameter in the ``gpd.sjoin()`` -function:-  ``"intersects"``-  ``"within"``-  ``"contains"``Sounds familiar? Yep, all of those spatial relationships were discussedin the [Point in Polygon lesson](point-in-polygon.ipynb), thus you should know how they work.Let's perform a spatial join between these two layers:- **Addresses:** the address-point Shapefile that we created trough geocoding - **Population grid:** a Polygon layer that is a 250m x 250m grid showing the amount of people living in the Helsinki Region.    - The population grid a dataset is produced by the **Helsinki Region EnvironmentalServices Authority (HSY)** (see [this page](https://www.hsy.fi/fi/asiantuntijalle/avoindata/Sivut/AvoinData.aspx?dataID=7) to access data from different years).    - For this lesson we will use the population grid for year 2017, which can be dowloaded as a shapefile [from this link](https://www.hsy.fi/sites/AvoinData/AvoinData/SYT/Tietoyhteistyoyksikko/Shape%20(Esri)/V%C3%A4est%C3%B6tietoruudukko/Vaestotietoruudukko_2017_SHP.zip) in the  [Helsinki Region Infroshare(HRI) open data portal](https://hri.fi/en_gb/)  Download and clean the data**Execute the following steps in a terminal window**- Navigate to the data folder```    $ cd data```- Download the population grid using wget:```    $ wget "https://www.hsy.fi/sites/AvoinData/AvoinData/SYT/Tietoyhteistyoyksikko/Shape%20(Esri)/V%C3%A4est%C3%B6tietoruudukko/Vaestotietoruudukko_2017_SHP.zip"```-  Unzip the file in Terminal into a folder called Pop17 (using -d flag)```    $ unzip Vaestotietoruudukko_2017_SHP.zip -d Pop17```You should now have a folder ``/data/Pop17`` containing the population grid shapefile.-  Let's read the data into memory and see what we have.<jupyter_code>import geopandas as gpd

# Filepath
fp = "data/Pop17/Vaestoruudukko_2017.shp"

# Read the data
pop = gpd.read_file(fp)

# See the first rows
pop.head()<jupyter_output><empty_output><jupyter_text>Okey so we have multiple columns in the dataset but the most importantone here is the column ``ASUKKAITA`` (*"population" in Finnish*) thattells the amount of inhabitants living under that polygon.-  Let's change the name of that columns into ``pop17`` so that it is   more intuitive. Changing column names is easy in Pandas / Geopandas   using a function called ``rename()`` where we pass a dictionary to a   parameter ``columns={'oldname': 'newname'}``.<jupyter_code># Change the name of a column
pop = pop.rename(columns={'ASUKKAITA': 'pop17'})

# See the column names and confirm that we now have a column called 'pop17'
pop.columns<jupyter_output><empty_output><jupyter_text>-  Let's also get rid of all unnecessary columns by selecting only   columns that we need i.e. ``pop17`` and ``geometry``<jupyter_code># Columns that will be sected
selected_cols = ['pop17', 'geometry']

# Select those columns
pop = pop[selected_cols]

# Let's see the last 2 rows
pop.head()<jupyter_output><empty_output><jupyter_text>Now we have cleaned the data and have only those columns that we needfor our analysis.  Join the layersNow we are ready to perform the spatial join between the two layers thatwe have. The aim here is to get information about **how many people livein a polygon that contains an individual address-point** . Thus, we wantto join attributes from the population layer we just modified into theaddresses point layer ``addresses.shp`` that we created trough gecoding in the previous section.-  Read the addresses layer into memory<jupyter_code># Addresses filpath
addr_fp = r"data/addresses.shp"

# Read data
addresses = gpd.read_file(addr_fp)

# Check the head of the file
addresses.head()<jupyter_output><empty_output><jupyter_text>In order to do a spatial join, the layers need to be in the same projection<jupyter_code># Are the layers in the same projection?
addresses.crs == pop.crs<jupyter_output><empty_output><jupyter_text>Let's re-project addresses to the projection of the population layer:<jupyter_code>addresses = addresses.to_crs(pop.crs)<jupyter_output><empty_output><jupyter_text>-  Let's make sure that the coordinate reference system of the layersare identical<jupyter_code># Check the crs of address points
print(addresses.crs)

# Check the crs of population layer
print(pop.crs)

# Do they match now?
addresses.crs == pop.crs<jupyter_output>{'proj': 'tmerc', 'lat_0': 0, 'lon_0': 25, 'k': 1, 'x_0': 25500000, 'y_0': 0, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}
{'proj': 'tmerc', 'lat_0': 0, 'lon_0': 25, 'k': 1, 'x_0': 25500000, 'y_0': 0, 'ellps': 'GRS80', 'units': 'm', 'no_defs': True}<jupyter_text>Now they should be identical. Thus, we can be sure that when doing spatialqueries between layers the locations match and we get the right resultse.g. from the spatial join that we are conducting here.-  Let's now join the attributes from ``pop`` GeoDataFrame into   ``addresses`` GeoDataFrame by using ``gpd.sjoin()`` -function<jupyter_code># Make a spatial join
join = gpd.sjoin(addresses, pop, how="inner", op="within")

# Let's check the result
join.head()<jupyter_output><empty_output><jupyter_text>Awesome! Now we have performed a successful spatial join where we gottwo new columns into our ``join`` GeoDataFrame, i.e. ``index_right``that tells the index of the matching polygon in the population grid and``pop17`` which is the population in the cell where the address-point islocated.-  Let's save this layer into a new Shapefile<jupyter_code># Output path
outfp = r"data/addresses_pop17_epsg3979.shp"

# Save to disk
join.to_file(outfp)<jupyter_output><empty_output><jupyter_text>Do the results make sense? Let's evaluate this a bit by plotting thepoints where color intensity indicates the population numbers.-  Plot the points and use the ``pop17`` column to indicate the color.   ``cmap`` -parameter tells to use a sequential colormap for the   values, ``markersize`` adjusts the size of a point, ``scheme`` parameter can be used to adjust the classification method based on [pysal](http://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html), and ``legend`` tells that we want to have a legend.<jupyter_code>import matplotlib.pyplot as plt

# Plot the points with population info
join.plot(column='pop17', cmap="Reds", markersize=7, scheme='quantiles', legend=True);

# Add title
plt.title("Amount of inhabitants living close the the point");

# Remove white space around the figure
plt.tight_layout()<jupyter_output>/opt/conda/lib/python3.6/site-packages/scipy/stats/stats.py:1713: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
  return np.add.reduce(sorted[indexer] * weights, axis=axis) / sumval
{'fraction_non_alphanumeric': 0.09818789358694255, 'fraction_numerical': 0.014779591312170673, 'mean_word_length': 5.447390223695112, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 39, 'https://': 4, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '419178', 'license': 'MIT', 'n_tokens_mistral': 2495, 'n_tokens_neox': 2353, 'n_words': 929, 'path': 'source/notebooks/L3/spatial-join.ipynb', 'repo_name': 'ropitz/spatial_teaching_fundamentals', 'size': '7781'}