help me to write r code


##############################################################################
# 1. map visualization: sf
##############################################################################

# for this lab we'll be using the following packages for static maps:
#install.packages('rgdal')
install.packages("sf")   # compile from source: no
install.packages("spData")

library(sf)          # classes and functions for vector data
library(spData)        # load geographic data


##### about the sf library & vector data:
# The geographic vector model is based on points located within a coordinate reference system. 
# Points can represent self-standing features (e.g. location) or they can be linked together to 
# form more complex geometries such as lines and polygons. 
# The coordinates of London, for example, can be represented by c(-0.1, 51.5). 
# This means that its location is -0.1 degrees east and 51.5 degrees north of the origin.
# Most point geometries contain only two dimensions, or they could include height to form a third dimension.

# The core geometry types are fully supported by the R package sf.
# sf can represent all common vector geometry types: points, lines, polygons and their 
# respective ‘multi’ versions (which group together features of the same type into a single feature).
# These vector data have discrete, well-defined borders, meaning that vector datasets 
# usually have a high level of precision.

# sf also supports geometry collections, which can contain multiple geometry types in a single object.

# The package documentaton can be accessed as follows: 

vignette('sf1')

# A vignette is a long-form guide to your package with function documentation. 
# A vignette is like a book chapter or an academic paper: 
# it can describe the problem that your package is designed to solve, and then show the reader how to solve it.

# As the first vignette explains, simple feature objects in R are stored in a data frame, 
# with geographic data occupying a special column, usually named ‘geom’ or ‘geometry’. 
# We will use the world dataset provided by the spData as an example.

# world is a spatial object containing spatial and attribute columns,
# check data type
class(world)  # both sf class and a data frame

# check column names 
names(world)

# the last column contains the geographic information.
# The contents of this geom column give sf objects their spatial powers: 
# world$geom is a ‘list column’ that contains all the coordinates of the country polygons.
world$geom 
# MULTIPOLYGON in the summary above refers to the geometry type of countries in the world object.


# you can check out columns in the data frame:
world$name_long
world$gdpPercap
# similar to regular data frames

# summarize information in a column
summary(world$gdpPercap)


# The sf package provides a plot() method for visualizing geographic data:
plot(world)

# Note that instead of creating a single map, the plot() command has created multiple maps, 
# one for each variable in the world datasets. 
# This is useful for exploring the spatial distribution of different variables 

# plot a single feature on map:
plot(world['gdpPercap'])  # GDP per capita
plot(world['lifeExp'])    # life expectancy

# plot a subset of the data
plot(world[3:6]) 

unique(world$continent)
# plot population on a single continent 
plot(world[world$continent=='Asia','pop'])

# Plots are added as layers to existing images by setting add = TRUE.
# We can plot the Asian continent over a map of the world.
# Note that the first plot must only have one facet for add = TRUE to work. 
# If the first plot has a key, reset = FALSE must be used
world_asia <- world[world$continent == "Asia", ]
# st_union return a single geometric object by combining several features; internal boundaries are ignored
asia <- st_union(world_asia)  
asia

plot(world["pop"], reset = F)  # reset = FALSE: keep the plot in a mode that allows adding further map elements
plot(asia, add = T, col = "red") # add = T: add more layers to the graph


# sf’s plot method also has arguments specific to geographic data. e.g. expandBB
# can be used plot an sf object in context: it takes a numeric vector of length four that expands the 
# bounding box of the plot relative to zero in the following order: bottom, left, top, right. 
# This is used to plot India in the context of its giant Asian neighbors, 
# with an emphasis on China to the east, in the following code chunk

india = world[world$name_long == "India", ]

plot(st_geometry(india), expandBB = c(0, 0.3, 0.4, 0.6), col = "gray", lwd = 3)
plot(world_asia[0], add = TRUE)

# Note the use of [0] to keep only the geometry column and lwd to emphasize India
# what happens if we increase the parameters in expandBB?


# There are various ways to modify maps with sf’s plot() method.
# e.g we can overlay circles, whose diameters are proportional to country populations, 
# on a map of the world. The circles are geographic centroids

cex <- sqrt(world$pop) /10000  # define circle size proportional to country population
world_cents <- st_centroid(world, of_largest = TRUE) 

# st_centroid() converts one geometry type (polygons) to points
# of_largest = TRUE: return centroid of the largest (sub)polygon of a MULTIPOLYGON 
# rather than of the whole MULTIPOLYGON

plot(world["continent"], reset = F)
plot(st_geometry(world_cents), add = T, cex = cex)
# st_geometry assigns geometry to a data.frame and creates an sf object. 
# here we're converting the data.frame of world_cents into points


##############################################################################
# 2. attribute manipulation
##############################################################################

# the world data includes geometric information and a data frame with attribute columns.
# this allows for standard data frame operations, such as merge, join, and select.

# select countries less than 10000 km2
sel_area <- world[world$area_km2 < 10000,]
sel_area$name_long # check country names

# subset() and filter() provide the same results:
sel_area <- subset(world, area_km2 < 10000) 
sel_area <- filter(world, area_km2 < 10000)

# now, select countries in Asia with populations larger than 100 mil
sel_area <- world %>% subset(continent == "Asia" & pop>=10^8)
plot(sel_area['pop'])

# generate population average by continent
# we'll use summaryBy from the doBy library
library(doBy)

# st_drop_geometry removes the geometry column and leaves only the data frame
# the geometry info can cause problems in data summary
dat <- data.frame(st_drop_geometry(world))
pop_ave <- summaryBy(pop ~ continent, data = dat, na.rm=T)
pop_ave

# the operation can also be achieved by a pipe operator. here we simply create a column of 
# population average for each continent:
world %>% 
  st_drop_geometry() %>% 
  dplyr::group_by(continent) %>%  # summarize population mean by continent groups
  dplyr::summarize(pop = mean(pop, na.rm = TRUE)) 

# additional attributes can be merged into the data from other sources
# here we introduce the coffee production data 
# "coffee_data" from the spData package
??coffee_data

# the coffee_data has 3 columns: name_long names major coffee-producing nations, coffee_production_2016 
# and coffee_production_2017 contain estimated values for coffee production in each year

# now we join the continent population average and coffee data with world. 
# we also create a new variable of population density:
world2 <- world %>% dplyr::left_join(pop_ave, by='continent') %>%
  dplyr::left_join(coffee_data, by='name_long') %>%
  dplyr::mutate(pop_dens = pop / area_km2)  # mutate creates a new column

names(world2)
# notice the new column, pop_dens, is at the second to the last position. 
# the last column is reserved for geometric info.

# plot population maps: continent population average, density
plot(world2["pop.mean"])
plot(world2["pop_dens"])

# coffee production in 2017
plot(world2["coffee_production_2017"])


####################################
# 3. spatial operations
####################################

# Spatial subsetting is the process of selecting features of a spatial object based on whether or not 
# they in some way relate in space to another object.
# An example of spatial subsetting is provided by the nz and nz_height datasets in spData. 
# These contain the 16 main regions and 101 highest points in New Zealand
# The following codes first creates an object representing Canterbury, 
# then uses spatial subsetting to return all high points in the region

canterbury <- nz %>% dplyr::filter(Name == "Canterbury")
canterbury_height <- nz_height[canterbury, ]
names(canterbury)
names(canterbury_height)
canterbury_height$elevation

# plot high points on Canterbury, New Zealand
plot(nz['Name'], col='white', reset=F, key.pos=NULL, main='')  # omit color keys
plot(st_union(canterbury), col = 'red', add=T)  # use st_union to consolidate the sf object. see above
plot(canterbury_height['elevation'], add=T)     # add high points in Canterbury


# suppose you want to find out the average height of high points in each region. 
# This is a good example of spatial aggregation: it is the geometry of the source 
# that defines how values in the target object are grouped. 
# to do this, you first join the nz data with height (nz_height),
# then summarize average elevation by regional names.

nz_avheight <- nz %>%
  st_join(nz_height) %>%
  dplyr::group_by(Name) %>%
  dplyr::summarize(elevation = mean(elevation, na.rm = TRUE))

names(nz_avheight)
plot(nz_avheight['elevation'])

# The resulting objects have the same geometry as the aggregating object nz but 
# with a new column representing the average height of each region.
# other summary functions such as median() and sd() can be used instead of mean()

# group sequence
sample(c(1, 2, 3, 5, 6, 7, 8, 9, 10, 16, 23))



####################################
# 4. week 10 Assignment: Due on Sunday midnight
####################################


# Task 3. Plot the population map of New Zealand, with colors for Marlborough, Auckland, and Gisborne only


# Task 4. Show geographical centroids on New Zealand, with size proportional to log population/10.
# color the map by median height of regions