Phelan12
Charlie Phelan
2024-06-14
library(ggplot2)
library(terra)## terra 1.7.78library(RColorBrewer)
library(sf)## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUElibrary(viridis)## Loading required package: viridisLitePart 1
in the zipfile countries.zip is a global shapefile with various socio-economic indicators for different countries. load this file into R and make plots of any of the two following variables; Median GDP, Population, or Income Group. try different color scales and transformations of the data to get the most informative maps
- Median GDP: (gdp_md_est)
- Population: (pop_est)
- Income group: (income_grp)
setwd("C:/Users/phela/OneDrive/Documents/GEOG5680/Mod12/countries")
list.files()## [1] "countries.dbf" "countries.geojson.txt" "countries.shp"
## [4] "countries.shx"countries = st_read("countries.shp")## Reading layer `countries' from data source
## `C:\Users\phela\OneDrive\Documents\GEOG5680\Mod12\countries\countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 177 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
## CRS: NAhead(countries)## Simple feature collection with 6 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825
## CRS: NA
## SP_ID scalerank featurecla labelrank sovereignt sov_a3
## 1 0 1 Admin-0 country 3 Afghanistan AFG
## 2 1 1 Admin-0 country 3 Angola AGO
## 3 2 1 Admin-0 country 6 Albania ALB
## 4 3 1 Admin-0 country 4 United Arab Emirates ARE
## 5 4 1 Admin-0 country 2 Argentina ARG
## 6 5 1 Admin-0 country 6 Armenia ARM
## adm0_dif level type admin adm0_a3 geou_dif
## 1 0 2 Sovereign country Afghanistan AFG 0
## 2 0 2 Sovereign country Angola AGO 0
## 3 0 2 Sovereign country Albania ALB 0
## 4 0 2 Sovereign country United Arab Emirates ARE 0
## 5 0 2 Sovereign country Argentina ARG 0
## 6 0 2 Sovereign country Armenia ARM 0
## geounit gu_a3 su_dif subunit su_a3 brk_diff
## 1 Afghanistan AFG 0 Afghanistan AFG 0
## 2 Angola AGO 0 Angola AGO 0
## 3 Albania ALB 0 Albania ALB 0
## 4 United Arab Emirates ARE 0 United Arab Emirates ARE 0
## 5 Argentina ARG 0 Argentina ARG 0
## 6 Armenia ARM 0 Armenia ARM 0
## name name_long brk_a3 brk_name
## 1 Afghanistan Afghanistan AFG Afghanistan
## 2 Angola Angola AGO Angola
## 3 Albania Albania ALB Albania
## 4 United Arab Emirates United Arab Emirates ARE United Arab Emirates
## 5 Argentina Argentina ARG Argentina
## 6 Armenia Armenia ARM Armenia
## brk_group abbrev postal formal_en formal_fr note_adm0
## 1 <NA> Afg. AF Islamic State of Afghanistan <NA> <NA>
## 2 <NA> Ang. AO People's Republic of Angola <NA> <NA>
## 3 <NA> Alb. AL Republic of Albania <NA> <NA>
## 4 <NA> U.A.E. AE United Arab Emirates <NA> <NA>
## 5 <NA> Arg. AR Argentine Republic <NA> <NA>
## 6 <NA> Arm. ARM Republic of Armenia <NA> <NA>
## note_brk name_sort name_alt mapcolor7 mapcolor8 mapcolor9
## 1 <NA> Afghanistan <NA> 5 6 8
## 2 <NA> Angola <NA> 3 2 6
## 3 <NA> Albania <NA> 1 4 1
## 4 <NA> United Arab Emirates <NA> 2 1 3
## 5 <NA> Argentina <NA> 3 1 3
## 6 <NA> Armenia <NA> 3 1 2
## mapcolor13 pop_est gdp_md_est pop_year lastcensus gdp_year
## 1 7 28400000 22270 -99 1979 -99
## 2 1 12799293 110300 -99 1970 -99
## 3 6 3639453 21810 -99 2001 -99
## 4 3 4798491 184300 -99 2010 -99
## 5 13 40913584 573900 -99 2010 -99
## 6 10 2967004 18770 -99 2001 -99
## economy income_grp wikipedia fips_10 iso_a2
## 1 7. Least developed region 5. Low income -99 <NA> AF
## 2 7. Least developed region 3. Upper middle income -99 <NA> AO
## 3 6. Developing region 4. Lower middle income -99 <NA> AL
## 4 6. Developing region 2. High income: nonOECD -99 <NA> AE
## 5 5. Emerging region: G20 3. Upper middle income -99 <NA> AR
## 6 6. Developing region 4. Lower middle income -99 <NA> AM
## iso_a3 iso_n3 un_a3 wb_a2 wb_a3 woe_id adm0_a3_is adm0_a3_us adm0_a3_un
## 1 AFG 004 004 AF AFG -99 AFG AFG -99
## 2 AGO 024 024 AO AGO -99 AGO AGO -99
## 3 ALB 008 008 AL ALB -99 ALB ALB -99
## 4 ARE 784 784 AE ARE -99 ARE ARE -99
## 5 ARG 032 032 AR ARG -99 ARG ARG -99
## 6 ARM 051 051 AM ARM -99 ARM ARM -99
## adm0_a3_wb continent region_un subregion region_wb
## 1 -99 Asia Asia Southern Asia South Asia
## 2 -99 Africa Africa Middle Africa Sub-Saharan Africa
## 3 -99 Europe Europe Southern Europe Europe & Central Asia
## 4 -99 Asia Asia Western Asia Middle East & North Africa
## 5 -99 South America Americas South America Latin America & Caribbean
## 6 -99 Asia Asia Western Asia Europe & Central Asia
## name_len long_len abbrev_len tiny homepart geometry
## 1 11 11 4 -99 1 MULTIPOLYGON (((61.21082 35...
## 2 6 6 4 -99 1 MULTIPOLYGON (((16.32653 -5...
## 3 7 7 4 -99 1 MULTIPOLYGON (((20.59025 41...
## 4 20 20 6 -99 1 MULTIPOLYGON (((51.57952 24...
## 5 9 9 4 -99 1 MULTIPOLYGON (((-65.5 -55.2...
## 6 7 7 4 -99 1 MULTIPOLYGON (((43.58275 41...ggplot() +
geom_sf(data = countries, aes(fill = gdp_md_est)) +
scale_fill_viridis(option = "cividis") +
coord_sf(xlim = c(-10,175), ylim = c(-45,75)) +
ggtitle("gpd_md_test Per Country - Europe, Asia, Africa and Austraila") +
theme_bw()
ggplot() +
geom_sf(data = countries, aes(fill = pop_est)) +
scale_fill_viridis(option = "H") +
coord_sf(xlim = c(-160,-25), ylim = c(-50,75)) +
ggtitle("pop_est Per Country - the Americas") +
theme_bw()
Part2
Using the NDVI values you calculated in the raster section, calculate the median NDVI for Bethel Island and Oakley. (Optional - for a slightly more challenging exercise, extract all NDVI values for these two places and compare using a t -test.)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionsetwd("C:/Users/phela/OneDrive/Documents/GEOG5680/Mod12/ca_places")
list.files()## [1] "ca_places.dbf" "ca_places.prj"
## [3] "ca_places.shp" "ca_places.shx"
## [5] "LC08_044034_20170614_B4.tif" "LC08_044034_20170614_B5.tif"places = st_read("ca_places.shp")## Reading layer `ca_places' from data source
## `C:\Users\phela\OneDrive\Documents\GEOG5680\Mod12\ca_places\ca_places.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1618 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.2695 ymin: 32.53432 xmax: -114.229 ymax: 41.99317
## Geodetic CRS: WGS 84str(places)## Classes 'sf' and 'data.frame': 1618 obs. of 17 variables:
## $ STATEFP : chr "06" "06" "06" "06" ...
## $ PLACEFP : chr "66140" "14190" "16560" "65042" ...
## $ PLACENS : chr "02411785" "02409487" "02410240" "02411779" ...
## $ GEOID : chr "0666140" "0614190" "0616560" "0665042" ...
## $ GEOIDFQ : chr "1600000US0666140" "1600000US0614190" "1600000US0616560" "1600000US0665042" ...
## $ NAME : chr "San Fernando" "Cloverdale" "Cotati" "San Buenaventura (Ventura)" ...
## $ NAMELSAD: chr "San Fernando city" "Cloverdale city" "Cotati city" "San Buenaventura (Ventura) city" ...
## $ LSAD : chr "25" "25" "25" "25" ...
## $ CLASSFP : chr "C1" "C1" "C1" "C1" ...
## $ PCICBSA : chr "N" "N" "N" "Y" ...
## $ MTFCC : chr "G4110" "G4110" "G4110" "G4110" ...
## $ FUNCSTAT: chr "A" "A" "A" "A" ...
## $ ALAND : num 6149607 8710251 4854614 56652007 22354266 ...
## $ AWATER : num 0 0 8380 26982227 42498 ...
## $ INTPTLAT: chr "+34.2886523" "+38.7947947" "+38.3284920" "+34.2677796" ...
## $ INTPTLON: chr "-118.4362428" "-123.0147032" "-122.7100491" "-119.2542062" ...
## $ geometry:sfc_MULTIPOLYGON of length 1618; first list element: List of 1
## ..$ :List of 1
## .. ..$ : num [1:174, 1:2] -118 -118 -118 -118 -118 ...
## ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:16] "STATEFP" "PLACEFP" "PLACENS" "GEOID" ...b5 = rast("LC08_044034_20170614_B5.tif")
b4 = rast("LC08_044034_20170614_B4.tif")
ndvi = (b5-b4) / (b5+b4)
bethel = places %>%
dplyr::filter(NAME == "Bethel Island")
bethel## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -121.6763 ymin: 38.00883 xmax: -121.6072 ymax: 38.04844
## Geodetic CRS: WGS 84
## STATEFP PLACEFP PLACENS GEOID GEOIDFQ NAME
## 1 06 06210 02407834 0606210 1600000US0606210 Bethel Island
## NAMELSAD LSAD CLASSFP PCICBSA MTFCC FUNCSTAT ALAND AWATER
## 1 Bethel Island CDP 57 U1 N G4210 S 13599384 1468307
## INTPTLAT INTPTLON geometry
## 1 +38.0281297 -121.6298509 MULTIPOLYGON (((-121.6763 3...oakley = places %>%
dplyr::filter(NAME == "Oakley")
oakley## Simple feature collection with 1 feature and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -121.7564 ymin: 37.96878 xmax: -121.6226 ymax: 38.02108
## Geodetic CRS: WGS 84
## STATEFP PLACEFP PLACENS GEOID GEOIDFQ NAME NAMELSAD LSAD
## 1 06 53070 02411294 0653070 1600000US0653070 Oakley Oakley city 25
## CLASSFP PCICBSA MTFCC FUNCSTAT ALAND AWATER INTPTLAT INTPTLON
## 1 C1 N G4110 A 41091677 789638 +37.9961485 -121.6936235
## geometry
## 1 MULTIPOLYGON (((-121.7562 3...bethel_ndvi = extract(ndvi, bethel, fun = median) ## Warning: [extract] transforming vector data to the CRS of the rasterbethel_median_ndvi = bethel_ndvi[, 2]
oakley_ndvi = extract(ndvi, oakley, fun = median)## Warning: [extract] transforming vector data to the CRS of the rasteroakley_median_ndvi = oakley_ndvi[, 2]
bethel_median_ndvi## [1] 0.4180438oakley_median_ndvi## [1] 0.3005387Extra credit
extract all NDVI values for these two places and compare using a t-test
bethel_values <- extract(ndvi, bethel)[, 2]## Warning: [extract] transforming vector data to the CRS of the rasteroakley_values <- extract(ndvi, oakley)[, 2]## Warning: [extract] transforming vector data to the CRS of the rastert_test_result <- t.test(bethel_values, oakley_values)
print(t_test_result)##
## Welch Two Sample t-test
##
## data: bethel_values and oakley_values
## t = 45.979, df = 26353, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.07298078 0.07948010
## sample estimates:
## mean of x mean of y
## 0.4124472 0.3362167