2 + 2 #Any part of a line in`R`that starts with a # will be ignored with the code is run.[1] 4#We use lines (or parts of lines) that start with # to make notes in the code.Introducing Geospatial Analysis in R
Introduction
In this tutorial, have a gentle introduction to coding in R and help you get familiar with how the system works.
It is important to remember that writing code in R is very tricky - every capital letter and punctuation mark needs to be in the right place, or the code won’t work.
Luckily, the easiest way to learn to code in R is to look at examples and then to switch out parts of those examples to see what happens. You can use that strategy to learn how the R language works for the tasks that you are trying to accomplish. That’s the strategy that I’m going to be using to teach you in this class.
The interesting thing about R is that pretty much everything you need can be written inRscripts. TheRscript is a series of commands forRthat it follows to conduct analysis, build maps, and so on.
Scripts are really useful because they are also a record of every step of your analysis. That means that if you wind up making a mistake that you don’t notice until much later, all you need to do is find the source of the problem in your script, fix it, and rerun your code from there. That’s a lot better than having to redo all your work from scratch!
I’m using a special format ofRscripts called a Quarto Document, which allows you to combine the script with HTML text. When I generate the HTML, it runs all the code blocks in the document in order and puts the results below. That means you can check what you’re getting when you run the code against what I got when I ran it.
Quarto Documents are handy for making class materials, but they can also be used to make general websites or even presentations. We’ll see how to use those in this class.
The Objectives
- Become familiar with the basics of the
Rlanguage - Become familiar with editing and running
Rcode in RStudio - Learn how to install and use basic geospatial package in
R - Learn how to plot vector objects
- Learn how to access vector data from websites and add them directly to
R
The Data
Most of the data in this tutorial will be downloaded directly from the web and read into R, but we do have two datasets that we will be using.
- jerash.gpkg contains data produced by Stott, et al. (2018). Mapping an ancient city with a century of remotely sensed data. PNAS, 115(34), E5450-E5458, who digitized excavation maps for the Graeco-Roman city of Jerash, in modern-day Jordan
- carleton_Angkor_temples.csv contains the locations and estimated construction dates of temples that have been idetified in the Greater Angkor complex, as reported in Carleton, et al. (2023). Bayesian regression versus machine learning for rapid age estimation of archaeological features identified with lidar at Angkor. Scientific Reports, 13, 17913.
These data are available in this zip file.
Introduction to RStudio
Ris a programming language, while RStudio is an environment to help you edit and runR, so when you say you have performed an analysis, you say you performed it inR, rather than in RStudio- Have we all successfully installed
Rand RStudio? - Tour of main features
- Introduction to Quarto documents
First Lines of Code
Let’s start with your first line of code. You can copy and paste the following into your script editor in RStudio to run it:
Once you have the code pasted in your script editor, you highlight the code with the cursor and then click “Run”, which is in the upper right-hand corner of this script editor window, to run it. Try running the code above now.
Congratulations! That’s your first code. This is pretty basic. Among many other things,Rworks as a simple calculator. It also obeys standard order of operations rules. Run the lines below for examples:
2 + (3 * 2)[1] 8(2 + 3) * 2[1] 104 / 2[1] 2(4 + 2) / 2[1] 34 + (2 / 2)[1] 5Okay, so you get the idea!
Object-Oriented Programming in R
R is an object-oriented programming language. What that means is thatRworks by creating snippets of information in the computer memory that we call “objects” and then doing things to those objects.
I often think about the data that I am working on almost like a physical structure that I’m usingRto build and manipulate. Another helpful metaphor might be to think of objects as little digital creatures that you useRto communicate with.
However it feels easiest to think about it, let’s introduce your first object by running the following line:
x <- 5The <- sign in that line is the assignment function. In R, functions are pieces of code that do things to objects. In those terms, all the mathematical symbols we used above would count as functions.
The assignment function takes whatever is on the right-hand side of the function and creates an object according to the instructions on the left-hand side that contains the thing on the right-hand side.
In this case, we created an object we called x and gave it the content of the number 5. We can then type in the name of the object to askRto print back its contents:
x[1] 5Importantly,Ris case sensitive. That means it cares whether or not letters are capitalized. Try running this:
XError in eval(expr, envir, enclos): object 'X' not foundR should have given you this message: Error: object 'X' not found. When R says an object is not found, it’s telling you that something that you think is read into the environment (that is, the collection of objects thatRhas in memory at the moment) isn’t actually there.
If you look in the Environment tab in RStudio, which is probably to the right of your script editing window, you can see a list of all the objects that you have created in the environment. You should see that there is an object called “x” there. Notice that “x” isn’t capitalized. That’s whyRcan’t find it. R considers x and X to be completely different.
You can pack pretty much anything thatRcan hold in memory into an object. We’ll be putting whole layers of spatial data into objects soon! For now, let me just show you a bit of how this works:
x <- c(1, 5, 10)In this line, we used the c, or concatenate function, which takes a list of items separated by commas and puts them into a single object called a vector (not to be confused with the idea of a vector layer in GIS), which is just a line of data:
x[1] 1 5 10The parentheses and commas are really, really important in R. Commas separate different inputs into a function. We call these separate inputs arguments. The parentheses tellRwhere the arguments for each function start and end.
See what happens if we mess this up:
x <- c(1 5, 10)
x <- c(1, 5) 10Error: <text>:1:10: unexpected numeric constant
1: x <- c(1 5
^Both of those situations causeRto say Error: unexpected . . .. When you see an error that starts that way, it’s R’s way of telling you that there’s something wrong with the way you typed your code and that things are out of order.
RStudio actually helps you keep track of coding problems. In the first of the two lines above, for example, you should see a red underline under 5. That’s RStudio telling you that something fishy is going on with your code there (in this case, you’re missing the comma after 1).
Okay, let’s get our x back in working order:
x <- c(1, 2, 3, 4, 5, 10)Basic Vector Operations
Once you’ve created a vector, you can use functions to perform operations on all the elements of the vector at once:
x <- c(1, 5, 10)
x + 2[1] 3 7 12x * 2[1] 2 10 20x ^ 2[1] 1 25 100x / 2[1] 0.5 2.5 5.0You can also use statistical functions on the whole vector. Here are some common handy examples:
min(x)[1] 1mean(x)[1] 5.333333median(x)[1] 5max(x)[1] 10sd(x) # Standard Deviation[1] 4.50925summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 3.000 5.000 5.333 7.500 10.000 Now that we have some sense of how R works, we’ll turn our attention to how to use GIS data in R.
The data
- Global country boundaries and basic socioeconomic data from Natural Earth, sourced via the
rnaturalearthpackage. - Archaeological GIS data from ancient Roman city of Jerash, in modern-day Jordan, from Stott, D., et al (2018). Mapping an ancient city with a century of remotely sensed data. Proceedings of the National Academies of Science, 115(24): E5450-E5458.
- Archaeological GIS data on temple locations in Greater Angkor, in modern-day Cambodia, from Carleton, W. C., et al. (2023). Bayesian regression versus machine learning for rapid age estimation of archaeological features identified at Angkor. Scientific Reports, 13, 17913.
Packages
Packages are additional chunks of code that extend R’s functionality
Today, we will be using:
tidyverse: Set of packages that facilitate easy data cleaning, transformation, and visualizationsf: Main package for working with vector data inRterra: Main package for working with raster data inRtidyterra: Package that allows you to easily plotterraobjects alongsidesfobjects on a maprnaturalearth: Package that allows us to import cartographic data from Natural Earth directly intoRggspatial: Package for plotting geospatial data usingtidyverse’sggplot2packagemaptiles: Package that allows you to extract map tiles from online services like OpenStreetMap and use them in your own maps
To use a package, it must first be installed on your computer: * require() will attempt to load a package and let you know if it needs to be installed * Using install_package() with the package name in quotations will install the named package
require(devtools) # Used to install a support package for rnatural earthLoading required package: devtoolsWarning: package 'devtools' was built under R version 4.3.3Loading required package: usethisWarning: package 'usethis' was built under R version 4.3.3devtools::install_github("ropensci/rnaturalearthhires") # The support packageSkipping install of 'rnaturalearthhires' from a github remote, the SHA1 (dd1e210c) has not changed since last install.
Use `force = TRUE` to force installationrequire(rnaturalearth)Loading required package: rnaturalearthWarning: package 'rnaturalearth' was built under R version 4.3.3require(tidyverse)Loading required package: tidyverseWarning: package 'ggplot2' was built under R version 4.3.3── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.0
✔ purrr 1.0.2 ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsrequire(sf)Loading required package: sf
Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUErequire(terra)Loading required package: terraWarning: package 'terra' was built under R version 4.3.3terra 1.7.78
Attaching package: 'terra'
The following object is masked from 'package:tidyr':
extractrequire(tidyterra)Loading required package: tidyterraWarning: package 'tidyterra' was built under R version 4.3.3
Attaching package: 'tidyterra'
The following object is masked from 'package:stats':
filterrequire(ggspatial)Loading required package: ggspatialrequire(maptiles)Loading required package: maptilesWarning: package 'maptiles' was built under R version 4.3.3Plotting vector data with R
- Starting with country boundaries data from Natural Earth
ne_countries()function automatically downloads Natural Earth country boundaries, along with some data on the countries, and imports them intoRas ansf(vector) object
countries <- ne_countries(scale = "large", returnclass = "sf")
countriesSimple feature collection with 258 features and 168 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
Geodetic CRS: WGS 84
First 10 features:
featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level
1 Admin-0 country 0 2 Indonesia IDN 0 2
2 Admin-0 country 0 3 Malaysia MYS 0 2
3 Admin-0 country 0 2 Chile CHL 0 2
4 Admin-0 country 0 3 Bolivia BOL 0 2
5 Admin-0 country 0 2 Peru PER 0 2
6 Admin-0 country 0 2 Argentina ARG 0 2
7 Admin-0 country 3 3 United Kingdom GB1 1 2
8 Admin-0 country 1 5 Cyprus CYP 0 2
9 Admin-0 country 0 2 India IND 0 2
10 Admin-0 country 0 2 China CH1 1 2
type tlc admin adm0_a3 geou_dif
1 Sovereign country 1 Indonesia IDN 0
2 Sovereign country 1 Malaysia MYS 0
3 Sovereign country 1 Chile CHL 0
4 Sovereign country 1 Bolivia BOL 0
5 Sovereign country 1 Peru PER 0
6 Sovereign country 1 Argentina ARG 0
7 Dependency 1 Dhekelia Sovereign Base Area ESB 0
8 Sovereign country 1 Cyprus CYP 0
9 Sovereign country 1 India IND 0
10 Country 1 China CHN 0
geounit gu_a3 su_dif subunit su_a3
1 Indonesia IDN 0 Indonesia IDN
2 Malaysia MYS 0 Malaysia MYS
3 Chile CHL 0 Chile CHL
4 Bolivia BOL 0 Bolivia BOL
5 Peru PER 0 Peru PER
6 Argentina ARG 0 Argentina ARG
7 Dhekelia Sovereign Base Area ESB 0 Dhekelia Sovereign Base Area ESB
8 Cyprus CYP 0 Cyprus CYP
9 India IND 0 India IND
10 China CHN 0 China CHN
brk_diff name name_long brk_a3 brk_name brk_group abbrev postal
1 0 Indonesia Indonesia IDN Indonesia <NA> Indo. INDO
2 0 Malaysia Malaysia MYS Malaysia <NA> Malay. MY
3 0 Chile Chile CHL Chile <NA> Chile CL
4 0 Bolivia Bolivia BOL Bolivia <NA> Bolivia BO
5 0 Peru Peru PER Peru <NA> Peru PE
6 0 Argentina Argentina ARG Argentina <NA> Arg. AR
7 0 Dhekelia Dhekelia ESB Dhekelia <NA> Dhek. DH
8 0 Cyprus Cyprus CYP Cyprus <NA> Cyp. CY
9 0 India India IND India <NA> India IND
10 0 China China CHN China <NA> China CN
formal_en formal_fr name_ciawf note_adm0 note_brk
1 Republic of Indonesia <NA> Indonesia <NA> <NA>
2 Malaysia <NA> Malaysia <NA> <NA>
3 Republic of Chile <NA> Chile <NA> <NA>
4 Plurinational State of Bolivia <NA> Bolivia <NA> <NA>
5 Republic of Peru <NA> Peru <NA> <NA>
6 Argentine Republic <NA> Argentina <NA> <NA>
7 <NA> <NA> <NA> U.K. U.K. Base
8 Republic of Cyprus <NA> Cyprus <NA> <NA>
9 Republic of India <NA> India <NA> <NA>
10 People's Republic of China <NA> China <NA> <NA>
name_sort name_alt mapcolor7 mapcolor8 mapcolor9
1 Indonesia <NA> 6 6 6
2 Malaysia <NA> 2 4 3
3 Chile <NA> 5 1 5
4 Bolivia <NA> 1 5 2
5 Peru <NA> 4 4 4
6 Argentina <NA> 3 1 3
7 Dhekelia Sovereign Base Area <NA> 6 6 6
8 Cyprus <NA> 1 2 3
9 India <NA> 1 3 2
10 China <NA> 4 4 4
mapcolor13 pop_est pop_rank pop_year gdp_md gdp_year
1 11 270625568 17 2019 1119190 2019
2 6 31949777 15 2019 364681 2019
3 9 18952038 14 2019 282318 2019
4 3 11513100 14 2019 40895 2019
5 11 32510453 15 2019 226848 2019
6 13 44938712 15 2019 445445 2019
7 3 7850 5 2013 314 2013
8 7 1198575 12 2019 24948 2019
9 2 1366417754 18 2019 2868929 2019
10 3 1397715000 18 2019 14342903 2019
economy income_grp fips_10 iso_a2 iso_a2_eh
1 4. Emerging region: MIKT 4. Lower middle income ID ID ID
2 6. Developing region 3. Upper middle income MY MY MY
3 5. Emerging region: G20 3. Upper middle income CI CL CL
4 5. Emerging region: G20 4. Lower middle income BL BO BO
5 5. Emerging region: G20 3. Upper middle income PE PE PE
6 5. Emerging region: G20 3. Upper middle income AR AR AR
7 2. Developed region: nonG7 2. High income: nonOECD -99 -99 -99
8 6. Developing region 2. High income: nonOECD CY CY CY
9 3. Emerging region: BRIC 4. Lower middle income IN IN IN
10 3. Emerging region: BRIC 3. Upper middle income CH CN CN
iso_a3 iso_a3_eh iso_n3 iso_n3_eh un_a3 wb_a2 wb_a3 woe_id woe_id_eh
1 IDN IDN 360 360 360 ID IDN 23424846 23424846
2 MYS MYS 458 458 458 MY MYS 23424901 23424901
3 CHL CHL 152 152 152 CL CHL 23424782 23424782
4 BOL BOL 068 068 068 BO BOL 23424762 23424762
5 PER PER 604 604 604 PE PER 23424919 23424919
6 ARG ARG 032 032 032 AR ARG 23424747 23424747
7 -99 -99 -99 -99 -099 -99 -99 -99 -99
8 CYP CYP 196 196 196 CY CYP -90 23424994
9 IND IND 356 356 356 IN IND 23424848 23424848
10 CHN CHN 156 156 156 CN CHN 23424781 23424781
woe_note adm0_iso adm0_diff adm0_tlc adm0_a3_us
1 Exact WOE match as country IDN <NA> IDN IDN
2 Exact WOE match as country MYS <NA> MYS MYS
3 Exact WOE match as country CHL <NA> CHL CHL
4 Exact WOE match as country BOL <NA> BOL BOL
5 Exact WOE match as country PER <NA> PER PER
6 Exact WOE match as country ARG <NA> ARG ARG
7 No WOE equivalent. -99 1 ESB ESB
8 WOE lists as subunit of united Cyprus CYP <NA> CYP CYP
9 Exact WOE match as country IND <NA> IND IND
10 Exact WOE match as country CHN <NA> CHN CHN
adm0_a3_fr adm0_a3_ru adm0_a3_es adm0_a3_cn adm0_a3_tw adm0_a3_in adm0_a3_np
1 IDN IDN IDN IDN IDN IDN IDN
2 MYS MYS MYS MYS MYS MYS MYS
3 CHL CHL CHL CHL CHL CHL CHL
4 BOL BOL BOL BOL BOL BOL BOL
5 PER PER PER PER PER PER PER
6 ARG ARG ARG ARG ARG ARG ARG
7 ESB ESB ESB ESB ESB ESB ESB
8 CYP CYP CYP CYP CYP CYP CYP
9 IND IND IND IND IND IND IND
10 CHN CHN CHN CHN TWN CHN CHN
adm0_a3_pk adm0_a3_de adm0_a3_gb adm0_a3_br adm0_a3_il adm0_a3_ps adm0_a3_sa
1 IDN IDN IDN IDN IDN IDN IDN
2 MYS MYS MYS MYS MYS MYS MYS
3 CHL CHL CHL CHL CHL CHL CHL
4 BOL BOL BOL BOL BOL BOL BOL
5 PER PER PER PER PER PER PER
6 ARG ARG ARG ARG ARG ARG ARG
7 ESB ESB ESB ESB ESB ESB ESB
8 CYP CYP CYP CYP CYP CYP CYP
9 IND IND IND IND IND IND IND
10 CHN CHN CHN CHN CHN CHN CHN
adm0_a3_eg adm0_a3_ma adm0_a3_pt adm0_a3_ar adm0_a3_jp adm0_a3_ko adm0_a3_vn
1 IDN IDN IDN IDN IDN IDN IDN
2 MYS MYS MYS MYS MYS MYS MYS
3 CHL CHL CHL CHL CHL CHL CHL
4 BOL BOL BOL BOL BOL BOL BOL
5 PER PER PER PER PER PER PER
6 ARG ARG ARG ARG ARG ARG ARG
7 ESB ESB ESB ESB ESB ESB ESB
8 CYP CYP CYP CYP CYP CYP CYP
9 IND IND IND IND IND IND IND
10 CHN CHN CHN CHN CHN CHN CHN
adm0_a3_tr adm0_a3_id adm0_a3_pl adm0_a3_gr adm0_a3_it adm0_a3_nl adm0_a3_se
1 IDN IDN IDN IDN IDN IDN IDN
2 MYS MYS MYS MYS MYS MYS MYS
3 CHL CHL CHL CHL CHL CHL CHL
4 BOL BOL BOL BOL BOL BOL BOL
5 PER PER PER PER PER PER PER
6 ARG ARG ARG ARG ARG ARG ARG
7 ESB ESB ESB ESB ESB ESB ESB
8 CYP CYP CYP CYP CYP CYP CYP
9 IND IND IND IND IND IND IND
10 CHN CHN CHN CHN CHN CHN CHN
adm0_a3_bd adm0_a3_ua adm0_a3_un adm0_a3_wb continent region_un
1 IDN IDN -99 -99 Asia Asia
2 MYS MYS -99 -99 Asia Asia
3 CHL CHL -99 -99 South America Americas
4 BOL BOL -99 -99 South America Americas
5 PER PER -99 -99 South America Americas
6 ARG ARG -99 -99 South America Americas
7 ESB ESB -99 -99 Asia Asia
8 CYP CYP -99 -99 Asia Asia
9 IND IND -99 -99 Asia Asia
10 CHN CHN -99 -99 Asia Asia
subregion region_wb name_len long_len abbrev_len
1 South-Eastern Asia East Asia & Pacific 9 9 5
2 South-Eastern Asia East Asia & Pacific 8 8 6
3 South America Latin America & Caribbean 5 5 5
4 South America Latin America & Caribbean 7 7 7
5 South America Latin America & Caribbean 4 4 4
6 South America Latin America & Caribbean 9 9 4
7 Western Asia Europe & Central Asia 8 8 5
8 Western Asia Europe & Central Asia 6 6 4
9 Southern Asia South Asia 5 5 5
10 Eastern Asia East Asia & Pacific 5 5 5
tiny homepart min_zoom min_label max_label label_x label_y ne_id
1 -99 1 0 1.7 6.7 101.89295 -0.954404 1159320845
2 -99 1 0 3.0 8.0 113.83708 2.528667 1159321083
3 -99 1 0 1.7 6.7 -72.31887 -38.151771 1159320493
4 -99 1 0 3.0 7.5 -64.59343 -16.666015 1159320439
5 -99 1 0 2.0 7.0 -72.90016 -12.976679 1159321163
6 -99 1 0 2.0 7.0 -64.17333 -33.501159 1159320331
7 3 -99 0 6.5 11.0 33.79058 35.011042 1159320709
8 -99 1 0 4.5 9.5 33.08418 34.913329 1159320533
9 -99 1 0 1.7 6.7 79.35810 22.686852 1159320847
10 -99 1 0 1.7 5.7 106.33729 32.498178 1159320471
wikidataid name_ar name_bn name_de
1 Q252 إندونيسيا ইন্দোনেশিয়া Indonesien
2 Q833 ماليزيا মালয়েশিয়া Malaysia
3 Q298 تشيلي চিলি Chile
4 Q750 بوليفيا বলিভিয়া Bolivien
5 Q419 بيرو পেরু Peru
6 Q414 الأرجنتين আর্জেন্টিনা Argentinien
7 Q9206745 ديكيليا كانتونمنت দেখেলিয়া ক্যান্টনমেন্ট Dekelia
8 Q229 قبرص সাইপ্রাস Republik Zypern
9 Q668 الهند ভারত Indien
10 Q148 الصين গণচীন Volksrepublik China
name_en name_es name_fa
1 Indonesia Indonesia اندونزی
2 Malaysia Malasia مالزی
3 Chile Chile شیلی
4 Bolivia Bolivia بولیوی
5 Peru Perú پرو
6 Argentina Argentina آرژانتین
7 Dhekelia Cantonment Dekelia دکلیا
8 Cyprus Chipre قبرس
9 India India هند
10 People's Republic of China China جمهوری خلق چین
name_fr name_el
1 Indonésie Ινδονησία
2 Malaisie Μαλαισία
3 Chili Χιλή
4 Bolivie Βολιβία
5 Pérou Περού
6 Argentine Αργεντινή
7 Dhekelia Ντεκέλια Κάντονμεντ
8 Chypre Κύπρος
9 Inde Ινδία
10 République populaire de Chine Λαϊκή Δημοκρατία της Κίνας
name_he name_hi name_hu
1 אינדונזיה इंडोनेशिया Indonézia
2 מלזיה मलेशिया Malajzia
3 צ'ילה चिली Chile
4 בוליביה बोलिविया Bolívia
5 פרו पेरू Peru
6 ארגנטינה अर्जेण्टीना Argentína
7 דקליה ढेकेलिया छावनी Dekélia
8 קפריסין साइप्रस Ciprus
9 הודו भारत India
10 הרפובליקה העממית של סין चीनी जनवादी गणराज्य Kína
name_id name_it name_ja name_ko
1 Indonesia Indonesia インドネシア 인도네시아
2 Malaysia Malaysia マレーシア 말레이시아
3 Chili Cile チリ 칠레
4 Bolivia Bolivia ボリビア 볼리비아
5 Peru Perù ペルー 페루
6 Argentina Argentina アルゼンチン 아르헨티나
7 Dhekelia Cantonment Base di Dhekelia デケリア 데켈리아 지역
8 Siprus Cipro キプロス 키프로스
9 India India インド 인도
10 Republik Rakyat Tiongkok Cina 中華人民共和国 중화인민공화국
name_nl name_pl name_pt
1 Indonesië Indonezja Indonésia
2 Maleisië Malezja Malásia
3 Chili Chile Chile
4 Bolivia Boliwia Bolívia
5 Peru Peru Peru
6 Argentinië Argentyna Argentina
7 Dhekelia Cantonment Dhekelia Deceleia
8 Cyprus Cypr Chipre
9 India Indie Índia
10 Volksrepubliek China Chińska Republika Ludowa China
name_ru name_sv name_tr
1 Индонезия Indonesien Endonezya
2 Малайзия Malaysia Malezya
3 Чили Chile Şili
4 Боливия Bolivia Bolivya
5 Перу Peru Peru
6 Аргентина Argentina Arjantin
7 Декелия Dhekelia Dhekelia Kantonu
8 Кипр Cypern Kıbrıs Cumhuriyeti
9 Индия Indien Hindistan
10 Китайская Народная Республика Kina Çin Halk Cumhuriyeti
name_uk name_ur name_vi
1 Індонезія انڈونیشیا Indonesia
2 Малайзія ملائیشیا Malaysia
3 Чилі چلی Chile
4 Болівія بولیویا Bolivia
5 Перу پیرو Peru
6 Аргентина ارجنٹائن Argentina
7 Муніципалітет Декелія دحیکیلیا کانتونمینٹ Căn cứ quân sự Dhekelia
8 Кіпр قبرص Cộng hòa Síp
9 Індія بھارت Ấn Độ
10 Китайська Народна Республіка عوامی جمہوریہ چین Trung Quốc
name_zh name_zht fclass_iso tlc_diff fclass_tlc
1 印度尼西亚 印度尼西亞 Admin-0 country <NA> Admin-0 country
2 马来西亚 馬來西亞 Admin-0 country <NA> Admin-0 country
3 智利 智利 Admin-0 country <NA> Admin-0 country
4 玻利维亚 玻利維亞 Admin-0 country <NA> Admin-0 country
5 秘鲁 秘魯 Admin-0 country <NA> Admin-0 country
6 阿根廷 阿根廷 Admin-0 country <NA> Admin-0 country
7 泽凯利亚军营 德凱利亞軍營 Unrecognized 1 Admin-0 dependency
8 塞浦路斯 賽普勒斯 Admin-0 country <NA> Admin-0 country
9 印度 印度 Admin-0 country <NA> Admin-0 country
10 中华人民共和国 中華人民共和國 Admin-0 country <NA> Admin-0 country
fclass_us fclass_fr fclass_ru fclass_es fclass_cn
1 <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA> <NA>
7 Admin-0 dependency Admin-0 dependency <NA> Admin-0 dependency <NA>
8 <NA> <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA> <NA>
10 <NA> <NA> <NA> <NA> <NA>
fclass_tw fclass_in fclass_np fclass_pk fclass_de
1 <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA> <NA>
7 <NA> <NA> <NA> <NA> Admin-0 dependency
8 <NA> <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA> <NA>
10 Unrecognized <NA> <NA> <NA> <NA>
fclass_gb fclass_br fclass_il fclass_ps fclass_sa fclass_eg
1 <NA> <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA> <NA> <NA>
7 Admin-0 dependency <NA> <NA> <NA> <NA> <NA>
8 <NA> <NA> <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA> <NA> <NA>
10 <NA> <NA> <NA> <NA> <NA> <NA>
fclass_ma fclass_pt fclass_ar fclass_jp fclass_ko
1 <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA> <NA>
7 <NA> Admin-0 dependency <NA> <NA> Admin-0 dependency
8 <NA> <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA> <NA>
10 <NA> <NA> <NA> <NA> <NA>
fclass_vn fclass_tr fclass_id fclass_pl fclass_gr
1 <NA> <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA> <NA>
7 <NA> Admin-0 dependency <NA> Admin-0 dependency Admin-0 dependency
8 <NA> <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA> <NA>
10 <NA> <NA> <NA> <NA> <NA>
fclass_it fclass_nl fclass_se fclass_bd
1 <NA> <NA> <NA> <NA>
2 <NA> <NA> <NA> <NA>
3 <NA> <NA> <NA> <NA>
4 <NA> <NA> <NA> <NA>
5 <NA> <NA> <NA> <NA>
6 <NA> <NA> <NA> <NA>
7 Admin-0 dependency Admin-0 dependency Admin-0 dependency <NA>
8 <NA> <NA> <NA> <NA>
9 <NA> <NA> <NA> <NA>
10 <NA> <NA> <NA> <NA>
fclass_ua geometry
1 <NA> MULTIPOLYGON (((117.7036 4....
2 <NA> MULTIPOLYGON (((117.7036 4....
3 <NA> MULTIPOLYGON (((-69.51009 -...
4 <NA> MULTIPOLYGON (((-69.51009 -...
5 <NA> MULTIPOLYGON (((-69.51009 -...
6 <NA> MULTIPOLYGON (((-67.1939 -2...
7 Admin-0 dependency MULTIPOLYGON (((33.78094 34...
8 <NA> MULTIPOLYGON (((33.78183 34...
9 <NA> MULTIPOLYGON (((77.80035 35...
10 <NA> MULTIPOLYGON (((78.91769 33...How do we interpret all the information that we’re seeing here? What’s important?
names(countries) [1] "featurecla" "scalerank" "labelrank" "sovereignt" "sov_a3"
[6] "adm0_dif" "level" "type" "tlc" "admin"
[11] "adm0_a3" "geou_dif" "geounit" "gu_a3" "su_dif"
[16] "subunit" "su_a3" "brk_diff" "name" "name_long"
[21] "brk_a3" "brk_name" "brk_group" "abbrev" "postal"
[26] "formal_en" "formal_fr" "name_ciawf" "note_adm0" "note_brk"
[31] "name_sort" "name_alt" "mapcolor7" "mapcolor8" "mapcolor9"
[36] "mapcolor13" "pop_est" "pop_rank" "pop_year" "gdp_md"
[41] "gdp_year" "economy" "income_grp" "fips_10" "iso_a2"
[46] "iso_a2_eh" "iso_a3" "iso_a3_eh" "iso_n3" "iso_n3_eh"
[51] "un_a3" "wb_a2" "wb_a3" "woe_id" "woe_id_eh"
[56] "woe_note" "adm0_iso" "adm0_diff" "adm0_tlc" "adm0_a3_us"
[61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
[66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
[71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
[76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
[81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
[86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
[91] "adm0_a3_un" "adm0_a3_wb" "continent" "region_un" "subregion"
[96] "region_wb" "name_len" "long_len" "abbrev_len" "tiny"
[101] "homepart" "min_zoom" "min_label" "max_label" "label_x"
[106] "label_y" "ne_id" "wikidataid" "name_ar" "name_bn"
[111] "name_de" "name_en" "name_es" "name_fa" "name_fr"
[116] "name_el" "name_he" "name_hi" "name_hu" "name_id"
[121] "name_it" "name_ja" "name_ko" "name_nl" "name_pl"
[126] "name_pt" "name_ru" "name_sv" "name_tr" "name_uk"
[131] "name_ur" "name_vi" "name_zh" "name_zht" "fclass_iso"
[136] "tlc_diff" "fclass_tlc" "fclass_us" "fclass_fr" "fclass_ru"
[141] "fclass_es" "fclass_cn" "fclass_tw" "fclass_in" "fclass_np"
[146] "fclass_pk" "fclass_de" "fclass_gb" "fclass_br" "fclass_il"
[151] "fclass_ps" "fclass_sa" "fclass_eg" "fclass_ma" "fclass_pt"
[156] "fclass_ar" "fclass_jp" "fclass_ko" "fclass_vn" "fclass_tr"
[161] "fclass_id" "fclass_pl" "fclass_gr" "fclass_it" "fclass_nl"
[166] "fclass_se" "fclass_bd" "fclass_ua" "geometry" See any variables we might want to map?
map_1 <- ggplot(data = countries) +
geom_sf(aes(fill=income_grp))
map_1
- Kinda’ ugly; try changing color ramp and theme:
map_2 <- ggplot(data = countries) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer(type = "qual")
map_2
- Can change palette to get different types of colors:
map_3 <- ggplot(data = countries) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3)
map_3
- Adjust theme to change overall format:
map_4 <- ggplot(data = countries) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3) +
theme_light()
map_4
- Remove uncategorized areas:
map_4 <- ggplot(data = countries[countries$income_grp!="-99",]) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3) +
theme_light()
map_4
- Layer is currently in WGS 1984 (latitude and longitude coordinates)
- Can change to a different CRS when plotting to improve output
- CRS options can be found at spatialreference.org
- Easiest to designate desired CRS with EPSG codes
- We’ll try Equal Earth, as an example:
map_5 <- ggplot(data = countries[countries$income_grp!="-99",]) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3) +
coord_sf(crs = st_crs(8857)) +
theme_light()
map_5
Have a look at some of the map projections listed on Wikipedia and see if you can find one on spatialreference.org to try. Play around with the following code and see what you can make:
map_6 <- ggplot(data = countries[countries$income_grp!="-99",]) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3) +
coord_sf(crs = st_crs(8857)) +
theme_light()
map_6
- Add titles, a scale bar, and a north arrow (not that you need one at the global scale - this is just for demonstration):
map_7 <- ggplot(data = countries[countries$income_grp!="-99",]) +
geom_sf(aes(fill=income_grp)) +
scale_fill_brewer("Income Grouping", type = "qual", palette = 3) +
coord_sf(crs = st_crs(8857)) +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) + #This line adds a scale
#bar. The location argument shows where to place it (br = bottom right;
#tr = top right; bl = bottom left; tl = top left). The
#pad_ arguments show how far you want it from the corner of the map.
#The unit() function specifies what measurement units to use for the
#offset.
annotation_north_arrow(location = "bl", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) + #Adds the north arrow.
#The syntax is pretty similar to the scalebar function.
labs(title = "Country Income Groupings",
subtitle = "Data from Natural Earth") +
theme_minimal()
map_7Scale on map varies by more than 10%, scale bar may be inaccurate
Notice the warning about the scale bar. Why is that showing up, and what is it telling you?
Reading in and plotting vector data in the sf package
sf package also allows you to read in vector data files: geopackage (gpkg), shapefile (shp), etc.
Setting the Working Directory
R reads in data using information on the paths to files on your computer. It’s helpful to set the working directory, which givesRa default folder for reading and writing data for the session.
You can set the working directory to the folder where your script and data are located using the below code. But how do you find the PATH_TO_FILE? What I recommend for beginners is to save yourRscript and all the data you are using into a single folder somewhere on your computer. Then, if you are in RStudio, you can click Session -> Set Working Directory -> To Source File Location. RStudio will automatically generate a line of code like the one below. Then you can copy and paste that code into your script:
setwd("G:/My Drive/My Classes/World Cities/Labs/Geospatial_in_R_Jerash_Angkor")Jerash
Graeco-Roman city with considerable surviving ruins, in present-day Jordan. Co-exists with a contemporary city across a rainy season riverbed called a wadi. Subject to significant archaeological excavations over the past century and more.
Reading in an external vector data file:
jerash <- st_read("jerash.gpkg")Reading layer `jerash_interpretation_20180102' from data source
`G:\My Drive\My Classes\World Cities\Labs\Geospatial_in_R_Jerash_Angkor\jerash.gpkg'
using driver `GPKG'
Simple feature collection with 2501 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 770027.2 ymin: 3572942 xmax: 773161.7 ymax: 3577973
Projected CRS: ETRS89 / UTM zone 36NjerashSimple feature collection with 2501 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 770027.2 ymin: 3572942 xmax: 773161.7 ymax: 3577973
Projected CRS: ETRS89 / UTM zone 36N
First 10 features:
id Source Class IntBy LastObs FirstObs Conf
1 NA 2015 LiDAR / AP FB DS 2015-03-27 <NA> 2
2 0 2015 LiDAR / AP WALL DS 2015-03-27 <NA> 0
3 NA 2015 LiDAR / AP WALL DS 2015-03-27 <NA> 0
4 NA 2015 LiDAR / AP BR DS 2015-03-27 <NA> 2
5 NA 2015 LiDAR / AP WALL DS 2015-03-27 <NA> 0
6 NA 2015 LiDAR / AP WALL DS 2015-03-27 <NA> 3
7 NA 2015 LiDAR / AP WALL DS 2015-03-27 <NA> 0
8 NA 2015 LiDAR / AP BR DS 2015-03-27 <NA> 0
9 NA 2015 LiDAR / AP BR DS 2015-03-27 <NA> 1
10 NA 2015 LiDAR / AP BR DS 2015-03-27 <NA> 1
geom
1 MULTIPOLYGON (((771952.1 35...
2 MULTIPOLYGON (((772096.9 35...
3 MULTIPOLYGON (((771967.3 35...
4 MULTIPOLYGON (((772010.4 35...
5 MULTIPOLYGON (((772131.3 35...
6 MULTIPOLYGON (((772720.7 35...
7 MULTIPOLYGON (((772530.8 35...
8 MULTIPOLYGON (((771938.1 35...
9 MULTIPOLYGON (((771961.7 35...
10 MULTIPOLYGON (((772029.3 35...Notice that this layer is projected:
st_crs(jerash)Coordinate Reference System:
User input: ETRS89 / UTM zone 36N
wkt:
PROJCRS["ETRS89 / UTM zone 36N",
BASEGEOGCRS["ETRS89",
ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
MEMBER["European Terrestrial Reference Frame 1989"],
MEMBER["European Terrestrial Reference Frame 1990"],
MEMBER["European Terrestrial Reference Frame 1991"],
MEMBER["European Terrestrial Reference Frame 1992"],
MEMBER["European Terrestrial Reference Frame 1993"],
MEMBER["European Terrestrial Reference Frame 1994"],
MEMBER["European Terrestrial Reference Frame 1996"],
MEMBER["European Terrestrial Reference Frame 1997"],
MEMBER["European Terrestrial Reference Frame 2000"],
MEMBER["European Terrestrial Reference Frame 2005"],
MEMBER["European Terrestrial Reference Frame 2014"],
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[0.1]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4258]],
CONVERSION["UTM zone 36N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",33,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
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["Engineering survey, topographic mapping."],
AREA["Europe between 30°E and 36°E: Finland - onshore and offshore; Norway including Svalbard - onshore and offshore."],
BBOX[61.73,30,84.7,36]],
ID["EPSG",25836]]The most important variable for us is Class which denotes the type of feature:
table(jerash$Class)
BR EW FB MOD NAT RB RD STA WALL Water
1534 8 88 1 1 484 42 215 11 117 Here’s what these mean:
| Value | Summary |
|---|---|
| BR | Building remains |
| RB | Rubble |
| FB | Field boundary |
| STA | Standing building |
| WALL | City wall |
| RD | Road/Street |
| NAT | Natural |
| Water | Water management |
Unfortunately, some of the codes in the data are not defined in the data’s readme file.
Let’s separate these out into different objects that we can layer on the map separately:
jerash_building <- jerash %>% #Pipeline operator
filter(Class %in% c("BR", "RB", "STA"))
jerash_wall <- jerash %>%
filter(Class == "WALL")
jerash_water <- jerash %>%
filter(Class == "Water")Now we can make a plot that draws these objects as separate layers:
jerash_map_1 <- ggplot() +
geom_sf(data = jerash_water, fill="lightblue") +
geom_sf(data = jerash_building, fill="gray90") +
geom_sf(data = jerash_wall, fill="gray20") +
annotation_scale(location = "bl", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal()
jerash_map_1
A bit ugly; we can clip to the bounding box of the wall:
jerash_map_2 <- ggplot() +
geom_sf(data = jerash_water, fill="lightblue") +
geom_sf(data = jerash_building, fill="gray90") +
geom_sf(data = jerash_wall, fill="gray20") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal() +
coord_sf(xlim = st_bbox(jerash_wall)[c("xmin","xmax")],
ylim = st_bbox(jerash_wall)[c("ymin", "ymax")])
jerash_map_2
Hard to tell difference b/n features b/c of lines; try turning lines same colors as fill:
jerash_map_3 <- ggplot() +
geom_sf(data = jerash_water, fill="lightblue", color="lightblue") +
geom_sf(data = jerash_building, fill="gray90", color="gray90") +
geom_sf(data = jerash_wall, fill="gray20", color="gray20") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal() +
coord_sf(xlim = st_bbox(jerash_wall)[c("xmin","xmax")],
ylim = st_bbox(jerash_wall)[c("ymin", "ymax")])
jerash_map_3
The lat/long aren’t very informative:
jerash_map_5 <- ggplot() +
geom_sf(data = jerash_water, fill="lightblue", color="lightblue") +
geom_sf(data = jerash_building, fill="gray90", color="gray90") +
geom_sf(data = jerash_wall, fill="gray20", color="gray20") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal() +
#Removing grid lines and axis labels:
theme(panel.grid = element_blank(),
axis.text = element_blank()) +
coord_sf(xlim = st_bbox(jerash_wall)[c("xmin","xmax")],
ylim = st_bbox(jerash_wall)[c("ymin", "ymax")])
jerash_map_5
Finish this off by adding some meaningful labels and your own stylistic choices.
Adding context with maptiler
maptilerallows us to take snapshots of online maps (called tiles) and use them as the basemap (contextual bottom layer upon which a map can be constructed)get_tiles()takes bothsfandterraobjects and returns basemaps for the corresponding areasmaptilerhas tons of source map options, which you can see by looking athelp("get_tiles")
jerash_terrain <- get_tiles(jerash_wall,
provider = "Esri.WorldImagery")
jerash_terrainclass : SpatRaster
dimensions : 523, 526, 3 (nrow, ncol, nlyr)
resolution : 4.04, 4.04 (x, y)
extent : 771356.6, 773481.6, 3573851, 3575964 (xmin, xmax, ymin, ymax)
coord. ref. : ETRS89 / UTM zone 36N (EPSG:25836)
source(s) : memory
colors RGB : 1, 2, 3
names : Esri.World~50_13276_1, Esri.World~50_13276_2, Esri.World~50_13276_3
min values : 4, 15, 0
max values : 255, 255, 253 - This is a SpatRaster, which is the raster data format used by the
terrapackage. That means we can useterraandtidyterrafunctions to add this layer to our maps
jerash_map_6 <- ggplot() +
#If we're using an image, we use spatraster_rgb:
geom_spatraster_rgb(data = jerash_terrain) +
geom_sf(data = jerash_water, fill="lightblue", color="lightblue") +
geom_sf(data = jerash_building, fill="gray90", color="gray90") +
geom_sf(data = jerash_wall, fill="gray20", color="gray20") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal() +
#Removing grid lines and axis labels:
theme(panel.grid = element_blank(),
axis.text = element_blank()) +
coord_sf(xlim = st_bbox(jerash_wall)[c("xmin","xmax")],
ylim = st_bbox(jerash_wall)[c("ymin", "ymax")])
jerash_map_6
Using ESRI’s global imagery data, we can clearly see the modern city across the wadi from the ruins of Jerash.
Try some other tile sources and color adjustments to see if you can make these images clearer.
Greater Angkor
You’re already fairly familiar with Greater Angkor, and we’ll be spending more time with it next week, but I can’t pass up a gratuitous temple picture:
- Data are point locations of temples across the Greater Angkor complex
- In a csv file with longitude and latitude coordinates
- Need to read in csv and then tell
Rwhich columns have coordinate information to convert tosfobject:
#Function to load a CSV file; note the stringsAsFactors bit - it
#makes sure`R`doesn't automatically recode text, which can mess things
#up later:
angkor <- read.csv("carleton_Angkor_temples.csv",
stringsAsFactors = FALSE)
head(angkor) id morph azimuth area trait_1 trait_2 trait_3 trait_4 trait_5 trait_6
1 876 square 90.42425 90.8603 FALSE FALSE TRUE FALSE TRUE TRUE
2 874 square 90.00000 121.1024 FALSE FALSE FALSE FALSE TRUE TRUE
3 878 square 90.00000 144.9094 FALSE FALSE TRUE FALSE TRUE TRUE
4 933 square 90.00000 182.4583 FALSE FALSE TRUE FALSE FALSE FALSE
5 973 square 90.65866 185.5123 FALSE FALSE NA NA NA NA
6 968 square 89.53753 205.1487 FALSE FALSE NA NA NA NA
trait_7 trait_8 name date
1 FALSE FALSE Unknown 995
2 FALSE FALSE Sâk Krâop (Pr.) 889
3 FALSE FALSE Tuol 889
4 FALSE FALSE Top (Pr.) 1115
5 NA NA Angkor Thom (Laterite slabs and sculpture debris) 1002
6 NA NA Angkor Thom (terrace) 1115
dating_notes
1 Date derived from art historical analyses of lintels (Polkinghorne 2007).
2 Temple is included in CCC in Period 2.
3 Temple is included in CCC in Period 2.
4 Temple is included in the CCC
5 Temple is included in the CCC during Period 3.
6 Temple is included in the CCC for Period 4.
xlong ylat date_type xeast ynorth date_emp
1 103.8582 13.42621 empirical 376391.1 1484555 995
2 NA NA empirical 376029.9 1484559 889
3 NA NA empirical 376568.1 1484474 889
4 NA NA empirical 377453.2 1486464 1115
5 NA NA empirical 377851.6 1486810 1002
6 NA NA empirical 377220.5 1487538 1115- All temples have some sort of UTM coordinate, but we don’t know which UTM zone; not all have longitude and latitude, but looking at the code for the paper, the authors also dropped all without longitude and latitude, so we can do the same
angkor <- angkor %>%
# ! = "not", is.na() returns a TRUE for any missing values, so
# !is.na() means "not missing":
filter(!is.na(xlong)) %>%
# this next line tells`R`to create an sf points object and where
# to look for the points (the x-coordinate must always be first):
st_as_sf(coords = c("xlong", "ylat"))
angkorSimple feature collection with 1247 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 103.6134 ymin: 13.20946 xmax: 104.2757 ymax: 13.80604
CRS: NA
First 10 features:
id morph azimuth area trait_1 trait_2 trait_3 trait_4 trait_5
1 876 square 90.42425 90.8603 FALSE FALSE TRUE FALSE TRUE
2 965 square 89.91356 215.0931 FALSE FALSE NA NA NA
3 956 square 87.16414 226.2985 FALSE FALSE NA NA NA
4 1015 square 89.72978 229.4486 TRUE FALSE TRUE FALSE FALSE
5 1069 square 91.18823 231.9674 FALSE FALSE TRUE FALSE TRUE
6 940 square 88.94559 247.1555 FALSE FALSE NA NA NA
7 958 square 82.44541 259.3759 FALSE FALSE NA NA NA
8 542 square 100.91121 259.9270 TRUE FALSE FALSE FALSE TRUE
9 912 square 82.87515 264.0127 FALSE FALSE TRUE FALSE FALSE
10 832 square 87.17983 287.4387 FALSE FALSE FALSE FALSE FALSE
trait_6 trait_7 trait_8 name date
1 TRUE FALSE FALSE Unknown 995.000
2 NA NA NA Angkor Thom (Terrace) 926.000
3 NA NA NA Angkor Thom (Remains of a prasat) 926.000
4 TRUE FALSE FALSE Unknown 958.000
5 TRUE FALSE FALSE Dot Sdach Kamlong (Pr.) 986.000
6 NA NA NA Terrace I 926.000
7 NA NA NA Angkor Thom (Terrace R) 926.000
8 TRUE FALSE FALSE Srah Khout 1132.733
9 FALSE FALSE FALSE Preah Pithu, Y 982.000
10 TRUE FALSE FALSE Ak Yum (Pr.) 1001.000
dating_notes
1 Date derived from art historical analyses of lintels (Polkinghorne 2007).
2 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
3 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
4 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
5 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
6 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
7 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
8 Date derived through multiple linear regression (Klassen et al. 2018).
9 Date derived through graph-based semi-supervised machine learning (Klassen et al. 2018).
10 The temple was still active into the 10th century based on inscriptions found here (K. 752). We place an endpoint on the habitation in this area with the construction of the West Baray in the 11th century (Pottier et al. 2001b). The population associated with CCC is reduced to zero and temple communites are instead counted independently as part of the AMA.
date_type xeast ynorth date_emp geometry
1 empirical 376391.1 1484555 995 POINT (103.8582 13.42621)
2 graphbased 377181.0 1487114 NA POINT (103.8654 13.44916)
3 graphbased 375082.2 1487158 NA POINT (103.8461 13.45018)
4 graphbased 376487.6 1484022 NA POINT (103.8591 13.42129)
5 graphbased 406681.0 1496672 NA POINT (104.1031 13.56599)
6 graphbased 375715.8 1485894 NA POINT (103.8518 13.43813)
7 graphbased 376009.9 1487201 NA POINT (103.854 13.45042)
8 regression 383400.9 1491293 NA POINT (103.9227 13.48739)
9 graphbased 376725.1 1487109 NA POINT (103.8612 13.44932)
10 empirical 367571.5 1484426 1001 POINT (103.7768 13.42456)Now we have an sf object we can plot!
angkor_map_1 <- ggplot(data = angkor) +
geom_sf(aes(color = date)) +
# This next line is setting a color palette for a continuous scale;
# before we were just working with categories:
scale_color_viridis_c("Estimated Date") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal()
angkor_map_1Using plotunit = 'm'
Map scale is in centimeters?! Why?
st_crs(angkor)Coordinate Reference System: NANeed to tell it it’s in WGS 1984! Look it up on spatialreference.org!
# If we put st_crs() on the left side of the assignment operator,
# we can specify the CRS (but this isn't reprojecting the object!)
st_crs(angkor) <- st_crs(4326)
angkor_map_2 <- ggplot(data = angkor) +
geom_sf(aes(color = date)) +
# This next line is setting a color palette for a continuous scale;
# before we were just working with categories:
scale_color_viridis_c("Estimated Date") +
annotation_scale(location = "br", pad_x = unit(0.25, "in"),
pad_y = unit(0.25, "in")) +
theme_minimal()
angkor_map_2
Well, that’s a lot bigger than Jerash, but we need some more context!
Finish up the map with a suitable basemap, titles, etc.