-
-
Notifications
You must be signed in to change notification settings - Fork 69
/
Copy path06-reproj.html
1811 lines (1766 loc) · 150 KB
/
06-reproj.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
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" lang="en" xml:lang="en"><head>
<meta charset="utf-8">
<meta name="generator" content="quarto-1.5.54">
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
<meta name="description" content="An introductory resource for working with geographic data in Python">
<title>6 Reprojecting geographic data – Geocomputation with Python</title>
<style>
code{white-space: pre-wrap;}
span.smallcaps{font-variant: small-caps;}
div.columns{display: flex; gap: min(4vw, 1.5em);}
div.column{flex: auto; overflow-x: auto;}
div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
ul.task-list{list-style: none;}
ul.task-list li input[type="checkbox"] {
width: 0.8em;
margin: 0 0.8em 0.2em -1em; /* quarto-specific, see https://github.com/quarto-dev/quarto-cli/issues/4556 */
vertical-align: middle;
}
/* CSS for syntax highlighting */
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
div.sourceCode { margin: 1em 0; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { display: inline-block; text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
{ counter-reset: source-line 0; }
pre.numberSource code > span
{ position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
{ content: counter(source-line);
position: relative; left: -1em; text-align: right; vertical-align: baseline;
border: none; display: inline-block;
-webkit-touch-callout: none; -webkit-user-select: none;
-khtml-user-select: none; -moz-user-select: none;
-ms-user-select: none; user-select: none;
padding: 0 4px; width: 4em;
}
pre.numberSource { margin-left: 3em; padding-left: 4px; }
div.sourceCode
{ }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
/* CSS for citations */
div.csl-bib-body { }
div.csl-entry {
clear: both;
margin-bottom: 0em;
}
.hanging-indent div.csl-entry {
margin-left:2em;
text-indent:-2em;
}
div.csl-left-margin {
min-width:2em;
float:left;
}
div.csl-right-inline {
margin-left:2em;
padding-left:1em;
}
div.csl-indent {
margin-left: 2em;
}</style>
<script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.5.1/jquery.min.js" integrity="sha512-bLT0Qm9VnAYZDflyKcBaQ2gg0hSYNQrJ8RilYldYQ1FxQYoCLtUjuuRuZo+fjqhx/qtq/1itJ0C2ejDxltZVFg==" crossorigin="anonymous"></script><script src="site_libs/quarto-nav/quarto-nav.js"></script>
<script src="site_libs/quarto-nav/headroom.min.js"></script>
<script src="site_libs/clipboard/clipboard.min.js"></script>
<script src="site_libs/quarto-search/autocomplete.umd.js"></script>
<script src="site_libs/quarto-search/fuse.min.js"></script>
<script src="site_libs/quarto-search/quarto-search.js"></script>
<meta name="quarto:offset" content="./">
<link href="./07-read-write.html" rel="next">
<link href="./05-raster-vector.html" rel="prev">
<link href="./favicon-32x32.png" rel="icon" type="image/png">
<script src="site_libs/quarto-html/quarto.js"></script>
<script src="site_libs/quarto-html/popper.min.js"></script>
<script src="site_libs/quarto-html/tippy.umd.min.js"></script>
<script src="site_libs/quarto-html/anchor.min.js"></script>
<link href="site_libs/quarto-html/tippy.css" rel="stylesheet">
<link href="site_libs/quarto-html/quarto-syntax-highlighting.css" rel="stylesheet" id="quarto-text-highlighting-styles">
<script src="site_libs/bootstrap/bootstrap.min.js"></script>
<link href="site_libs/bootstrap/bootstrap-icons.css" rel="stylesheet">
<link href="site_libs/bootstrap/bootstrap.min.css" rel="stylesheet" id="quarto-bootstrap" data-mode="light">
<script id="quarto-search-options" type="application/json">{
"location": "sidebar",
"copy-button": false,
"collapse-after": 3,
"panel-placement": "start",
"type": "textbox",
"limit": 50,
"keyboard-shortcut": [
"f",
"/",
"s"
],
"show-item-context": false,
"language": {
"search-no-results-text": "No results",
"search-matching-documents-text": "matching documents",
"search-copy-link-title": "Copy link to search",
"search-hide-matches-text": "Hide additional matches",
"search-more-match-text": "more match in this document",
"search-more-matches-text": "more matches in this document",
"search-clear-button-title": "Clear",
"search-text-placeholder": "",
"search-detached-cancel-button-title": "Cancel",
"search-submit-button-title": "Submit",
"search-label": "Search"
}
}</script>
<script async="" src="https://www.googletagmanager.com/gtag/js?id=G-ZEMGTY4VV3"></script>
<script type="text/javascript">
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'G-ZEMGTY4VV3', { 'anonymize_ip': true});
</script>
<script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" integrity="sha512-c3Nl8+7g4LMSTdrm621y7kf9v3SDPnhxLNhcjFJbKECVnmZHTdo+IRO05sNLTH/D3vA6u1X32ehoLC7WFVdheg==" crossorigin="anonymous"></script>
<script type="application/javascript">define('jquery', [],function() {return window.jQuery;})</script>
<script src="https://cdnjs.cloudflare.com/polyfill/v3/polyfill.min.js?features=es6"></script>
<script src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-chtml-full.js" type="text/javascript"></script>
<script type="text/javascript">
const typesetMath = (el) => {
if (window.MathJax) {
// MathJax Typeset
window.MathJax.typeset([el]);
} else if (window.katex) {
// KaTeX Render
var mathElements = el.getElementsByClassName("math");
var macros = [];
for (var i = 0; i < mathElements.length; i++) {
var texText = mathElements[i].firstChild;
if (mathElements[i].tagName == "SPAN") {
window.katex.render(texText.data, mathElements[i], {
displayMode: mathElements[i].classList.contains('display'),
throwOnError: false,
macros: macros,
fleqn: false
});
}
}
}
}
window.Quarto = {
typesetMath
};
</script>
<link rel="stylesheet" href="helpers/mystyle.css">
</head>
<body class="nav-sidebar floating">
<div id="quarto-search-results"></div>
<header id="quarto-header" class="headroom fixed-top">
<nav class="quarto-secondary-nav">
<div class="container-fluid d-flex">
<button type="button" class="quarto-btn-toggle btn" data-bs-toggle="collapse" role="button" data-bs-target=".quarto-sidebar-collapse-item" aria-controls="quarto-sidebar" aria-expanded="false" aria-label="Toggle sidebar navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
<i class="bi bi-layout-text-sidebar-reverse"></i>
</button>
<nav class="quarto-page-breadcrumbs" aria-label="breadcrumb"><ol class="breadcrumb"><li class="breadcrumb-item"><a href="./06-reproj.html"><span class="chapter-number">6</span> <span class="chapter-title">Reprojecting geographic data</span></a></li></ol></nav>
<a class="flex-grow-1" role="navigation" data-bs-toggle="collapse" data-bs-target=".quarto-sidebar-collapse-item" aria-controls="quarto-sidebar" aria-expanded="false" aria-label="Toggle sidebar navigation" onclick="if (window.quartoToggleHeadroom) { window.quartoToggleHeadroom(); }">
</a>
<button type="button" class="btn quarto-search-button" aria-label="Search" onclick="window.quartoOpenSearch();">
<i class="bi bi-search"></i>
</button>
</div>
</nav>
</header>
<!-- content -->
<div id="quarto-content" class="quarto-container page-columns page-rows-contents page-layout-article">
<!-- sidebar -->
<nav id="quarto-sidebar" class="sidebar collapse collapse-horizontal quarto-sidebar-collapse-item sidebar-navigation floating overflow-auto">
<div class="pt-lg-2 mt-2 text-left sidebar-header">
<div class="sidebar-title mb-0 py-0">
<a href="./">Geocomputation with Python</a>
<div class="sidebar-tools-main">
<a href="https://github.com/geocompx/geocompy/" title="Source Code" class="quarto-navigation-tool px-1" aria-label="Source Code"><i class="bi bi-github"></i></a>
<div class="dropdown">
<a href="" title="Share" id="quarto-navigation-tool-dropdown-0" class="quarto-navigation-tool dropdown-toggle px-1" data-bs-toggle="dropdown" aria-expanded="false" role="link" aria-label="Share"><i class="bi bi-share"></i></a>
<ul class="dropdown-menu" aria-labelledby="quarto-navigation-tool-dropdown-0">
<li>
<a class="dropdown-item sidebar-tools-main-item" href="https://twitter.com/intent/tweet?url=|url|">
<i class="bi bi-twitter pe-1"></i>
Twitter
</a>
</li>
<li>
<a class="dropdown-item sidebar-tools-main-item" href="https://www.facebook.com/sharer/sharer.php?u=|url|">
<i class="bi bi-facebook pe-1"></i>
Facebook
</a>
</li>
<li>
<a class="dropdown-item sidebar-tools-main-item" href="https://www.linkedin.com/sharing/share-offsite/?url=|url|">
<i class="bi bi-linkedin pe-1"></i>
LinkedIn
</a>
</li>
</ul>
</div>
</div>
</div>
</div>
<div class="mt-2 flex-shrink-0 align-items-center">
<div class="sidebar-search">
<div id="quarto-search" class="" title="Search"></div>
</div>
</div>
<div class="sidebar-menu-container">
<ul class="list-unstyled mt-1">
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./index.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Welcome</span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./preface.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">Preface</span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./01-spatial-data.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">1</span> <span class="chapter-title">Geographic data in Python</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./02-attribute-operations.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">2</span> <span class="chapter-title">Attribute data operations</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./03-spatial-operations.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">3</span> <span class="chapter-title">Spatial data operations</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./04-geometry-operations.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">4</span> <span class="chapter-title">Geometry operations</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./05-raster-vector.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">5</span> <span class="chapter-title">Raster-vector interactions</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./06-reproj.html" class="sidebar-item-text sidebar-link active">
<span class="menu-text"><span class="chapter-number">6</span> <span class="chapter-title">Reprojecting geographic data</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./07-read-write.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">7</span> <span class="chapter-title">Geographic data I/O</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./08-mapping.html" class="sidebar-item-text sidebar-link">
<span class="menu-text"><span class="chapter-number">8</span> <span class="chapter-title">Making maps with Python</span></span></a>
</div>
</li>
<li class="sidebar-item">
<div class="sidebar-item-container">
<a href="./references.html" class="sidebar-item-text sidebar-link">
<span class="menu-text">References</span></a>
</div>
</li>
</ul>
</div>
</nav>
<div id="quarto-sidebar-glass" class="quarto-sidebar-collapse-item" data-bs-toggle="collapse" data-bs-target=".quarto-sidebar-collapse-item"></div>
<!-- margin-sidebar -->
<div id="quarto-margin-sidebar" class="sidebar margin-sidebar">
<nav id="TOC" role="doc-toc" class="toc-active">
<h2 class="anchored">First Edition</h2>
<ul>
<li><a href="https://geocompx.org/" class="nav-link active" data-scroll-target="https\://geocompx.org/">Visit the geocompx website 🌐</a></li>
<li><a href="https://github.com/geocompx/geocompy/issues" class="nav-link" data-scroll-target="https\://github.com/geocompx/geocompy/issues">Open an issue ❔</a></li>
<li><a href="https://discord.gg/PMztXYgNxp" class="nav-link" data-scroll-target="https\://discord.gg/PMztXYgNxp">Chat on Discord 📣</a></li>
<li><a href="https://donate.stripe.com/4gweWl94Q9E35AQ6oo" class="nav-link" data-scroll-target="https\://donate.stripe.com/4gweWl94Q9E35AQ6oo">Support this project 💸</a></li>
</ul>
<hr>
<h2 id="toc-title">On this page</h2>
<ul>
<li><a href="#prerequisites" id="toc-prerequisites" class="nav-link" data-scroll-target="#prerequisites">Prerequisites</a></li>
<li><a href="#introduction" id="toc-introduction" class="nav-link" data-scroll-target="#introduction"><span class="header-section-number">6.1</span> Introduction</a></li>
<li><a href="#sec-coordinate-reference-systems" id="toc-sec-coordinate-reference-systems" class="nav-link" data-scroll-target="#sec-coordinate-reference-systems"><span class="header-section-number">6.2</span> Coordinate Reference Systems</a></li>
<li><a href="#sec-querying-and-setting-coordinate-systems" id="toc-sec-querying-and-setting-coordinate-systems" class="nav-link" data-scroll-target="#sec-querying-and-setting-coordinate-systems"><span class="header-section-number">6.3</span> Querying and setting coordinate systems</a></li>
<li><a href="#sec-geometry-operations-on-projected-and-unprojected-data" id="toc-sec-geometry-operations-on-projected-and-unprojected-data" class="nav-link" data-scroll-target="#sec-geometry-operations-on-projected-and-unprojected-data"><span class="header-section-number">6.4</span> Geometry operations on projected and unprojected data</a></li>
<li><a href="#sec-when-to-reproject" id="toc-sec-when-to-reproject" class="nav-link" data-scroll-target="#sec-when-to-reproject"><span class="header-section-number">6.5</span> When to reproject?</a></li>
<li><a href="#sec-which-crs-to-use" id="toc-sec-which-crs-to-use" class="nav-link" data-scroll-target="#sec-which-crs-to-use"><span class="header-section-number">6.6</span> Which CRS to use?</a></li>
<li><a href="#sec-reprojecting-vector-geometries" id="toc-sec-reprojecting-vector-geometries" class="nav-link" data-scroll-target="#sec-reprojecting-vector-geometries"><span class="header-section-number">6.7</span> Reprojecting vector geometries</a></li>
<li><a href="#sec-reprojecting-raster-geometries" id="toc-sec-reprojecting-raster-geometries" class="nav-link" data-scroll-target="#sec-reprojecting-raster-geometries"><span class="header-section-number">6.8</span> Reprojecting raster geometries</a></li>
<li><a href="#sec-custom-map-projections" id="toc-sec-custom-map-projections" class="nav-link" data-scroll-target="#sec-custom-map-projections"><span class="header-section-number">6.9</span> Custom map projections</a></li>
</ul>
<div class="toc-actions"><ul><li><a href="https://github.com/geocompx/geocompy/edit/main/06-reproj.qmd" class="toc-action"><i class="bi bi-github"></i>Edit this page</a></li></ul></div></nav>
</div>
<!-- main -->
<main class="content" id="quarto-document-content">
<header id="title-block-header" class="quarto-title-block default">
<div class="quarto-title">
<h1 class="title"><span id="sec-reproj-geo-data" class="quarto-section-identifier"><span class="chapter-number">6</span> <span class="chapter-title">Reprojecting geographic data</span></span></h1>
</div>
<!--
-->
<div class="quarto-title-meta">
</div>
</header>
<section id="prerequisites" class="level2 unnumbered">
<h2 class="unnumbered anchored" data-anchor-id="prerequisites">Prerequisites</h2>
<p>This chapter requires importing the following packages:</p>
<div id="1e03caec" class="cell" data-execution_count="3">
<div class="sourceCode cell-code" id="cb1"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb1-1"><a href="#cb1-1" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> shutil</span>
<span id="cb1-2"><a href="#cb1-2" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> math</span>
<span id="cb1-3"><a href="#cb1-3" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> numpy <span class="im">as</span> np</span>
<span id="cb1-4"><a href="#cb1-4" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> matplotlib.pyplot <span class="im">as</span> plt</span>
<span id="cb1-5"><a href="#cb1-5" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> shapely</span>
<span id="cb1-6"><a href="#cb1-6" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> pyproj</span>
<span id="cb1-7"><a href="#cb1-7" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> geopandas <span class="im">as</span> gpd</span>
<span id="cb1-8"><a href="#cb1-8" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> rasterio</span>
<span id="cb1-9"><a href="#cb1-9" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> rasterio.plot</span>
<span id="cb1-10"><a href="#cb1-10" aria-hidden="true" tabindex="-1"></a><span class="im">import</span> rasterio.warp</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>It also relies on the following data files:</p>
<div id="e46c9320" class="cell" data-execution_count="4">
<div class="sourceCode cell-code" id="cb2"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb2-1"><a href="#cb2-1" aria-hidden="true" tabindex="-1"></a>src_srtm <span class="op">=</span> rasterio.<span class="bu">open</span>(<span class="st">'data/srtm.tif'</span>)</span>
<span id="cb2-2"><a href="#cb2-2" aria-hidden="true" tabindex="-1"></a>src_nlcd <span class="op">=</span> rasterio.<span class="bu">open</span>(<span class="st">'data/nlcd.tif'</span>)</span>
<span id="cb2-3"><a href="#cb2-3" aria-hidden="true" tabindex="-1"></a>zion <span class="op">=</span> gpd.read_file(<span class="st">'data/zion.gpkg'</span>)</span>
<span id="cb2-4"><a href="#cb2-4" aria-hidden="true" tabindex="-1"></a>world <span class="op">=</span> gpd.read_file(<span class="st">'data/world.gpkg'</span>)</span>
<span id="cb2-5"><a href="#cb2-5" aria-hidden="true" tabindex="-1"></a>cycle_hire_osm <span class="op">=</span> gpd.read_file(<span class="st">'data/cycle_hire_osm.gpkg'</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
</section>
<section id="introduction" class="level2" data-number="6.1">
<h2 data-number="6.1" class="anchored" data-anchor-id="introduction"><span class="header-section-number">6.1</span> Introduction</h2>
<p><a href="01-spatial-data.html#sec-coordinate-reference-systems-intro" class="quarto-xref"><span>Section 1.4</span></a> introduced coordinate reference systems (CRSs), with a focus on the two major types: geographic (‘lon/lat’, with units in degrees longitude and latitude) and projected (typically with units of meters from a datum) coordinate systems. This chapter builds on that knowledge and goes further. It demonstrates how to set and transform geographic data from one CRS to another and, furthermore, highlights specific issues that can arise due to ignoring CRSs that you should be aware of, especially if your data is stored with lon/lat coordinates.</p>
<p>It is important to know if your data is in a projected or geographic coordinate system, and the consequences of this for geometry operations. However, if you know the CRS of your data and the consequences for geometry operations (covered in the next section), CRSs should just work behind the scenes: people often suddenly need to learn about CRSs when things go wrong. Having a clearly defined project CRS that all project data is in, plus understanding how and why to use different CRSs, can ensure that things do not go wrong. Furthermore, learning about coordinate systems will deepen your knowledge of geographic datasets and how to use them effectively.</p>
<p>This chapter teaches the fundamentals of CRSs, demonstrates the consequences of using different CRSs (including what can go wrong), and how to ‘reproject’ datasets from one coordinate system to another. In the next section we introduce CRSs in Python, followed by <a href="#sec-querying-and-setting-coordinate-systems" class="quarto-xref"><span>Section 6.3</span></a> which shows how to get and set CRSs associated with spatial objects. <a href="#sec-geometry-operations-on-projected-and-unprojected-data" class="quarto-xref"><span>Section 6.4</span></a> demonstrates the importance of knowing what CRS your data is in with reference to a worked example of creating buffers. We tackle questions of when to reproject and which CRS to use in <a href="#sec-when-to-reproject" class="quarto-xref"><span>Section 6.5</span></a> and <a href="#sec-which-crs-to-use" class="quarto-xref"><span>Section 6.6</span></a>, respectively. Finally, we cover reprojecting vector and raster objects in <a href="#sec-reprojecting-vector-geometries" class="quarto-xref"><span>Section 6.7</span></a> and <a href="#sec-reprojecting-raster-geometries" class="quarto-xref"><span>Section 6.8</span></a> and using custom projections in <a href="#sec-custom-map-projections" class="quarto-xref"><span>Section 6.9</span></a>.</p>
</section>
<section id="sec-coordinate-reference-systems" class="level2" data-number="6.2">
<h2 data-number="6.2" class="anchored" data-anchor-id="sec-coordinate-reference-systems"><span class="header-section-number">6.2</span> Coordinate Reference Systems</h2>
<p>Most modern geographic tools that require CRS conversions, including Python packages and desktop GIS software such as QGIS, interface with PROJ, an open source C++ library that ‘transforms coordinates from one coordinate reference system (CRS) to another’. CRSs can be described in many ways, including the following:</p>
<ul>
<li>Simple, yet potentially ambiguous, statements, such as ‘it’s in lon/lat coordinates’</li>
<li>Formalized, yet now outdated, ‘proj-strings’, such as <code>+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs</code></li>
<li>With an identifying ‘authority:code’ text string, such as <code>EPSG:4326</code></li>
</ul>
<p>Each refers to the same thing: the ‘WGS84’ coordinate system that forms the basis of Global Positioning System (GPS) coordinates and many other datasets. But which one is correct?</p>
<p>The short answer is that the third way to identify CRSs is correct: <code>EPSG:4326</code> is understood by <strong>geopandas</strong> and <strong>rasterio</strong> packages covered in this book, plus many other software projects for working with geographic data including QGIS and PROJ. <code>EPSG:4326</code> is future-proof. Furthermore, although it is machine readable, unlike the proj-string representation <code>EPSG:4326</code> is short, easy to remember and highly ‘findable’ online (searching for <code>EPSG:4326</code> yields a dedicated page on the website epsg.io<a href="#fn1" class="footnote-ref" id="fnref1" role="doc-noteref"><sup>1</sup></a>, for example). The more concise identifier <code>4326</code> is also understood by <strong>geopandas</strong> and <strong>rasterio</strong>.</p>
<p>The longer answer is that none of the three descriptions is sufficient, and more detail is needed for unambiguous CRS handling and transformations: due to the complexity of CRSs, it is not possible to capture all relevant information about them in such short text strings. For this reason, the Open Geospatial Consortium (OGC, which also developed the Simple Features specification that the <strong>geopandas</strong> package implements) developed an open standard format for describing CRSs that is called WKT (Well Known Text). This is detailed in a 100+ page document that ‘defines the structure and content of a text string implementation of the abstract model for coordinate reference systems described in ISO 19111:2019’ <span class="citation" data-cites="opengeospatialconsortium_wellknown_2019">(<a href="references.html#ref-opengeospatialconsortium_wellknown_2019" role="doc-biblioref">Open Geospatial Consortium 2019</a>)</span>. The WKT representation of the WGS84 CRS, which has the identifier <code>EPSG:4326</code> is as follows.</p>
<div id="a667c876" class="cell" data-execution_count="5">
<div class="sourceCode cell-code" id="cb3"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb3-1"><a href="#cb3-1" aria-hidden="true" tabindex="-1"></a>crs <span class="op">=</span> pyproj.CRS.from_string(<span class="st">'EPSG:4326'</span>) <span class="co"># or '.from_epsg(4326)'</span></span>
<span id="cb3-2"><a href="#cb3-2" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(crs.to_wkt(pretty<span class="op">=</span><span class="va">True</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]</code></pre>
</div>
</div>
<p>The output of the command shows how the CRS identifier (also known as a Spatial Reference Identifier, or SRID) works: it is simply a look-up, providing a unique identifier associated with a more complete WKT representation of the CRS. This raises the question: what happens if there is a mismatch between the identifier and the longer WKT representation of a CRS? On this point Open Geospatial Consortium <span class="citation" data-cites="opengeospatialconsortium_wellknown_2019">(<a href="references.html#ref-opengeospatialconsortium_wellknown_2019" role="doc-biblioref">Open Geospatial Consortium 2019</a>)</span> is clear, the verbose WKT representation takes precedence over the identifier:</p>
<blockquote class="blockquote">
<p>Should any attributes or values given in the cited identifier be in conflict with attributes or values given explicitly in the WKT description, the WKT values shall prevail.</p>
</blockquote>
<p>The convention of referring to CRSs identifiers in the form <code>AUTHORITY:CODE</code> allows a wide range of formally defined coordinate systems to be referred to. The most commonly used authority in CRS identifiers is EPSG, an acronym for the European Petroleum Survey Group which published a standardized list of CRSs. Other authorities can be used in CRS identifiers. <code>ESRI:54030</code>, for example, refers to ESRI’s implementation of the Robinson projection, which has the following WKT string.</p>
<div id="9fb95ce8" class="cell" data-execution_count="6">
<div class="sourceCode cell-code" id="cb5"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb5-1"><a href="#cb5-1" aria-hidden="true" tabindex="-1"></a>crs <span class="op">=</span> pyproj.CRS.from_string(<span class="st">'ESRI:54030'</span>)</span>
<span id="cb5-2"><a href="#cb5-2" aria-hidden="true" tabindex="-1"></a><span class="bu">print</span>(crs.to_wkt(pretty<span class="op">=</span><span class="va">True</span>))</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stdout">
<pre><code>PROJCRS["World_Robinson",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["World_Robinson",
METHOD["Robinson"],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Not known."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["ESRI",54030]]</code></pre>
</div>
</div>
<p>WKT strings are exhaustive, detailed, and precise, allowing for unambiguous CRSs storage and transformations. They contain all relevant information about any given CRS, including its datum and ellipsoid, prime meridian, projection, and units.</p>
<p>Recent PROJ versions (6+) still allow use of proj-strings to define coordinate operations, but some proj-string keys (<code>+nadgrids</code>, <code>+towgs84</code>, <code>+k</code>, <code>+init=epsg:</code>) are either no longer supported or are discouraged. Additionally, only three datums (i.e., WGS84, NAD83, and NAD27) can be directly set in proj-string. Longer explanations of the evolution of CRS definitions and the PROJ library can be found in <span class="citation" data-cites="bivand_progress_2021">Bivand (<a href="references.html#ref-bivand_progress_2021" role="doc-biblioref">2021</a>)</span>, Chapter 2 of <span class="citation" data-cites="pebesma_spatial_2022">Pebesma and Bivand (<a href="references.html#ref-pebesma_spatial_2022" role="doc-biblioref">2022</a>)</span>, and a blog post by Floris Vanderhaeghe<a href="#fn2" class="footnote-ref" id="fnref2" role="doc-noteref"><sup>2</sup></a>.</p>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>As outlined in the PROJ documentation, there are different versions of the WKT CRS format including WKT1 and two variants of WKT2, the latter of which (WKT2, 2018 specification) corresponds to the ISO 19111:2019 <span class="citation" data-cites="opengeospatialconsortium_wellknown_2019">(<a href="references.html#ref-opengeospatialconsortium_wellknown_2019" role="doc-biblioref">Open Geospatial Consortium 2019</a>)</span>.</p>
</div>
</div>
</section>
<section id="sec-querying-and-setting-coordinate-systems" class="level2" data-number="6.3">
<h2 data-number="6.3" class="anchored" data-anchor-id="sec-querying-and-setting-coordinate-systems"><span class="header-section-number">6.3</span> Querying and setting coordinate systems</h2>
<p>Let’s see how CRSs are stored in Python spatial objects and how they can be queried and set. First, we will look at getting and setting CRSs in vector geographic data objects. Consider the <code>GeoDataFrame</code> object named <code>world</code>, imported from a file <code>world.gpkg</code> that represents countries worldwide. Its CRS can be retrieved using the <code>.crs</code> property.</p>
<div id="42257a20" class="cell" data-execution_count="7">
<div class="sourceCode cell-code" id="cb7"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb7-1"><a href="#cb7-1" aria-hidden="true" tabindex="-1"></a>world.crs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="7">
<pre><code><Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich</code></pre>
</div>
</div>
<p>The output specifies the following pieces of information:</p>
<ol type="1">
<li>The CRS type (<code>Geographic 2D CRS</code>) and SRID code (<code>EPSG:4326</code>)</li>
<li>The CRS name (<code>WGS 84</code>)</li>
<li>The axes (<code>latitude</code>, <code>longitude</code>) and their units (<code>degree</code>)</li>
<li>The applicable area name (<code>World</code>) and bounding box (<code>(-180.0, -90.0, 180.0, 90.0)</code>)</li>
<li>The datum (<code>WGS 84</code>)</li>
</ol>
<p>The WKT representation, which is internally used when saving the object to a file or doing any coordinate operations, can be extracted using <code>.crs.to_wkt()</code> as shown above (<a href="#sec-coordinate-reference-systems" class="quarto-xref"><span>Section 6.2</span></a>). We can also see that the <code>world</code> object has the WGS84 ellipsoid, the latitude and longitude axis order, and uses the Greenwich prime meridian. We also have the suitable area specification for the use of this CRS, and CRS identifier: <code>EPSG:4326</code>.</p>
<p>The CRS specification object, such as <code>world.crs</code>, has several useful properties and methods to explicitly retrieve information about the used CRS. For example, we can check whether the CRS is geographic with the <code>.is_geographic</code> property.</p>
<div id="36cae9b4" class="cell" data-execution_count="8">
<div class="sourceCode cell-code" id="cb9"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb9-1"><a href="#cb9-1" aria-hidden="true" tabindex="-1"></a>world.crs.is_geographic</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="8">
<pre><code>True</code></pre>
</div>
</div>
<p>CRS units of both axes (typically identical) can be retrieved with the <code>.axis_info</code> property.</p>
<div id="a38a46c2" class="cell" data-execution_count="9">
<div class="sourceCode cell-code" id="cb11"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb11-1"><a href="#cb11-1" aria-hidden="true" tabindex="-1"></a>world.crs.axis_info[<span class="dv">0</span>].unit_name, world.crs.axis_info[<span class="dv">1</span>].unit_name</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="9">
<pre><code>('degree', 'degree')</code></pre>
</div>
</div>
<p><code>AUTHORITY</code> and <code>CODE</code> strings may be obtained with the <code>.to_authority()</code> method.</p>
<div id="c44cd61d" class="cell" data-execution_count="10">
<div class="sourceCode cell-code" id="cb13"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb13-1"><a href="#cb13-1" aria-hidden="true" tabindex="-1"></a>world.crs.to_authority()</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="10">
<pre><code>('EPSG', '4326')</code></pre>
</div>
</div>
<p>In cases when a coordinate reference system (CRS) is missing or the wrong CRS is set, the <code>.set_crs</code> method can be used on a <code>GeoSeries</code> or a <code>GeoDataFrame</code> to set it. The CRS can be specified using an EPSG code as the first argument. In case the object already has a different CRS definition, we must also specify <code>allow_override=True</code> to replace it (otherwise we get an error). In the first example we set the <code>EPSG:4326</code> CRS, which has no effect because <code>world</code> already has that exact CRS definition, while the second example replaces the existing CRS with a new definition of <code>EPSG:3857</code>.</p>
<div id="65a2697d" class="cell" data-execution_count="11">
<div class="sourceCode cell-code" id="cb15"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb15-1"><a href="#cb15-1" aria-hidden="true" tabindex="-1"></a>world2 <span class="op">=</span> world.set_crs(<span class="dv">4326</span>)</span>
<span id="cb15-2"><a href="#cb15-2" aria-hidden="true" tabindex="-1"></a>world3 <span class="op">=</span> world.set_crs(<span class="dv">3857</span>, allow_override<span class="op">=</span><span class="va">True</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>The provided number is interpreted as an <code>EPSG</code> code. We can also use strings, as in <code>'EPSG:4326'</code>, which is useful to make the code more clear and when using other authorities than <code>EPSG</code>.</p>
<div id="9c42fbc9" class="cell" data-execution_count="12">
<div class="sourceCode cell-code" id="cb16"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb16-1"><a href="#cb16-1" aria-hidden="true" tabindex="-1"></a>world4 <span class="op">=</span> world.set_crs(<span class="st">'ESRI:54009'</span>, allow_override<span class="op">=</span><span class="va">True</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>In <strong>rasterio</strong>, the CRS information is stored as part of a raster file connection metadata (<a href="01-spatial-data.html#sec-using-rasterio" class="quarto-xref"><span>Section 1.3.1</span></a>). Replacing the CRS definition for a <strong>rasterio</strong> file connection is typically not necessary, because it is not considered in any operation; only the transformation matrix and coordinates are. One exception is when writing the raster, in which case we need to construct the metadata of the raster file to be written, and therein specify the CRS anyway (<a href="01-spatial-data.html#sec-raster-from-scratch" class="quarto-xref"><span>Section 1.3.2</span></a>). However, if we, for some reason, need to change the CRS definition in the file connection metadata, we can do that when opening the file in <code>r+</code> (reading and writing) mode. To demonstrate, we will create a copy of the <code>nlcd.tif</code> file, named <code>nlcd_modified_crs.tif</code>,</p>
<div id="74c32c37" class="cell" data-execution_count="13">
<div class="sourceCode cell-code" id="cb17"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb17-1"><a href="#cb17-1" aria-hidden="true" tabindex="-1"></a>shutil.copy(<span class="st">'data/nlcd.tif'</span>, <span class="st">'output/nlcd_modified_crs.tif'</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="13">
<pre><code>'output/nlcd_modified_crs.tif'</code></pre>
</div>
</div>
<p>and examine its existing CRS.</p>
<div id="e6709032" class="cell" data-execution_count="14">
<div class="sourceCode cell-code" id="cb19"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb19-1"><a href="#cb19-1" aria-hidden="true" tabindex="-1"></a>src_nlcd2 <span class="op">=</span> rasterio.<span class="bu">open</span>(<span class="st">'output/nlcd_modified_crs.tif'</span>, <span class="st">'r+'</span>)</span>
<span id="cb19-2"><a href="#cb19-2" aria-hidden="true" tabindex="-1"></a>src_nlcd2.crs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="14">
<pre><code>CRS.from_epsg(26912)</code></pre>
</div>
</div>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>The <code>rasterio.open</code> function <code>mode</code>s generally follows Python’s standard file connection modes, with possible arguments being <code>'r'</code> (read), <code>'w'</code> (write), <code>'r+'</code> (read/write), and <code>'w+'</code> (write/read) (the <code>'a'</code> ‘append’ mode is irrelevant for raster files). In the book, and in general, the most commonly used modes are <code>'r'</code> (read) and <code>'w'</code> (write). <code>'r+'</code>, used in the last example, means ‘read/write’. Unlike with <code>'w'</code>, <code>'r+'</code> does not delete the existing content on open, making <code>'r+'</code> suitable for making changes in an existing file (such as here, replacing the CRS).</p>
</div>
</div>
<p>To replace the definition with a new one, such as <code>EPSG:3857</code>, we can use the <code>.crs</code> method, as shown below.</p>
<div id="5a480e86" class="cell" data-execution_count="15">
<div class="sourceCode cell-code" id="cb21"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb21-1"><a href="#cb21-1" aria-hidden="true" tabindex="-1"></a>src_nlcd2.crs <span class="op">=</span> <span class="dv">3857</span></span>
<span id="cb21-2"><a href="#cb21-2" aria-hidden="true" tabindex="-1"></a>src_nlcd2.close()</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Next, examining the file connection demonstrates that the CRS was indeed changed.</p>
<div id="f3b67a33" class="cell" data-execution_count="16">
<div class="sourceCode cell-code" id="cb22"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb22-1"><a href="#cb22-1" aria-hidden="true" tabindex="-1"></a>rasterio.<span class="bu">open</span>(<span class="st">'output/nlcd_modified_crs.tif'</span>).crs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="16">
<pre><code>CRS.from_epsg(3857)</code></pre>
</div>
</div>
<p>Importantly, the <code>.set_crs</code> (for vector layers) or the assignment to <code>.crs</code> (for rasters), as shown above, do not alter coordinates’ values or geometries. Their role is only to set a metadata information about the object CRS. Consequently, the objects we created, <code>world3</code>, <code>world4</code>, and <code>src_nlcd2</code> are ‘incorrect’, in the sense that the geometries are in fact given in a different CRS than specified in the associated CRS definition.</p>
<p>In some cases, the CRS of a geographic object is unknown, as is the case in the London dataset created in the code chunk below, building on the example of London introduced in <a href="01-spatial-data.html#sec-vector-layer-from-scratch" class="quarto-xref"><span>Section 1.2.6</span></a>.</p>
<div id="9513ff56" class="cell" data-execution_count="17">
<div class="sourceCode cell-code" id="cb24"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb24-1"><a href="#cb24-1" aria-hidden="true" tabindex="-1"></a>lnd_point <span class="op">=</span> shapely.Point(<span class="op">-</span><span class="fl">0.1</span>, <span class="fl">51.5</span>)</span>
<span id="cb24-2"><a href="#cb24-2" aria-hidden="true" tabindex="-1"></a>lnd_geom <span class="op">=</span> gpd.GeoSeries([lnd_point])</span>
<span id="cb24-3"><a href="#cb24-3" aria-hidden="true" tabindex="-1"></a>lnd_layer <span class="op">=</span> gpd.GeoDataFrame({<span class="st">'geometry'</span>: lnd_geom})</span>
<span id="cb24-4"><a href="#cb24-4" aria-hidden="true" tabindex="-1"></a>lnd_layer</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="17">
<div>
<table class="dataframe caption-top table table-sm table-striped small" data-quarto-postprocess="true" data-border="1">
<thead>
<tr class="header">
<th data-quarto-table-cell-role="th"></th>
<th data-quarto-table-cell-role="th">geometry</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td data-quarto-table-cell-role="th">0</td>
<td>POINT (-0.1 51.5)</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
<p>Querying the <code>.crs</code> of such a layer returns <code>None</code>, therefore nothing is printed.</p>
<div id="083a8cc9" class="cell" data-execution_count="18">
<div class="sourceCode cell-code" id="cb25"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb25-1"><a href="#cb25-1" aria-hidden="true" tabindex="-1"></a>lnd_layer.crs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>This implies that <strong>geopandas</strong> does not know what the CRS is and is unwilling to guess. Unless a CRS is manually specified or is loaded from a source that has CRS metadata, <strong>geopandas</strong> does not make any explicit assumptions about which coordinate systems, other than to say ‘I don’t know’. This behavior makes sense given the diversity of available CRSs but differs from some approaches, such as the GeoJSON file format specification, which makes the simplifying assumption that all coordinates have a lon/lat CRS: <code>EPSG:4326</code>.</p>
<p>A CRS can be added to <code>GeoSeries</code> or <code>GeoDataFrame</code> objects using the <code>.set_crs</code> method, as mentioned above.</p>
<div id="7712c8dc" class="cell" data-execution_count="19">
<div class="sourceCode cell-code" id="cb26"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb26-1"><a href="#cb26-1" aria-hidden="true" tabindex="-1"></a>lnd_layer <span class="op">=</span> lnd_layer.set_crs(<span class="dv">4326</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>When working with <strong>geopandas</strong> and <strong>rasterio</strong>, datasets without a specified CRS are not an issue in most workflows, since only the coordinates are considered. It is up to the user to make sure that, when working with more than one layer, all of the coordinates are given in the same CRS (whether specified or not). When exporting the results, though, it is important to keep the CRS definition in place, because other software typically <em>do</em> use, and require, the CRS definition in calculations. It should also be mentioned that, in some cases the CRS specification is left unspecified on purpose, for example when working with layers in arbitrary or non-geographic space (simulations, internal building plans, analysis of plot-scale ecological patterns, etc.).</p>
</section>
<section id="sec-geometry-operations-on-projected-and-unprojected-data" class="level2" data-number="6.4">
<h2 data-number="6.4" class="anchored" data-anchor-id="sec-geometry-operations-on-projected-and-unprojected-data"><span class="header-section-number">6.4</span> Geometry operations on projected and unprojected data</h2>
<p>The <strong>geopandas</strong> package, through its dependency <strong>shapely</strong>, assumes planar geometry and works with distance/area values assumed to be in CRS units. In fact, the CRS definition is typically ignored, and the respective functions (such as in plotting and distance calculations) are applied on the ‘bare’ <strong>shapely</strong> geometries. Accordingly, it is crucial to make sure that:</p>
<ul>
<li>Geometric calculations are only applied in projected CRS</li>
<li>If there is more than one layer involved—all layers have to be in the same (projected) CRS</li>
<li>Distance and area values, are passed, and returned, in CRS units</li>
</ul>
<p>For example, to calculate a buffer of 100 <span class="math inline">\(km\)</span> around London, we need to work with a layer representing London in a projected CRS (e.g., <code>EPSG:27700</code>) and pass the distance value in the CRS units (e.g., <code>100000</code> <span class="math inline">\(m\)</span>).</p>
<p>In the following code chunk we create, from scratch, a point layer <code>lnd_layer_proj</code> with a point representing London (compare it to <code>lnd_layer</code>, in a geographical CRS which we created above, see <a href="#sec-querying-and-setting-coordinate-systems" class="quarto-xref"><span>Section 6.3</span></a>).</p>
<div id="39a11f68" class="cell" data-execution_count="20">
<div class="sourceCode cell-code" id="cb27"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb27-1"><a href="#cb27-1" aria-hidden="true" tabindex="-1"></a>lnd_point_proj <span class="op">=</span> shapely.Point(<span class="dv">530000</span>, <span class="dv">180000</span>)</span>
<span id="cb27-2"><a href="#cb27-2" aria-hidden="true" tabindex="-1"></a>lnd_geom_proj <span class="op">=</span> gpd.GeoSeries([lnd_point_proj], crs<span class="op">=</span><span class="dv">27700</span>)</span>
<span id="cb27-3"><a href="#cb27-3" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj <span class="op">=</span> gpd.GeoDataFrame({<span class="st">'geometry'</span>: lnd_geom_proj})</span>
<span id="cb27-4"><a href="#cb27-4" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="20">
<div>
<table class="dataframe caption-top table table-sm table-striped small" data-quarto-postprocess="true" data-border="1">
<thead>
<tr class="header">
<th data-quarto-table-cell-role="th"></th>
<th data-quarto-table-cell-role="th">geometry</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td data-quarto-table-cell-role="th">0</td>
<td>POINT (530000 180000)</td>
</tr>
</tbody>
</table>
</div>
</div>
</div>
<p>Now, we can use the <code>.buffer</code> method (<a href="04-geometry-operations.html#sec-buffers" class="quarto-xref"><span>Section 4.2.3</span></a>) to calculate the buffer of 100 <span class="math inline">\(km\)</span> around London.</p>
<div id="9a8430eb" class="cell" data-execution_count="21">
<div class="sourceCode cell-code" id="cb28"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb28-1"><a href="#cb28-1" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj_buff <span class="op">=</span> lnd_layer_proj.<span class="bu">buffer</span>(<span class="dv">100000</span>)</span>
<span id="cb28-2"><a href="#cb28-2" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj_buff</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="21">
<pre><code>0 POLYGON ((630000 180000, 629518...
dtype: geometry</code></pre>
</div>
</div>
<p>The resulting buffer is shown in the left panel of <a href="#fig-reprojection-geo-proj" class="quarto-xref">Figure <span>6.1</span></a>.</p>
<p>Calculating a 100-<span class="math inline">\(km\)</span> buffer directly for <code>lnd_layer</code>, which is in a geographical CRS, is impossible. Since the <code>lnd_layer</code> is in decimal degrees, the closest thing to a 100-<span class="math inline">\(km\)</span> buffer would be to use a distance of 1 degree, which is roughly equivalent to 100 <span class="math inline">\(km\)</span> (1 degree is about 111 <span class="math inline">\(km\)</span> at the equator):</p>
<div id="2f8a430b" class="cell" data-execution_count="22">
<div class="sourceCode cell-code" id="cb30"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb30-1"><a href="#cb30-1" aria-hidden="true" tabindex="-1"></a>lnd_layer_buff <span class="op">=</span> lnd_layer.<span class="bu">buffer</span>(<span class="dv">1</span>)</span>
<span id="cb30-2"><a href="#cb30-2" aria-hidden="true" tabindex="-1"></a>lnd_layer_buff</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stderr">
<pre><code>/tmp/ipykernel_363/855451079.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
lnd_layer_buff = lnd_layer.buffer(1)</code></pre>
</div>
<div class="cell-output cell-output-display" data-execution_count="22">
<pre><code>0 POLYGON ((0.9 51.5, 0.89518 51....
dtype: geometry</code></pre>
</div>
</div>
<p>However, this is incorrect, as told by the warning message and shown in the right panel of <a href="#fig-reprojection-geo-proj" class="quarto-xref">Figure <span>6.1</span></a>. The association between degrees and true distance varies over the surface of the earth and we cannot assume it is fixed.</p>
<div class="sourceCode cell-code" id="cb33"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb33-1"><a href="#cb33-1" aria-hidden="true" tabindex="-1"></a>uk <span class="op">=</span> world[world[<span class="st">'name_long'</span>] <span class="op">==</span> <span class="st">'United Kingdom'</span>]</span>
<span id="cb33-2"><a href="#cb33-2" aria-hidden="true" tabindex="-1"></a>uk_proj <span class="op">=</span> uk.to_crs(<span class="dv">27700</span>)</span>
<span id="cb33-3"><a href="#cb33-3" aria-hidden="true" tabindex="-1"></a><span class="co"># Around projected point</span></span>
<span id="cb33-4"><a href="#cb33-4" aria-hidden="true" tabindex="-1"></a>base <span class="op">=</span> uk_proj.plot(color<span class="op">=</span><span class="st">'none'</span>, edgecolor<span class="op">=</span><span class="st">'darkgrey'</span>, linewidth<span class="op">=</span><span class="fl">0.5</span>)</span>
<span id="cb33-5"><a href="#cb33-5" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj_buff.plot(color<span class="op">=</span><span class="st">'grey'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, ax<span class="op">=</span>base)</span>
<span id="cb33-6"><a href="#cb33-6" aria-hidden="true" tabindex="-1"></a>lnd_layer_proj.plot(color<span class="op">=</span><span class="st">'red'</span>, ax<span class="op">=</span>base)<span class="op">;</span></span>
<span id="cb33-7"><a href="#cb33-7" aria-hidden="true" tabindex="-1"></a><span class="co"># Around point in lon/lat</span></span>
<span id="cb33-8"><a href="#cb33-8" aria-hidden="true" tabindex="-1"></a>base <span class="op">=</span> uk.plot(color<span class="op">=</span><span class="st">'none'</span>, edgecolor<span class="op">=</span><span class="st">'darkgrey'</span>, linewidth<span class="op">=</span><span class="fl">0.5</span>)</span>
<span id="cb33-9"><a href="#cb33-9" aria-hidden="true" tabindex="-1"></a>lnd_layer_buff.plot(color<span class="op">=</span><span class="st">'grey'</span>, edgecolor<span class="op">=</span><span class="st">'black'</span>, alpha<span class="op">=</span><span class="fl">0.5</span>, ax<span class="op">=</span>base)</span>
<span id="cb33-10"><a href="#cb33-10" aria-hidden="true" tabindex="-1"></a>lnd_layer.plot(color<span class="op">=</span><span class="st">'red'</span>, ax<span class="op">=</span>base)<span class="op">;</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div id="fig-reprojection-geo-proj" class="quarto-layout-panel">
<figure class="quarto-float quarto-float-fig figure">
<div aria-describedby="fig-reprojection-geo-proj-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
<div class="quarto-layout-row">
<div class="cell-output cell-output-display quarto-layout-cell-subref quarto-layout-cell" data-ref-parent="fig-reprojection-geo-proj" style="flex-basis: 50.0%;justify-content: flex-start;">
<div id="fig-reprojection-geo-proj-1" class="quarto-float quarto-figure quarto-figure-center anchored">
<figure class="quarto-float quarto-subfloat-fig figure">
<div aria-describedby="fig-reprojection-geo-proj-1-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
<img src="06-reproj_files/figure-html/fig-reprojection-geo-proj-output-1.png" data-ref-parent="fig-reprojection-geo-proj" width="285" height="425" class="figure-img">
</div>
<figcaption class="quarto-float-caption-bottom quarto-subfloat-caption quarto-subfloat-fig" id="fig-reprojection-geo-proj-1-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
(a) Around a projected point and distance of 100 <span class="math inline">\(km\)</span>
</figcaption>
</figure>
</div>
</div>
<div class="cell-output cell-output-display quarto-layout-cell-subref quarto-layout-cell" data-ref-parent="fig-reprojection-geo-proj" style="flex-basis: 50.0%;justify-content: flex-start;">
<div id="fig-reprojection-geo-proj-2" class="quarto-float quarto-figure quarto-figure-center anchored">
<figure class="quarto-float quarto-subfloat-fig figure">
<div aria-describedby="fig-reprojection-geo-proj-2-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
<img src="06-reproj_files/figure-html/fig-reprojection-geo-proj-output-2.png" data-ref-parent="fig-reprojection-geo-proj" width="292" height="411" class="figure-img">
</div>
<figcaption class="quarto-float-caption-bottom quarto-subfloat-caption quarto-subfloat-fig" id="fig-reprojection-geo-proj-2-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
(b) Around a point in lon/lat using distance of 1 degree (incorrectly approximating 100 <span class="math inline">\(km\)</span>)
</figcaption>
</figure>
</div>
</div>
</div>
</div>
<figcaption class="quarto-float-caption-bottom quarto-float-caption quarto-float-fig" id="fig-reprojection-geo-proj-caption-0ceaefa1-69ba-4598-a22c-09a6ac19f8ca">
Figure 6.1: Buffers around London
</figcaption>
</figure>
</div>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>The distance between two lines of longitude, called meridians, is around 111 <span class="math inline">\(km\)</span> at the equator (execute <code>import geopy.distance;geopy.distance.geodesic((0,0),(0,1))</code> to find the precise distance). This shrinks to zero at the poles. At the latitude of London, for example, meridians are less than 70 <span class="math inline">\(km\)</span> apart (challenge: execute code that verifies this). Lines of latitude, by contrast, are equidistant from each other irrespective of latitude: they are always around 111 <span class="math inline">\(km\)</span> apart, including at the equator and near the poles.</p>
</div>
</div>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>The <strong>spherely</strong><a href="#fn3" class="footnote-ref" id="fnref3" role="doc-noteref"><sup>3</sup></a> package, in early stages of development at the time of writing, is aimed at providing a spherical-geometry counterpart to <strong>shapely</strong>, so that true distances (in <span class="math inline">\(m\)</span>) and areas (in <span class="math inline">\(m^2\)</span>) can be directly calculated on geometries in geographic CRS.</p>
</div>
</div>
</section>
<section id="sec-when-to-reproject" class="level2" data-number="6.5">
<h2 data-number="6.5" class="anchored" data-anchor-id="sec-when-to-reproject"><span class="header-section-number">6.5</span> When to reproject?</h2>
<p>The previous section showed how to set the CRS manually, with an expression such as <code>lnd_layer.set_crs(4326)</code>. In real-world applications, however, CRSs are usually set automatically when data is read-in. Thus, in many projects the main CRS-related task is to transform objects, from one CRS into another. But when should data be transformed? And into which CRS? There are no clear-cut answers to these questions and CRS selection always involves trade-offs <span class="citation" data-cites="maling_coordinate_1992">(<a href="references.html#ref-maling_coordinate_1992" role="doc-biblioref">Maling 1992</a>)</span>. However, there are some general principles provided in this section that can help you decide.</p>
<p>First, it’s worth considering when to transform. In some cases, transformation to a geographic CRS is essential, such as when publishing data online (for example, a Leaflet-based map using Python package <strong>folium</strong>). Another case is when two objects with different CRSs must be compared or combined, as shown when we try to find the distance between two objects with different CRSs.</p>
<div id="fd8df35b" class="cell" data-execution_count="24">
<div class="sourceCode cell-code" id="cb34"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb34-1"><a href="#cb34-1" aria-hidden="true" tabindex="-1"></a>lnd_layer.distance(lnd_layer_proj)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-stderr">
<pre><code>/tmp/ipykernel_363/2145313019.py:1: UserWarning: Geometry is in a geographic CRS. Results from 'distance' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
lnd_layer.distance(lnd_layer_proj)
/tmp/ipykernel_363/2145313019.py:1: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.
Left CRS: EPSG:4326
Right CRS: EPSG:27700
lnd_layer.distance(lnd_layer_proj)</code></pre>
</div>
<div class="cell-output cell-output-display" data-execution_count="24">
<pre><code>0 559715.614087
dtype: float64</code></pre>
</div>
</div>
<p>Here, we got a meaningless distance value of <code>559715</code>, and a warning.</p>
<p>To make the <code>lnd_layer</code> and <code>lnd_layer_proj</code> objects geographically comparable, one of them must be transformed into the CRS of the other. But which CRS to use? The answer depends on context: many projects, especially those involving web mapping, require outputs in <code>EPSG:4326</code>, in which case it is worth transforming the projected object. If, however, the project requires geometric calculations, implying planar geometry, e.g., calculating buffers (<a href="#sec-geometry-operations-on-projected-and-unprojected-data" class="quarto-xref"><span>Section 6.4</span></a>), it is necessary to transform data with a geographic CRS into an equivalent object with a projected CRS, such as the British National Grid (<code>EPSG:27700</code>). That is the subject of <a href="#sec-which-crs-to-use" class="quarto-xref"><span>Section 6.6</span></a>.</p>
</section>
<section id="sec-which-crs-to-use" class="level2" data-number="6.6">
<h2 data-number="6.6" class="anchored" data-anchor-id="sec-which-crs-to-use"><span class="header-section-number">6.6</span> Which CRS to use?</h2>
<p>The question of which CRS is tricky, and there is rarely a ‘right’ answer: ‘There exist no all-purpose projections, all involve distortion when far from the center of the specified frame’ <span class="citation" data-cites="bivand_applied_2013">(<a href="references.html#ref-bivand_applied_2013" role="doc-biblioref">Bivand, Pebesma, and Gómez-Rubio 2013</a>)</span>. Additionally, you should not be attached just to one projection for every task. It is possible to use one projection for some part of the analysis, another projection for a different part, and even some other for visualization. Always try to pick the CRS that serves your goal best!</p>
<p>When selecting <em>geographic</em> CRSs, the answer is often WGS84. It is used not only for web mapping, but also because GPS datasets and thousands of raster and vector datasets are provided in this CRS by default. WGS84 is the most common CRS in the world, so it is worth knowing its EPSG code: <code>4326</code>. This ‘magic number’ can be used to convert objects with unusual projected CRSs into something that is widely understood.</p>
<p>What about when a <em>projected</em> CRS is required? In some cases, it is not something that we are free to decide: ‘often the choice of projection is made by a public mapping agency’ <span class="citation" data-cites="bivand_applied_2013">(<a href="references.html#ref-bivand_applied_2013" role="doc-biblioref">Bivand, Pebesma, and Gómez-Rubio 2013</a>)</span>. This means that when working with local data sources, it is likely preferable to work with the CRS in which the data was provided, to ensure compatibility, even if the official CRS is not the most accurate. The example of London was easy to answer because the British National Grid (with its associated EPSG code <code>27700</code>) is well known, and the original dataset (<code>lnd_layer</code>) already had that CRS.</p>
<p>A commonly used default is Universal Transverse Mercator (UTM), a set of CRSs that divide the Earth into 60 longitudinal wedges and 20 latitudinal segments. The transverse Mercator projection used by UTM CRSs is conformal but distorts areas and distances with increasing severity with distance from the center of the UTM zone. Documentation from the GIS software Manifold therefore suggests restricting the longitudinal extent of projects using UTM zones to 6 degrees from the central meridian<a href="#fn4" class="footnote-ref" id="fnref4" role="doc-noteref"><sup>4</sup></a>. Therefore, we recommend using UTM only when your focus is on preserving angles for a relatively small area!</p>
<p>Almost every place on Earth has a UTM code, such as <code>'60H'</code> which refers, among others, to northern New Zealand. UTM EPSG codes run sequentially from <code>32601</code> to <code>32660</code> for northern hemisphere locations and from <code>32701</code> to <code>32760</code> for southern hemisphere locations.</p>
<p>To show how the system works, let’s create a function, <code>lonlat2UTM</code> to calculate the EPSG code associated with any point on the planet.</p>
<div id="d843a599" class="cell" data-execution_count="25">
<div class="sourceCode cell-code" id="cb37"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb37-1"><a href="#cb37-1" aria-hidden="true" tabindex="-1"></a><span class="kw">def</span> lonlat2UTM(lon, lat):</span>
<span id="cb37-2"><a href="#cb37-2" aria-hidden="true" tabindex="-1"></a> utm <span class="op">=</span> (math.floor((lon <span class="op">+</span> <span class="dv">180</span>) <span class="op">/</span> <span class="dv">6</span>) <span class="op">%</span> <span class="dv">60</span>) <span class="op">+</span> <span class="dv">1</span></span>
<span id="cb37-3"><a href="#cb37-3" aria-hidden="true" tabindex="-1"></a> <span class="cf">if</span> lat <span class="op">></span> <span class="dv">0</span>:</span>
<span id="cb37-4"><a href="#cb37-4" aria-hidden="true" tabindex="-1"></a> utm <span class="op">+=</span> <span class="dv">32600</span></span>
<span id="cb37-5"><a href="#cb37-5" aria-hidden="true" tabindex="-1"></a> <span class="cf">else</span>:</span>
<span id="cb37-6"><a href="#cb37-6" aria-hidden="true" tabindex="-1"></a> utm <span class="op">+=</span> <span class="dv">32700</span></span>
<span id="cb37-7"><a href="#cb37-7" aria-hidden="true" tabindex="-1"></a> <span class="cf">return</span> utm</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>The following command uses this function to identify the UTM zone and associated EPSG code for Auckland.</p>
<div id="f031484e" class="cell" data-execution_count="26">
<div class="sourceCode cell-code" id="cb38"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb38-1"><a href="#cb38-1" aria-hidden="true" tabindex="-1"></a>lonlat2UTM(<span class="fl">174.7</span>, <span class="op">-</span><span class="fl">36.9</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="26">
<pre><code>32760</code></pre>
</div>
</div>
<p>Here is another example for London (where we ‘unpack’ the coordinates of the 1<sup>st</sup> geometry in <code>lnd_layer</code> into the <code>lonlat2UTM</code> function arguments).</p>
<div id="c739c97a" class="cell" data-execution_count="27">
<div class="sourceCode cell-code" id="cb40"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb40-1"><a href="#cb40-1" aria-hidden="true" tabindex="-1"></a>lonlat2UTM(<span class="op">*</span>lnd_layer.geometry.iloc[<span class="dv">0</span>].coords[<span class="dv">0</span>])</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="27">
<pre><code>32630</code></pre>
</div>
</div>
<p>Currently, we also have tools helping us to select a proper CRS. For example, the webpage <a href="https://crs-explorer.proj.org/" class="uri">https://crs-explorer.proj.org/</a> lists CRSs based on selected location and type. Important note: while these tools are helpful in many situations, you need to be aware of the properties of the recommended CRS before you apply it.</p>
<p>In cases where an appropriate CRS is not immediately clear, the choice of CRS should depend on the properties that are most important to preserve in the subsequent maps and analysis. All CRSs are either equal-area, equidistant, conformal (with shapes remaining unchanged), or some combination of compromises of those (<a href="01-spatial-data.html#sec-projected-coordinate-reference-systems" class="quarto-xref"><span>Section 1.4.2</span></a>). Custom CRSs with local parameters can be created for a region of interest and multiple CRSs can be used in projects when no single CRS suits all tasks. ‘Geodesic calculations’ can provide a fall-back if no CRSs are appropriate<a href="#fn5" class="footnote-ref" id="fnref5" role="doc-noteref"><sup>5</sup></a>. Regardless of the projected CRS used, the results may not be accurate for geometries covering hundreds of kilometers.</p>
<p>When deciding on a custom CRS, we recommend the following:</p>
<ul>
<li>A Lambert azimuthal equal-area (LAEA) projection for a custom local projection (set latitude and longitude of origin to the center of the study area), which is an equal-area projection at all locations but distorts shapes beyond thousands of kilometers</li>
<li>Azimuthal equidistant (AEQD) projections for a specifically accurate straight-line distance between a point and the center point of the local projection</li>
<li>Lambert conformal conic (LCC) projections for regions covering thousands of kilometers, with the cone set to keep distance and area properties reasonable between the secant lines</li>
<li>Stereographic (STERE) projections for polar regions, but taking care not to rely on area and distance calculations thousands of kilometers from the center</li>
</ul>
<p>One possible approach to automatically select a projected CRS specific to a local dataset is to create an azimuthal equidistant (AEQD) projection for the center-point of the study area. This involves creating a custom CRS (with no EPSG code) with units of meters based on the center point of a dataset. Note that this approach should be used with caution: no other datasets will be compatible with the custom CRS created and results may not be accurate when used on extensive datasets covering hundreds of kilometers.</p>
<p>The principles outlined in this section apply equally to vector and raster datasets. Some features of CRS transformation however are unique to each geographic data model. We will cover the particularities of vector data transformation in <a href="#sec-reprojecting-vector-geometries" class="quarto-xref"><span>Section 6.7</span></a> and those of raster transformation in <a href="#sec-reprojecting-raster-geometries" class="quarto-xref"><span>Section 6.8</span></a>. The last section, <a href="#sec-custom-map-projections" class="quarto-xref"><span>Section 6.9</span></a>, shows how to create custom map projections.</p>
</section>
<section id="sec-reprojecting-vector-geometries" class="level2" data-number="6.7">
<h2 data-number="6.7" class="anchored" data-anchor-id="sec-reprojecting-vector-geometries"><span class="header-section-number">6.7</span> Reprojecting vector geometries</h2>
<p><a href="01-spatial-data.html#sec-vector-data" class="quarto-xref"><span>Section 1.2</span></a> demonstrated how vector geometries are made-up of points, and how points form the basis of more complex objects such as lines and polygons. Reprojecting vectors thus consists of transforming the coordinates of these points, which form the vertices of lines and polygons.</p>
<p><a href="#sec-geometry-operations-on-projected-and-unprojected-data" class="quarto-xref"><span>Section 6.4</span></a> contains an example in which at a <code>GeoDataFrame</code> had to be transformed into an equivalent object, with a different CRS, to calculate the distance between two objects. Reprojection of vector layers is done using the <code>.to_crs</code> method.</p>
<div id="50d3e9c3" class="cell" data-execution_count="28">
<div class="sourceCode cell-code" id="cb42"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb42-1"><a href="#cb42-1" aria-hidden="true" tabindex="-1"></a>lnd_layer2 <span class="op">=</span> lnd_layer.to_crs(<span class="dv">27700</span>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Now that a transformed version of <code>lnd_layer</code> has been created, the distance between the two representations of London can be found using the <code>.distance</code> method.</p>
<div id="55a2956f" class="cell" data-execution_count="29">
<div class="sourceCode cell-code" id="cb43"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb43-1"><a href="#cb43-1" aria-hidden="true" tabindex="-1"></a>lnd_layer2.distance(lnd_layer_proj)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="29">
<pre><code>0 2017.949587
dtype: float64</code></pre>
</div>
</div>
<p>It may come as a surprise that <code>lnd_layer</code> and <code>lnd_layer2</code> are just over 2 <span class="math inline">\(km\)</span> apart! The difference in location between the two points is not due to imperfections in the transforming operation (which is in fact very accurate) but the low precision of the manually specified coordinates when creating <code>lnd_layer</code> and <code>lnd_layer_proj</code>.</p>
<p>Reprojecting to a different CRS is also demonstrated below using <code>cycle_hire_osm</code>, a point layer that represents ‘docking stations’ where you can hire bicycles in London. The contents of the CRS object associated with a given geometry column are changed when the object’s CRS is transformed. In the code chunk below, we create a new version of <code>cycle_hire_osm</code> with a projected CRS.</p>
<div id="aa15e463" class="cell" data-execution_count="30">
<div class="sourceCode cell-code" id="cb45"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb45-1"><a href="#cb45-1" aria-hidden="true" tabindex="-1"></a>cycle_hire_osm_projected <span class="op">=</span> cycle_hire_osm.to_crs(<span class="dv">27700</span>)</span>
<span id="cb45-2"><a href="#cb45-2" aria-hidden="true" tabindex="-1"></a>cycle_hire_osm_projected.crs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="30">
<pre><code><Projected CRS: EPSG:27700>
Name: OSGB36 / British National Grid
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore.
- bounds: (-9.01, 49.75, 2.01, 61.01)
Coordinate Operation:
- name: British National Grid
- method: Transverse Mercator
Datum: Ordnance Survey of Great Britain 1936
- Ellipsoid: Airy 1830
- Prime Meridian: Greenwich</code></pre>
</div>
</div>
<p>The resulting object has a new CRS according to the EPSG code <code>27700</code>. How to find out more details about this EPSG code, or any code? One option is to search for it online. Another option is to create a standalone CRS object within the Python environment (using <code>pyproj.CRS.from_string</code> or <code>pyproj.CRS.from_epsg</code>, see <a href="#sec-coordinate-reference-systems" class="quarto-xref"><span>Section 6.2</span></a>), and then query its properties, such as <code>.name</code> and <code>.to_wkt()</code>.</p>
<div id="42b04989" class="cell" data-execution_count="31">
<div class="sourceCode cell-code" id="cb47"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb47-1"><a href="#cb47-1" aria-hidden="true" tabindex="-1"></a>crs_lnd_new <span class="op">=</span> pyproj.CRS.from_epsg(<span class="dv">27700</span>)</span>
<span id="cb47-2"><a href="#cb47-2" aria-hidden="true" tabindex="-1"></a>crs_lnd_new.name, crs_lnd_new.to_wkt()</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="31">
<pre><code>('OSGB36 / British National Grid',
'PROJCRS["OSGB36 / British National Grid",BASEGEOGCRS["OSGB36",DATUM["Ordnance Survey of Great Britain 1936",ELLIPSOID["Airy 1830",6377563.396,299.3249646,LENGTHUNIT["metre",1]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4277]],CONVERSION["British National Grid",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",49,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-2,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996012717,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",400000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",-100000,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Engineering survey, topographic mapping."],AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45\'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],BBOX[49.75,-9.01,61.01,2.01]],ID["EPSG",27700]]')</code></pre>
</div>
</div>
<p>The result shows that the EPSG code <code>27700</code> represents the British National Grid, a result that could have been found by searching online for ‘EPSG 27700’.</p>
</section>
<section id="sec-reprojecting-raster-geometries" class="level2" data-number="6.8">
<h2 data-number="6.8" class="anchored" data-anchor-id="sec-reprojecting-raster-geometries"><span class="header-section-number">6.8</span> Reprojecting raster geometries</h2>
<p>The CRSs concepts described in the previous section apply equally to rasters. However, there are important differences in reprojection of vectors and rasters: transforming a vector object involves changing the coordinates of every vertex, but this does not apply to raster data. Rasters are composed of rectangular cells of the same size (expressed by map units, such as degrees or meters), so it is usually impracticable to transform coordinates of pixels separately. Raster reprojection involves creating a new raster object in the destination CRS, often with a different number of columns and rows than the original. The attributes must subsequently be re-estimated, allowing the new pixels to be ‘filled’ with appropriate values. In other words, raster reprojection can be thought of as two separate spatial operations: a vector reprojection of the raster extent to another CRS (<a href="#sec-reprojecting-vector-geometries" class="quarto-xref"><span>Section 6.7</span></a>), and computation of new pixel values through resampling (<a href="04-geometry-operations.html#sec-raster-resampling" class="quarto-xref"><span>Section 4.3.3</span></a>). Due to this additional complexity, in most cases when both raster and vector data are used, it is better to avoid reprojecting rasters and reproject vectors instead.</p>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>Reprojection of regular rasters is also known as warping. Additionally, there is a second similar operation called ‘transformation’. Instead of resampling all of the values, it leaves all values intact but recomputes new coordinates for every raster cell, changing the grid geometry. For example, it could convert the input raster (a regular grid) into a curvilinear grid. The <strong>rasterio</strong>, like common raster file formats (such as GeoTIFF), does not support curvilinear grids. The <strong>xarray</strong> package, for instance, can be used to work with curvilinear grids.</p>
</div>
</div>
<p>The raster reprojection process is done using two functions from the <code>rasterio.warp</code> sub-package:</p>
<ol type="1">
<li><code>rasterio.warp.calculate_default_transform</code>, used to calculate the new transformation matrix in the destination CRS, according to the source raster dimensions and bounds. Alternatively, the destination transformation matrix can be obtained from an existing raster; this is common practice when we need to align one raster with another, for instance to be able to combine them in raster algebra operations (<a href="03-spatial-operations.html#sec-raster-local-operations" class="quarto-xref"><span>Section 3.3.3</span></a>) (see below)</li>
<li><code>rasterio.warp.reproject</code>, introduced in <a href="04-geometry-operations.html#sec-raster-resampling" class="quarto-xref"><span>Section 4.3.3</span></a>, calculates cell values in the destination grid, using the user-selected resampling method (such as nearest neighbor, or bilinear)</li>
</ol>
<p>Let’s take a look at two examples of raster transformation: using categorical and continuous data. Land cover data are usually represented by categorical maps. The <code>nlcd.tif</code> file provides information for a small area in Utah, USA obtained from National Land Cover Database 2011 in the NAD83 / UTM zone 12N CRS. We already created a connection to the <code>nlcd.tif</code> file at the beginning of this chapter, named <code>src_nlcd</code>.</p>
<div id="c2f6b47d" class="cell" data-execution_count="32">
<div class="sourceCode cell-code" id="cb49"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb49-1"><a href="#cb49-1" aria-hidden="true" tabindex="-1"></a>src_nlcd</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="32">
<pre><code><open DatasetReader name='data/nlcd.tif' mode='r'></code></pre>
</div>
</div>
<p>Recall from previous chapters that the raster transformation matrix and dimensions are accessible from the file connection using <code>src_nlcd.transform</code>, <code>src_nlcd.width</code>, <code>src_nlcd.height</code>, and <code>src_nlcd.bounds</code>, respectively.</p>
<p>This information will be required to calculate the destination transformation matrix.</p>
<p>First, let’s define the destination CRS. In this case, we choose WGS84 (EPSG code <code>4326</code>).</p>
<div id="cee6a18b" class="cell" data-execution_count="33">
<div class="sourceCode cell-code" id="cb51"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb51-1"><a href="#cb51-1" aria-hidden="true" tabindex="-1"></a>dst_crs <span class="op">=</span> <span class="st">'EPSG:4326'</span></span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Now, we are ready to calculate the destination raster transformation matrix (<code>dst_transform</code>), and the destination dimensions (<code>dst_width</code>, <code>dst_height</code>), using <code>rasterio.warp.calculate_default_transform</code>, as follows:</p>
<div id="e042e39b" class="cell" data-execution_count="34">
<div class="sourceCode cell-code" id="cb52"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb52-1"><a href="#cb52-1" aria-hidden="true" tabindex="-1"></a>dst_transform, dst_width, dst_height <span class="op">=</span> rasterio.warp.calculate_default_transform(</span>
<span id="cb52-2"><a href="#cb52-2" aria-hidden="true" tabindex="-1"></a> src_nlcd.crs,</span>
<span id="cb52-3"><a href="#cb52-3" aria-hidden="true" tabindex="-1"></a> dst_crs,</span>
<span id="cb52-4"><a href="#cb52-4" aria-hidden="true" tabindex="-1"></a> src_nlcd.width,</span>
<span id="cb52-5"><a href="#cb52-5" aria-hidden="true" tabindex="-1"></a> src_nlcd.height,</span>
<span id="cb52-6"><a href="#cb52-6" aria-hidden="true" tabindex="-1"></a> <span class="op">*</span>src_nlcd.bounds</span>
<span id="cb52-7"><a href="#cb52-7" aria-hidden="true" tabindex="-1"></a>)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Here is the result.</p>
<div id="e3f76798" class="cell" data-execution_count="35">
<div class="sourceCode cell-code" id="cb53"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb53-1"><a href="#cb53-1" aria-hidden="true" tabindex="-1"></a>dst_transform</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="35">
<pre><code>Affine(0.00031506316853514724, 0.0, -113.24138811813536,
0.0, -0.00031506316853514724, 37.51912722777022)</code></pre>
</div>
</div>
<div id="6974a926" class="cell" data-execution_count="36">
<div class="sourceCode cell-code" id="cb55"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb55-1"><a href="#cb55-1" aria-hidden="true" tabindex="-1"></a>dst_width</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="36">
<pre><code>1244</code></pre>
</div>
</div>
<div id="6d1ae733" class="cell" data-execution_count="37">
<div class="sourceCode cell-code" id="cb57"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb57-1"><a href="#cb57-1" aria-hidden="true" tabindex="-1"></a>dst_height</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="37">
<pre><code>1246</code></pre>
</div>
</div>
<div class="callout callout-style-default callout-note callout-titled">
<div class="callout-header d-flex align-content-center">
<div class="callout-icon-container">
<i class="callout-icon"></i>
</div>
<div class="callout-title-container flex-fill">
Note
</div>
</div>
<div class="callout-body-container callout-body">
<p>The <code>*</code> syntax in Python is known as variable-length ‘<em>positional</em> arguments’. It is used to pass a <code>list</code> or <code>tuple</code> (or other iterables object) to positional arguments of a function.</p>
<p>For example, in the last code block, <code>*</code>, in <code>*src_nlcd.bounds</code>, is used to unpack <code>src_nlcd.bounds</code> (an iterable of length 4) to four separate arguments (<code>left</code>, <code>bottom</code>, <code>right</code>, and <code>top</code>), which <code>rasterio.warp.calculate_default_transform</code> requires in that order. In other words, the expression from the last example:</p>
<pre><code>rasterio.warp.calculate_default_transform(
src_nlcd.crs,
dst_crs,
src_nlcd.width,
src_nlcd.height,
*src_nlcd.bounds
)</code></pre>
<p>is a shortcut of:</p>
<pre><code>rasterio.warp.calculate_default_transform(
src_nlcd.crs,
dst_crs,
src_nlcd.width,
src_nlcd.height,
src_nlcd.bounds[0],
src_nlcd.bounds[1],
src_nlcd.bounds[2],
src_nlcd.bounds[3]
)</code></pre>
<p>‘<em>Keyword</em> arguments’ is a related technique; see note in <a href="04-geometry-operations.html#sec-raster-agg-disagg" class="quarto-xref"><span>Section 4.3.2</span></a>.</p>
</div>
</div>
<p>Recall from <a href="04-geometry-operations.html#sec-raster-resampling" class="quarto-xref"><span>Section 4.3.3</span></a> that resampling using <code>rasterio.warp.reproject</code> can take place directly into a ‘destination’ raster file connection. Therefore, our next step is to create the metadata file used for writing the reprojected raster to file. For convenience, we are taking the metadata of the source raster (<code>src_nlcd.meta</code>), making a copy (<code>dst_kwargs</code>), and then updating those specific properties that need to be changed. Note that the reprojection process typically creates ‘No Data’ pixels, even when there were none in the input raster, since the raster orientation changes and the edges need to be ‘filled’ to get back a rectangular extent. For example, a reprojected raster may appear as a ‘tilted’ rectangle, inside a larger straight rectangular extent, whereas the margins around the tilted rectangle are inevitably filled with ‘No Data’ (e.g., the white stripes surrounding the edges in <a href="#fig-raster-reproject-nlcd" class="quarto-xref">Figure <span>6.2</span></a> (b) are ‘No Data’ pixels created as a result of reprojection). We need to specify a ‘No Data’ value of our choice, if there is no existing definition, or keep the existing source raster ‘No Data’ setting, such as <code>255</code> in this case.</p>
<div id="efd3d2d9" class="cell" data-execution_count="38">
<div class="sourceCode cell-code" id="cb61"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb61-1"><a href="#cb61-1" aria-hidden="true" tabindex="-1"></a>dst_kwargs <span class="op">=</span> src_nlcd.meta.copy()</span>
<span id="cb61-2"><a href="#cb61-2" aria-hidden="true" tabindex="-1"></a>dst_kwargs.update({</span>
<span id="cb61-3"><a href="#cb61-3" aria-hidden="true" tabindex="-1"></a> <span class="st">'crs'</span>: dst_crs,</span>
<span id="cb61-4"><a href="#cb61-4" aria-hidden="true" tabindex="-1"></a> <span class="st">'transform'</span>: dst_transform,</span>
<span id="cb61-5"><a href="#cb61-5" aria-hidden="true" tabindex="-1"></a> <span class="st">'width'</span>: dst_width,</span>
<span id="cb61-6"><a href="#cb61-6" aria-hidden="true" tabindex="-1"></a> <span class="st">'height'</span>: dst_height</span>
<span id="cb61-7"><a href="#cb61-7" aria-hidden="true" tabindex="-1"></a>})</span>
<span id="cb61-8"><a href="#cb61-8" aria-hidden="true" tabindex="-1"></a>dst_kwargs</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
<div class="cell-output cell-output-display" data-execution_count="38">
<pre><code>{'driver': 'GTiff',
'dtype': 'uint8',
'nodata': 255.0,
'width': 1244,
'height': 1246,
'count': 1,
'crs': 'EPSG:4326',
'transform': Affine(0.00031506316853514724, 0.0, -113.24138811813536,
0.0, -0.00031506316853514724, 37.51912722777022)}</code></pre>
</div>
</div>
<p>Now, we are ready to create the reprojected raster. Here, reprojection takes place between two file connections, meaning that the raster value arrays are not being read into memory at once. (It is also possible to reproject into an in-memory <code>ndarray</code> object.)</p>
<p>To write the reprojected raster, we first create a destination file connection <code>dst_nlcd</code>, pointing at the output file path of our choice (<code>'output/nlcd_4326.tif'</code>), using the updated metadata object created earlier (<code>dst_kwargs</code>):</p>
<div id="fdf4787a" class="cell" data-execution_count="39">
<div class="sourceCode cell-code" id="cb63"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb63-1"><a href="#cb63-1" aria-hidden="true" tabindex="-1"></a>dst_nlcd <span class="op">=</span> rasterio.<span class="bu">open</span>(<span class="st">'output/nlcd_4326.tif'</span>, <span class="st">'w'</span>, <span class="op">**</span>dst_kwargs)</span></code><button title="Copy to Clipboard" class="code-copy-button"><i class="bi"></i></button></pre></div>
</div>
<p>Then, we use the <code>rasterio.warp.reproject</code> function to calculate and write the reprojection result into the <code>dst_nlcd</code> file connection.</p>
<div id="a5a3c9ce" class="cell" data-execution_count="40">
<div class="sourceCode cell-code" id="cb64"><pre class="sourceCode python code-with-copy"><code class="sourceCode python"><span id="cb64-1"><a href="#cb64-1" aria-hidden="true" tabindex="-1"></a>rasterio.warp.reproject(</span>
<span id="cb64-2"><a href="#cb64-2" aria-hidden="true" tabindex="-1"></a> source<span class="op">=</span>rasterio.band(src_nlcd, <span class="dv">1</span>),</span>