-
Notifications
You must be signed in to change notification settings - Fork 28
/
Copy pathLesson3-spatial-join.html
573 lines (445 loc) · 32.3 KB
/
Lesson3-spatial-join.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
<meta charset="utf-8">
<script type="text/javascript">
var _gaq = _gaq || [];
_gaq.push(['_setAccount', 'UA-88382509-1']);
_gaq.push(['_trackPageview']);
(function() {
var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true;
ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js';
var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s);
})();
</script>
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Spatial join — GeoPython - AutoGIS 1 documentation</title>
<link rel="stylesheet" href="_static/css/theme.css" type="text/css" />
<link rel="index" title="Index"
href="genindex.html"/>
<link rel="search" title="Search" href="search.html"/>
<link rel="top" title="GeoPython - AutoGIS 1 documentation" href="index.html"/>
<link rel="next" title="Exercise 3" href="Exercise-3.html"/>
<link rel="prev" title="Point in Polygon & Intersect" href="Lesson3-point-in-polygon.html"/>
<script src="_static/js/modernizr.min.js"></script>
</head>
<body class="wy-body-for-nav" role="document">
<div class="wy-grid-for-nav">
<nav data-toggle="wy-nav-shift" class="wy-nav-side">
<div class="wy-side-scroll">
<div class="wy-side-nav-search">
<a href="index.html" class="icon icon-home"> GeoPython - AutoGIS
<img src="_static/logo_hy_geo_135.png" class="logo" />
</a>
<div class="version">
2016 Autumn
</div>
<div role="search">
<form id="rtd-search-form" class="wy-form" action="search.html" method="get">
<input type="text" name="q" placeholder="Search docs" />
<input type="hidden" name="check_keywords" value="yes" />
<input type="hidden" name="area" value="default" />
</form>
</div>
</div>
<div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
<p class="caption"><span class="caption-text">Course information</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="course-info.html">General info</a></li>
<li class="toctree-l1"><a class="reference internal" href="Installing_Anacondas_GIS.html">Installing Python + GIS</a></li>
<li class="toctree-l1"><a class="reference internal" href="License-terms.html">License and terms of usage</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 1</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson1-Intro-Python-GIS.html">Introduction to Python GIS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson1-Geometric-Objects.html">Geometric Objects - Spatial Data Model</a></li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-1.html">Exercise 1</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 2</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson2-overview.html">Lesson 2 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson2-download-data.html">Download datasets</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson2-overview-pandas-geopandas.html">Pandas and Geopandas -modules</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson2-pandas-intro.html">Introduction to Pandas</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson2-geopandas-basics.html">Introduction to Geopandas</a></li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-2.html">Exercise 2</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 3</span></p>
<ul class="current">
<li class="toctree-l1"><a class="reference internal" href="Lesson3-overview.html">Lesson 3 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson3-geocoding.html">Geocoding</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson3-table-join.html">Table join</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson3-projections.html">Re-projecting data</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson3-point-in-polygon.html">Point in Polygon & Intersect</a></li>
<li class="toctree-l1 current"><a class="current reference internal" href="#">Spatial join</a><ul>
<li class="toctree-l2"><a class="reference internal" href="#download-and-clean-the-data">Download and clean the data</a></li>
<li class="toctree-l2"><a class="reference internal" href="#join-the-layers">Join the layers</a></li>
</ul>
</li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-3.html">Exercise 3</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 4</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson4-overview.html">Lesson 4 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson4-download-data.html">Download datasets</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson4-geometric-operations.html">Geometric operations</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson4-reclassify.html">Data reclassification</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson4-nearest-neighbour.html">Nearest Neighbour Analysis</a></li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-4.html">Exercise 4</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 5</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-overview.html">Lesson 5 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-download-data.html">Download datasets</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-static-maps.html">Static maps</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-interactive-map-bokeh.html">Interactive maps with Bokeh</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-share-on-github.html">Sharing interactive plots on GitHub</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-interactive-map-folium.html">Interactive maps on Leaflet</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson5-World-3D.html">Inspiration: World 3D</a></li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-5.html">Exercise 5</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 6</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson6-overview.html">Lesson 6 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson6-arcpy.html">Python in ArcGIS</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson6-toolbox.html">ArcGIS Toolbox</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson6-arcpy-script.html">Writing arcpy scripts</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson6-run-the-tool.html">Running the Python script from ArcGIS</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 7</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson7-overview.html">Lesson 7 Overview</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson7-download.html">Download data</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson7-read-raster.html">Reading raster files with GDAL</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson7-read-raster-array.html">Reading raster as a numerical array</a></li>
<li class="toctree-l1"><a class="reference internal" href="Lesson7-gdal-utilities.html">GDAL command line tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="Exercise-7.html">Exercise 7</a></li>
</ul>
<p class="caption"><span class="caption-text">Lesson 8</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Lesson8-network-analysis.html">Network analysis in Python</a></li>
</ul>
<p class="caption"><span class="caption-text">Final Assignment</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="Final-assignment.html">Final assignment</a></li>
</ul>
<p class="caption"><span class="caption-text">Map Challenge 2016</span></p>
<ul>
<li class="toctree-l1"><a class="reference internal" href="map-challenge.html">Map Challenge 2016</a></li>
</ul>
</div>
</div>
</nav>
<section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">
<nav class="wy-nav-top" role="navigation" aria-label="top navigation">
<i data-toggle="wy-nav-top" class="fa fa-bars"></i>
<a href="index.html">GeoPython - AutoGIS</a>
</nav>
<div class="wy-nav-content">
<div class="rst-content">
<div role="navigation" aria-label="breadcrumbs navigation">
<ul class="wy-breadcrumbs">
<li><a href="index.html">Docs</a> »</li>
<li>Spatial join</li>
<li class="wy-breadcrumbs-aside">
<a href="https://github.com/Automating-GIS-processes/2016/blob/master/source/Lesson3-spatial-join.rst" class="fa fa-github"> Edit on GitHub</a>
</li>
</ul>
<hr/>
</div>
<div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
<div itemprop="articleBody">
<div class="section" id="spatial-join">
<h1>Spatial join<a class="headerlink" href="#spatial-join" title="Permalink to this headline">¶</a></h1>
<p><strong>Sources</strong></p>
<p><em>Following materials are partly based on documentation of</em> <a class="reference external" href="http://geopandas.org">Geopandas</a>.</p>
<p><a class="reference external" href="http://wiki.gis.com/wiki/index.php/Spatial_Join">Spatial join</a> is
yet another classic GIS problem. Getting attributes from one layer and
transferring them into another layer based on their spatial relationship
is something you most likely need to do on a regular basis.</p>
<p>The previous materials focused on learning how to perform a <a class="reference external" href="Lesson3-point-in-polygon.html#how-to-check-if-point-is-inside-a-polygon">Point in Polygon query</a>.
We could now apply those techniques and create our
own function to perform a spatial join between two layers based on their
spatial relationship. We could for example join the attributes of a
polygon layer into a point layer where each point would get the
attributes of a polygon that <code class="docutils literal"><span class="pre">contains</span></code> the point.</p>
<p>Luckily, <a class="reference external" href="http://geopandas.org/mergingdata.html#spatial-joins">spatial join</a>
(<code class="docutils literal"><span class="pre">gpd.sjoin()</span></code> -function) is already implemented in Geopandas, thus we
do not need to create it ourselves. There are three possible types of
join that can be applied in spatial join that are determined with <code class="docutils literal"><span class="pre">op</span></code>
-parameter:</p>
<ul class="simple">
<li><code class="docutils literal"><span class="pre">"intersects"</span></code></li>
<li><code class="docutils literal"><span class="pre">"within"</span></code></li>
<li><code class="docutils literal"><span class="pre">"contains"</span></code></li>
</ul>
<p>Sounds familiar? Yep, all of those spatial relationships were discussed
in the <a class="reference external" href="Lesson3-point-in-polygon.html">previous materials</a>, thus you should know how they work.</p>
<p>Let’s perform a spatial join between the address-point Shapefile that we
<a class="reference external" href="Lesson3-table-join.html">created</a> and then <a class="reference external" href="Lesson3-projections.html">reprojected</a>
and a Polygon layer that is a
250m x 250m grid showing the amount of people living in Helsinki Region.</p>
<div class="section" id="download-and-clean-the-data">
<h2>Download and clean the data<a class="headerlink" href="#download-and-clean-the-data" title="Permalink to this headline">¶</a></h2>
<p>For this lesson we will be using publicly available population data from
Helsinki that can be downloaded from <a class="reference external" href="http://www.hri.fi/en/dataset/vaestotietoruudukko">Helsinki Region Infroshare
(HRI)</a> which is an
excellent source that provides all sorts of open data from Helsinki,
Finland.</p>
<p>From HRI <strong>download a</strong> <a class="reference external" href="https://www.hsy.fi/sites/AvoinData/AvoinData/SYT/Tietoyhteistyoyksikko/Shape%20(Esri)/V%C3%A4est%C3%B6tietoruudukko/Vaestotietoruudukko_2015.zip">Population grid for year
2015</a>
that is a dataset (.shp) produced by Helsinki Region Environmental
Services Authority (HSY) (see <a class="reference external" href="https://www.hsy.fi/fi/asiantuntijalle/avoindata/Sivut/AvoinData.aspx?dataID=7">this
page</a>
to access data from different years).</p>
<ul class="simple">
<li>Unzip the file in Terminal into a folder called Pop15 (using -d flag)</li>
</ul>
<div class="code bash highlight-default"><div class="highlight"><pre><span></span>$ cd
$ unzip Vaestotietoruudukko_2015.zip -d Pop15
$ ls Pop15
Vaestotietoruudukko_2015.dbf Vaestotietoruudukko_2015.shp
Vaestotietoruudukko_2015.prj Vaestotietoruudukko_2015.shx
</pre></div>
</div>
<p>You should now have a folder <code class="docutils literal"><span class="pre">/home/geo/Pop15</span></code> with files listed
above.</p>
<ul class="simple">
<li>Let’s read the data into memory and see what we have.</li>
</ul>
<div class="code python highlight-default"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">geopandas</span> <span class="k">as</span> <span class="nn">gpd</span>
<span class="c1"># Filepath</span>
<span class="n">fp</span> <span class="o">=</span> <span class="s2">"/home/geo/Pop15/Vaestotietoruudukko_2015.shp"</span>
<span class="c1"># Read the data</span>
<span class="n">pop</span> <span class="o">=</span> <span class="n">gpd</span><span class="o">.</span><span class="n">read_file</span><span class="p">(</span><span class="n">fp</span><span class="p">)</span>
</pre></div>
</div>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># See the first rows</span>
<span class="gp">In [1]: </span><span class="n">pop</span><span class="o">.</span><span class="n">head</span><span class="p">()</span>
<span class="gh">Out[1]: </span><span class="go"></span>
<span class="go"> INDEX ASUKKAITA ASVALJYYS IKA0_9 IKA10_19 IKA20_29 IKA30_39 \</span>
<span class="go">0 688 8 31.0 99 99 99 99 </span>
<span class="go">1 703 6 42.0 99 99 99 99 </span>
<span class="go">2 710 8 44.0 99 99 99 99 </span>
<span class="go">3 711 7 64.0 99 99 99 99 </span>
<span class="go">4 715 19 23.0 99 99 99 99 </span>
<span class="go"> IKA40_49 IKA50_59 IKA60_69 IKA70_79 IKA_YLI80 \</span>
<span class="go">0 99 99 99 99 99 </span>
<span class="go">1 99 99 99 99 99 </span>
<span class="go">2 99 99 99 99 99 </span>
<span class="go">3 99 99 99 99 99 </span>
<span class="go">4 99 99 99 99 99 </span>
<span class="go"> geometry </span>
<span class="go">0 POLYGON ((25472499.99532626 6689749.005069185,... </span>
<span class="go">1 POLYGON ((25472499.99532626 6685998.998064222,... </span>
<span class="go">2 POLYGON ((25472499.99532626 6684249.004130407,... </span>
<span class="go">3 POLYGON ((25472499.99532626 6683999.004997005,... </span>
<span class="go">4 POLYGON ((25472499.99532626 6682998.998461431,... </span>
</pre></div>
</div>
<p>Okey so we have multiple columns in the dataset but the most important
one here is the column <code class="docutils literal"><span class="pre">ASUKKAITA</span></code> (<em>population in Finnish</em>) that
tells the amount of inhabitants living under that polygon.</p>
<ul class="simple">
<li>Let’s change the name of that columns into <code class="docutils literal"><span class="pre">pop15</span></code> so that it is
more intuitive. Changing column names is easy in Pandas / Geopandas
using a function called <code class="docutils literal"><span class="pre">rename()</span></code> where we pass a dictionary to a
parameter <code class="docutils literal"><span class="pre">columns={'oldname':</span> <span class="pre">'newname'}</span></code>.</li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># Change the name of a column</span>
<span class="gp">In [2]: </span><span class="n">pop</span> <span class="o">=</span> <span class="n">pop</span><span class="o">.</span><span class="n">rename</span><span class="p">(</span><span class="n">columns</span><span class="o">=</span><span class="p">{</span><span class="s1">'ASUKKAITA'</span><span class="p">:</span> <span class="s1">'pop15'</span><span class="p">})</span>
<span class="go"># See the column names and confirm that we now have a column called 'pop15'</span>
<span class="gp">In [3]: </span><span class="n">pop</span><span class="o">.</span><span class="n">columns</span>
<span class="gh">Out[3]: </span><span class="go"></span>
<span class="go">Index(['INDEX', 'pop15', 'ASVALJYYS', 'IKA0_9', 'IKA10_19', 'IKA20_29',</span>
<span class="go"> 'IKA30_39', 'IKA40_49', 'IKA50_59', 'IKA60_69', 'IKA70_79', 'IKA_YLI80',</span>
<span class="go"> 'geometry'],</span>
<span class="go"> dtype='object')</span>
</pre></div>
</div>
<ul class="simple">
<li>Let’s also get rid of all unnecessary columns by selecting only
columns that we need i.e. <code class="docutils literal"><span class="pre">pop15</span></code> and <code class="docutils literal"><span class="pre">geometry</span></code></li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># Columns that will be sected</span>
<span class="gp">In [4]: </span><span class="n">selected_cols</span> <span class="o">=</span> <span class="p">[</span><span class="s1">'pop15'</span><span class="p">,</span> <span class="s1">'geometry'</span><span class="p">]</span>
<span class="go"># Select those columns</span>
<span class="gp">In [5]: </span><span class="n">pop</span> <span class="o">=</span> <span class="n">pop</span><span class="p">[</span><span class="n">selected_cols</span><span class="p">]</span>
<span class="go"># Let's see the last 2 rows</span>
<span class="gp">In [6]: </span><span class="n">pop</span><span class="o">.</span><span class="n">tail</span><span class="p">(</span><span class="mi">2</span><span class="p">)</span>
<span class="gh">Out[6]: </span><span class="go"></span>
<span class="go"> pop15 geometry</span>
<span class="go">5782 9 POLYGON ((25513499.99632164 6685498.999797418,...</span>
<span class="go">5783 30244 POLYGON ((25513999.999929 6659998.998172711, 2...</span>
</pre></div>
</div>
<p>Now we have cleaned the data and have only those columns that we need
for our analysis.</p>
</div>
<div class="section" id="join-the-layers">
<h2>Join the layers<a class="headerlink" href="#join-the-layers" title="Permalink to this headline">¶</a></h2>
<p>Now we are ready to perform the spatial join between the two layers that
we have. The aim here is to get information about <strong>how many people live
in a polygon that contains an individual address-point</strong> . Thus, we want
to join attributes from the population layer we just modified into the
addresses point layer <code class="docutils literal"><span class="pre">addresses_epsg3879.shp</span></code>.</p>
<ul class="simple">
<li>Read the addresses layer into memory</li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># Addresses filpath</span>
<span class="gp">In [7]: </span><span class="n">addr_fp</span> <span class="o">=</span> <span class="s2">r"/home/geo/addresses_epsg3879.shp"</span>
<span class="go"># Read data</span>
<span class="gp">In [8]: </span><span class="n">addresses</span> <span class="o">=</span> <span class="n">gpd</span><span class="o">.</span><span class="n">read_file</span><span class="p">(</span><span class="n">addr_fp</span><span class="p">)</span>
<span class="go"># Check the head of the file</span>
<span class="gp">In [9]: </span><span class="n">addresses</span><span class="o">.</span><span class="n">head</span><span class="p">(</span><span class="mi">2</span><span class="p">)</span>
<span class="gh">Out[9]: </span><span class="go"></span>
<span class="go"> address id \</span>
<span class="go">0 Kampinkuja 1, 00100 Helsinki, Finland 1001 </span>
<span class="go">1 Kaivokatu 8, 00101 Helsinki, Finland 1002 </span>
<span class="go"> geometry </span>
<span class="go">0 POINT (25496123.30852197 6672833.941567578) </span>
<span class="go">1 POINT (25496774.28242895 6672999.698581985) </span>
</pre></div>
</div>
<ul class="simple">
<li>Let’s make sure that the coordinate reference system of the layers
are identical</li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># Check the crs of address points</span>
<span class="gp">In [10]: </span><span class="n">addresses</span><span class="o">.</span><span class="n">crs</span>
<span class="gh">Out[10]: </span><span class="go"></span>
<span class="go">{'ellps': 'GRS80',</span>
<span class="go"> 'k': 1,</span>
<span class="go"> 'lat_0': 0,</span>
<span class="go"> 'lon_0': 25,</span>
<span class="go"> 'no_defs': True,</span>
<span class="go"> 'proj': 'tmerc',</span>
<span class="go"> 'units': 'm',</span>
<span class="go"> 'x_0': 25500000,</span>
<span class="go"> 'y_0': 0}</span>
<span class="go"># Check the crs of population layer</span>
<span class="gp">In [11]: </span><span class="n">pop</span><span class="o">.</span><span class="n">crs</span>
<span class="go"> Out[11]: </span>
<span class="go">{'ellps': 'GRS80',</span>
<span class="go"> 'k': 1,</span>
<span class="go"> 'lat_0': 0,</span>
<span class="go"> 'lon_0': 25,</span>
<span class="go"> 'no_defs': True,</span>
<span class="go"> 'proj': 'tmerc',</span>
<span class="go"> 'units': 'm',</span>
<span class="go"> 'x_0': 25500000,</span>
<span class="go"> 'y_0': 0}</span>
<span class="go"># Do they match? - We can test that</span>
<span class="gp">In [12]: </span><span class="n">addresses</span><span class="o">.</span><span class="n">crs</span> <span class="o">==</span> <span class="n">pop</span><span class="o">.</span><span class="n">crs</span>
<span class="go"> Out[12]: True</span>
</pre></div>
</div>
<p>Indeed they are identical. Thus, we can be sure that when doing spatial
queries between layers the locations match and we get the right results
e.g. from the spatial join that we are conducting here.</p>
<ul class="simple">
<li>Let’s now join the attributes from <code class="docutils literal"><span class="pre">pop</span></code> GeoDataFrame into
<code class="docutils literal"><span class="pre">addresses</span></code> GeoDataFrame by using <code class="docutils literal"><span class="pre">gpd.sjoin()</span></code> -function</li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="go"># Make a spatial join</span>
<span class="gp">In [13]: </span><span class="n">join</span> <span class="o">=</span> <span class="n">gpd</span><span class="o">.</span><span class="n">sjoin</span><span class="p">(</span><span class="n">addresses</span><span class="p">,</span> <span class="n">pop</span><span class="p">,</span> <span class="n">how</span><span class="o">=</span><span class="s2">"inner"</span><span class="p">,</span> <span class="n">op</span><span class="o">=</span><span class="s2">"within"</span><span class="p">)</span>
<span class="go"># Let's check the result</span>
<span class="gp">In [14]: </span><span class="n">join</span><span class="o">.</span><span class="n">head</span><span class="p">()</span>
<span class="gh">Out[14]: </span><span class="go"></span>
<span class="go"> address id \</span>
<span class="go">16 Malminkartanontie 17, 00410 Helsinki, Finland 1017 </span>
<span class="go">19 Pitäjänmäentie 15, 00370 Helsinki, Finland 1020 </span>
<span class="go">12 Trumstigen 8, 00420 Helsingfors, Finland 1013 </span>
<span class="go">11 Kuparitie 8, 00440 Helsinki, Finland 1012 </span>
<span class="go">15 Kylätie 23, 00320 Helsinki, Finland 1016 </span>
<span class="go"> geometry index_right pop15 </span>
<span class="go">16 POINT (25492349.68368251 6681772.551210108) 2684 74 </span>
<span class="go">19 POINT (25492292.61413005 6679039.264838208) 2691 241 </span>
<span class="go">12 POINT (25493207.62503373 6680836.727432437) 2852 577 </span>
<span class="go">11 POINT (25493575.10327127 6679775.868274149) 2949 562 </span>
<span class="go">15 POINT (25494077.1680778 6678341.639159317) 3036 414 </span>
</pre></div>
</div>
<p>Awesome! Now we have performed a successful spatial join where we got
two new columns into our <code class="docutils literal"><span class="pre">join</span></code> GeoDataFrame, i.e. <code class="docutils literal"><span class="pre">index_right</span></code>
that tells the index of the matching polygon in the <code class="docutils literal"><span class="pre">pop</span></code> layer and
<code class="docutils literal"><span class="pre">pop15</span></code> which is the population in the cell where the address-point is
located.</p>
<ul class="simple">
<li>Let’s save this layer into a new Shapefile</li>
</ul>
<div class="code python highlight-default"><div class="highlight"><pre><span></span><span class="c1"># Output path</span>
<span class="n">outfp</span> <span class="o">=</span> <span class="s2">r"/home/geo/addresses_pop15_epsg3979.shp"</span>
<span class="c1"># Save to disk</span>
<span class="n">join</span><span class="o">.</span><span class="n">to_file</span><span class="p">(</span><span class="n">outfp</span><span class="p">)</span>
</pre></div>
</div>
<p>Do the results make sense? Let’s evaluate this a bit by plotting the
points where color intensity indicates the population numbers.</p>
<ul class="simple">
<li>Plot the points and use the <code class="docutils literal"><span class="pre">pop15</span></code> column to indicate the color.
<code class="docutils literal"><span class="pre">cmap</span></code> -parameter tells to use a sequential colormap for the
values, <code class="docutils literal"><span class="pre">markersize</span></code> adjusts the size of a point, <code class="docutils literal"><span class="pre">scheme</span></code> parameter can be used to adjust the classification method based on <a class="reference external" href="http://pysal.readthedocs.io/en/latest/library/esda/mapclassify.html">pysal</a>, and <code class="docutils literal"><span class="pre">legend</span></code> tells that we want to have a legend.</li>
</ul>
<div class="highlight-ipython"><div class="highlight"><pre><span></span><span class="gp">In [15]: </span><span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="kn">as</span> <span class="nn">plt</span>
<span class="go"># Plot the points with population info</span>
<span class="gp">In [16]: </span><span class="n">join</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">column</span><span class="o">=</span><span class="s1">'pop15'</span><span class="p">,</span> <span class="n">cmap</span><span class="o">=</span><span class="s2">"Reds"</span><span class="p">,</span> <span class="n">markersize</span><span class="o">=</span><span class="mi">7</span><span class="p">,</span> <span class="n">scheme</span><span class="o">=</span><span class="s1">'natural_breaks'</span><span class="p">,</span> <span class="n">legend</span><span class="o">=</span><span class="bp">True</span><span class="p">);</span>
<span class="go"># Add title</span>
<span class="gp">In [17]: </span><span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="s2">"Amount of inhabitants living close the the point"</span><span class="p">);</span>
<span class="go"># Remove white space around the figure</span>
<span class="gp">In [18]: </span><span class="n">plt</span><span class="o">.</span><span class="n">tight_layout</span><span class="p">()</span>
</pre></div>
</div>
<a class="reference internal image-reference" href="_images/population_points.png"><img alt="_images/population_points.png" src="_images/population_points.png" style="width: 7in;" /></a>
<p>By knowing approximately how population is distributed in Helsinki, it
seems that the results do make sense as the points with highest
population are located in the south where the city center of Helsinki
is.</p>
</div>
</div>
</div>
</div>
<footer>
<div class="rst-footer-buttons" role="navigation" aria-label="footer navigation">
<a href="Exercise-3.html" class="btn btn-neutral float-right" title="Exercise 3" accesskey="n">Next <span class="fa fa-arrow-circle-right"></span></a>
<a href="Lesson3-point-in-polygon.html" class="btn btn-neutral" title="Point in Polygon & Intersect" accesskey="p"><span class="fa fa-arrow-circle-left"></span> Previous</a>
</div>
<hr/>
<div role="contentinfo">
<p>
© Copyright 2016, Henrikki Tenkanen.
Last updated on Feb 20, 2017.
</p>
</div>
Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/snide/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>.
</footer>
</div>
</div>
</section>
</div>
<script type="text/javascript">
var DOCUMENTATION_OPTIONS = {
URL_ROOT:'./',
VERSION:'1',
COLLAPSE_INDEX:false,
FILE_SUFFIX:'.html',
HAS_SOURCE: true
};
</script>
<script type="text/javascript" src="_static/jquery.js"></script>
<script type="text/javascript" src="_static/underscore.js"></script>
<script type="text/javascript" src="_static/doctools.js"></script>
<script type="text/javascript" src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
<script type="text/javascript" src="_static/js/theme.js"></script>
<script type="text/javascript">
jQuery(function () {
SphinxRtdTheme.StickyNav.enable();
});
</script>
<li><a href="http://project.invalid/">Project Homepage</a> »</li>
<div class="footer">
<img src="../img/GPLv3_Logo.svg">
</div>
</body>
</html>