R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

This document is to summarise some key features and geo-processing in Geocomputation with R

Packages used

library(sf)

package 㤼㸱sf㤼㸲 was built under R version 3.6.3Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

library(spData)

package 㤼㸱spData㤼㸲 was built under R version 3.6.3To access larger datasets in this package, install the spDataLarge package with:
`install.packages('spDataLarge', repos='https://nowosad.github.io/drat/', type='source')`

library(raster)

package 㤼㸱raster㤼㸲 was built under R version 3.6.3

library(ggplot2)

package 㤼㸱ggplot2㤼㸲 was built under R version 3.6.3RStudio Community is a great place to get help: https://community.rstudio.com/c/tidyverse

This section is dedicated to the use of sf package which covers a wide range of geo-processing tools for vector data. it provides a consistent and easy-to-use procedure to analyse and process vector data as well as working well with dataframe. If you like, please have a look here

vignette(package = “sf”) # see which vignettes are available vignette(“sf1”) # an introduction to the package

Loading and checking data

world data contain different features of the world country such as population, area and country names. This dataset is stored as sf object. If we want to convert sp or other shapefiles to sf, we can use function st_as_sf()

names(world) # Names of different countries 

 [1] "iso_a2"    "name_long" "continent" "region_un" "subregion" "type"      "area_km2"  "pop"       "lifeExp"   "gdpPercap"
[11] "geom"     

head(world) # Check first few rows of dataf 

Simple feature collection with 6 features and 10 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -180 ymin: -18.28799 xmax: 180 ymax: 83.23324
CRS:            EPSG:4326

plot(world)

plotting the first 9 out of 10 attributes; use max.plot = 10 to plot all

*Choosing to plot only population variable and Asia

sf_asia<-world[world$continent=="Asia",]

plot(sf_asia["pop"])

• st_union function is used to form a single boundary of an area. For example, we create asia

sf_asia_union<-st_union(sf_asia)

# Plot world map

plot(world["pop"], reset=F)

plot(sf_asia_union,add=T,col="2")

• Extract centroids from polygons

plot(world["continent"], reset = FALSE)
cex = sqrt(world$pop) / 10000
world_cents = st_centroid(world, of_largest = TRUE)

st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

plot(st_geometry(world_cents), add = TRUE, cex = cex)

Create points and dataframe using sf

lnd_point<-rbind(c(0.1,51.5),c(0.2,52.3))

lnd_point = st_multipoint(lnd_point)                 # sfg object

lnd_geom = st_sfc(lnd_point, crs = 4326)           # sfc object

lnd_attrib = data.frame(                           # data.frame object
  name = c("London","HN"), temperature = c(25,27),
  date = as.Date(c("2017-06-21","2018-02-01"))
  )
lnd_sf = st_sf(lnd_attrib, geometry = lnd_geom)    # sf object

row names were found from a short variable and have been discarded

lnd_sf

Simple feature collection with 2 features and 3 fields
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 0.1 ymin: 51.5 xmax: 0.2 ymax: 52.3
CRS:            EPSG:4326
    name temperature       date                       geometry
1 London          25 2017-06-21 MULTIPOINT ((0.1 51.5), (0....
2     HN          27 2018-02-01 MULTIPOINT ((0.1 51.5), (0....

# Coordinate Reference System CRS -Get all coordinate reference system

library(rgdal)

package 㤼㸱rgdal㤼㸲 was built under R version 3.6.3rgdal: version: 1.4-8, (SVN revision 845)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/DELL/Documents/R/win-library/3.6/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/DELL/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.4-1 

crs_data = rgdal::make_EPSG()

View(crs_data)

Get crs of the object

st_crs(world)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

If our object is missing crs or wrong, we can set crs with function st_set_crs(object,crs)

sf_newcrs<-st_set_crs(sf_asia_union,4326)

st_crs(sf_newcrs)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

library(ggplot2)
ggplot(data = world) + geom_sf(aes(fill = pop)) +scale_fill_viridis_c(option = "plasma", trans = "sqrt")

If we want to calculate the area of country or region, we can use st_area)(). This will result in square meters.

library(lwgeom)

package 㤼㸱lwgeom㤼㸲 was built under R version 3.6.3Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.6.1, PROJ 4.9.3

library(units)

package 㤼㸱units㤼㸲 was built under R version 3.6.3udunits system database from C:/Users/DELL/Documents/R/win-library/3.6/units/share/udunits

st_area(sf_asia_union)/1000000 # Square kilometer 

31252459 [m^2]

units::set_units(st_area(sf_asia_union), km^2)

31252459 [km^2]

If we want to label all country names, we can use text to name it

sf_points<-st_centroid(world, of_largest=T)

st_centroid assumes attributes are constant over geometries of xst_centroid does not give correct centroids for longitude/latitude data

sf_coord<-data.frame(st_coordinates(sf_points))

df<-data.frame(Names=world$name_long,x=sf_coord$X,y=sf_coord$Y)
               
plot(world["name_long"],main="World Map")

text(x=df[,2],y=df[,3],labels=df$Names)

Attribute data associated with sptial geometry

sf package supported following methods.

methods(class = "sf") # methods for sf objects, first 12 shown

 [1] $<-                        [                          [[<-                       aggregate                 
 [5] as.data.frame              cbind                      coerce                     dbDataType                
 [9] dbWriteTable               extent                     extract                    filter                    
[13] identify                   initialize                 mask                       merge                     
[17] plot                       print                      raster                     rasterize                 
[21] rbind                      select                     show                       slotsFromS3               
[25] st_agr                     st_agr<-                   st_area                    st_as_sf                  
[29] st_bbox                    st_boundary                st_buffer                  st_cast                   
[33] st_centroid                st_collection_extract      st_convex_hull             st_coordinates            
[37] st_crop                    st_crs                     st_crs<-                   st_difference             
[41] st_filter                  st_force_polygon_cw        st_geometry                st_geometry<-             
[45] st_interpolate_aw          st_intersection            st_intersects              st_is                     
[49] st_is_polygon_cw           st_join                    st_line_merge              st_linesubstring          
[53] st_m_range                 st_make_valid              st_minimum_bounding_circle st_nearest_points         
[57] st_node                    st_normalize               st_point_on_surface        st_polygonize             
[61] st_precision               st_reverse                 st_sample                  st_segmentize             
[65] st_set_precision           st_shift_longitude         st_simplify                st_snap                   
[69] st_snap_to_grid            st_split                   st_subdivide               st_sym_difference         
[73] st_transform               st_transform_proj          st_triangulate             st_union                  
[77] st_voronoi                 st_wrap_dateline           st_write                   st_z_range                
[81] st_zm                      transform                 
see '?methods' for accessing help and source code

Aggregation

-Nonsptial aggregation. Report the total of population in each continent

continent_pop<-aggregate(pop~continent, FUN=sum,data=world,na.rm=T)
continent_pop

-Spatial Aggregation. Report the total of population in each continent with spatial objects

continent_pop<-aggregate(world["pop"],by=list(world$continent),FUN=sum,na.rm=T)
continent_pop

Simple feature collection with 8 features and 2 fields
Attribute-geometry relationship: 0 constant, 1 aggregate, 1 identity
geometry type:  GEOMETRY
dimension:      XY
bbox:           xmin: -180 ymin: -90 xmax: 180 ymax: 83.64513
CRS:            EPSG:4326
                  Group.1        pop                       geometry
1                  Africa 1154946633 MULTIPOLYGON (((49.54352 -1...
2              Antarctica          0 MULTIPOLYGON (((-163.7129 -...
3                    Asia 4311408059 MULTIPOLYGON (((120.295 -10...
4                  Europe  669036256 MULTIPOLYGON (((-51.6578 4....
5           North America  565028684 MULTIPOLYGON (((-61.68 10.7...
6                 Oceania   37757833 MULTIPOLYGON (((169.6678 -4...
7 Seven seas (open ocean)          0 POLYGON ((68.935 -48.625, 6...
8           South America  412060811 MULTIPOLYGON (((-66.95992 -...

• A better way of aggregating the total of each continent population

continent_pop<- world %>% group_by(continent) %>% summarise(n=n(), total=sum(pop,na.rm=T))

Error in group_by(., continent) : could not find function "group_by"

• We can find the most populated continents

populated_continent<-world %>% dplyr::select(pop, continent) %>% group_by(continent) %>% summarize(pop = sum(pop, na.rm = TRUE), 
                                                                              n_countries = n()) %>%  top_n(n = 4, wt = pop) %>%
  arrange(desc(pop)) %>%  st_drop_geometry()
populated_continent

Vector attribute Joining

we can see the full list of joining functions vignette("two-table")

# Left join
world_coffee<-left_join(world,coffee_data, by="name_long")

plot(world_coffee["coffee_production_2016"])

• Create non-spatial object and force it to be spatial object

coffee_world<-left_join(coffee_data,world,by="name_long")

world_coffee_map<-st_as_sf(coffee_world) # Coarse the geometry columns to the sptial object. The coffee_data has only 47 observations
plot(world_coffee_map)

plotting the first 10 out of 12 attributes; use max.plot = 12 to plot all

Remove spatial informaiton or attributes

-This process is similar to dataframe, but it is noted that if we want to remove geometry column, we need to use st_drop_geometry()

-Creating a new spatial and give it a crs

st_crs(st_English)

Coordinate Reference System:
  User input: EPSG:27700 
  wkt:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","27700"]]

Identify the UTM code any points in the Planet

UTM_code = function(lonlat) {
  utm = (floor((lonlat[1] + 180) / 6) %% 60) + 1
  if(lonlat[2] > 0) {
    utm + 32600
  } else{
    utm + 32700
  }
}

Warning message:
In png_available() : internal error -3 in R_decompress1

# Type your long and lat and get your UTM EPSG
UTM_code(c(105.83,21.02)) # Hanoi Vietnam

[1] 32648

Writing shapefile or saving shapefiles

It can be saved in various formats ESRI Shapefile, KML, gpkg Geopackage

library(rgdal)
st_write(world,dsn = "WorldMap.gpkg", append = F)

Deleting layer `WorldMap' using driver `GPKG'
Writing layer `WorldMap' to data source `WorldMap.gpkg' using driver `GPKG'
Writing 177 features with 10 fields and geometry type Multi Polygon.