Set options and load packages
install.packages("rgdal")
install.packages("maptools")
install.packages("rworldmap")
install.packages("sp")
install.packages("maps")
install.packages("rworldmap")
install.packages("rgdal")
install.packages("rgeos")
install.packages("raster")
install.packages("plyr")
require(devtools)
install_github("riatelab/cartography")
install.packages("maptools")
install.packages("PBSmapping")
install.packages("ggmap")
install.packages("geosphere")
install.packages("RgoogleMaps")
install.packages("OpenStreetMap")
install.packages("tripack")
options(digits=3, width=70)
#install.packages("rgdal")
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.2-15, (SVN revision 691)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /usr/share/gdal/2.1
## 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: (autodetected)
## Linking to sp version: 1.2-5
library(maptools)
## Checking rgeos availability: TRUE
Proj <- CRS("+proj=longlat +datum=WGS84")
download.file("http://freakonometrics.free.fr/CAN1.zip",destfile="CAN.zip")
unzip("CAN.zip")
CAN_shp <- readShapePoly("CAN_adm1.shp", verbose=TRUE, proj4string=Proj)
## Warning: use rgdal::readOGR or sf::st_read
## Shapefile type: Polygon, (5), # of Shapes: 13
plot(CAN_shp)
new_CAN_shp <- spTransform(CAN_shp, CRS("+init=epsg:26978"))
plot(new_CAN_shp)
library(sp)
library(rworldmap)
## ### Welcome to rworldmap ###
## For a short introduction type : vignette('rworldmap')
data(countriesLow)
crs.WGS84 <- CRS("+proj=longlat +datum=WGS84")
countries.WGS84 <- spTransform(countriesLow, crs.WGS84)
plot(countries.WGS84,col="light green")
crs.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
countries.laea <- spTransform(countriesLow, crs.laea)
plot(countries.laea,col="light green")
library(maps)
data("world.cities")
world.cities2 = world.cities[world.cities$pop>100000,]
plot(world.cities2$lon,world.cities2$lat,pch=19,cex=.7,axes=FALSE,xlab="",ylab="")
loc="https://raw.githubusercontent.com/jjrom/hydre/master/mapserver/geodata/rivers/GRDC_687_rivers"
ext=c(".dbf",".prj",".sbn",".sbx",".shp",".shp.xml",".shx")
for(i in 1:length(ext)){
loc1=paste(loc,ext[i],"?raw=true",sep="")
loc2=paste("GRDC_687_rivers",ext[i],sep="")
download.file(loc1,loc2)
}
library(sp)
shape <- readShapeLines("GRDC_687_rivers.shp")
## Warning: use rgdal::readOGR or sf::st_read
plot(shape,col="blue")
#crs.laea <- CRS("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs")
#shape.nouveau = spTransform(shape,crs.laea)
#plot(shape.nouveau,col="blue")
#download.file("http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_roads.zip", "routes.zip")
unzip("routes.zip")
## Warning in unzip("routes.zip"): erreur 1 lors de l'extraction d'un
## fichier zip
shape <- readShapeLines("ne_10m_roads.shp")
## Warning: use rgdal::readOGR or sf::st_read
plot(shape,pch=19,lwd=.7)
download.file("http://freakonometrics.free.fr/Italy.zip",destfile="Italy.zip")
unzip("Italy.zip")
library(rgdal)
ita1<-readOGR(dsn="/home/arthur/Italy/",layer="ITA_adm1",verbose=FALSE)
plot(ita1,col="light green")
ita1@data[,"NAME_1"][1:6]
## [1] Abruzzo Apulia Basilicata Calabria
## [5] Campania Emilia-Romagna
## 20 Levels: Abruzzo Apulia Basilicata Calabria ... Veneto
pos<-which(ita1@data[,"NAME_1"]=="Sicily")
Poly_Sicily<-ita1[pos,]
plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)
plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)
abline(v=5:20,col="light blue")
abline(h=35:50,col="light blue")
axis(1)
axis(2)
pos<-which(ita1@data[,"NAME_1"] %in% c("Abruzzo","Apulia","Basilicata","Calabria","Campania","Lazio","Molise","Sardegna","Sicily"))
ita_south <- ita1[pos,]
ita_north <- ita1[-pos,]
plot(ita1)
plot(ita_south,col="red",add=TRUE)
plot(ita_north,col="blue",add=TRUE)
library(rgeos)
## rgeos version: 0.3-26, (SVN revision 560)
## GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246
## Linking to sp version: 1.2-5
## Polygon checking: TRUE
ita_s<-gUnionCascaded(ita_south)
ita_n<-gUnionCascaded(ita_north)
plot(ita1)
plot(ita_s,col="red",add=TRUE)
plot(ita_n,col="blue",add=TRUE)
plot(ita1,col="light green")
plot(Poly_Sicily,col="yellow",add=TRUE)
gCentroid(Poly_Sicily,byid=TRUE)
## SpatialPoints:
## x y
## 14 14.1 37.6
## Coordinate Reference System (CRS) arguments: +proj=longlat
## +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
points(gCentroid(Poly_Sicily,byid=TRUE),pch=19,col="red")
G <- as.data.frame(gCentroid(ita1,byid=TRUE))
plot(ita1,col="light green")
text(G$x,G$y,1:20)
c1 <-G[c(17,20,6,16,8,12,2),]
c2 <- G[c(13,11,1,5,3,4),]
L1<-SpatialLines(list(Lines(list(Line(c1)),"Line1")))
L2<-SpatialLines(list(Lines(list(Line(c2)),"Line2")))
L1@proj4string<-ita1@proj4string
L2@proj4string<-ita1@proj4string
cross_road <-gIntersection(L1,L2)
cross_road@coords
## x y
## 1 11.1 44.0
## 1 14.2 41.7
plot(ita1,col="light green")
plot(L1,col="red",lwd=2,add=TRUE)
plot(L2,col="blue",lwd=2,add=TRUE)
plot(cross_road,pch=19,cex=1.5,add=TRUE)
library(sp)
library(raster)
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/GBR_adm2.rds","GBR_adm2.rds")
UK=readRDS("GBR_adm2.rds")
UK@data[159,"HASC_2"]="GB.NR"
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/IRL_adm0.rds","IRL_adm0.rds")
IRL=readRDS("IRL_adm0.rds")
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")
plot(IRL,add=TRUE)
download.file("http://biogeo.ucdavis.edu/data/gadm2.8/rds/FRA_adm0.rds","FRA_adm0.rds")
FR=readRDS("FRA_adm0.rds")
plot(UK, xlim = c(-4,-2), ylim = c(50, 59), main="UK areas")
plot(IRL,add=TRUE)
plot(FR,add=TRUE)
loc="https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/12/EU-referendum-result-data.csv"
referendum=read.csv(loc,header=TRUE,dec=".",sep=",",stringsAsFactors = FALSE)
referendum=referendum[c(3,6,13,14)]
library(plyr)
##
## Attaching package: 'plyr'
## The following object is masked from 'package:maps':
##
## ozone
referendum=ddply(referendum,.(Region,HASC_code),summarise,Remain=sum(Remain),Leave=sum(Leave))
sum(referendum$Leave)/(sum(referendum$Leave)+sum(referendum$Remain))
## [1] 0.519
referendum=referendum[c(referendum$Region!="Northern Ireland"),]
referendum=referendum[c(referendum$HASC_code!="Gibraltar"),]
row.names(referendum)=seq(1,nrow(referendum),1)
leave_or_remain=cbind(referendum,"Brexit?"=0)
leave_or_remain[,"Brexit?"]=ifelse(leave_or_remain$Remain<leave_or_remain$Leave,rgb(1,0,0,.7),rgb(0,0,1,.7))
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Brexit"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Brexit"]="blue"}else{
map_data[i,"Brexit"]=as.character(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],'Brexit?'])
}
}
plot(UK, col = map_data$Brexit, border = "gray1", xlim = c(-4,-2), ylim = c(50, 59), main="How the UK has voted?", bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)
legend(-1,59,c("Leave","Remain"),fill=c("red","blue"),bty="n")
#require(devtools)
#install_github("riatelab/cartography")
library(cartography)
## Loading required package: sf
## Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
cols <- carto.pal(pal1 = "red.pal",n1 = 5,pal2="green.pal",n2=5)
map_data=data.frame(UK@data)
map_data=cbind(map_data,"Percentage_Remain"=0)
for (i in 1:nrow(map_data)){
if(map_data[i,"NAME_1"]=="Northern Ireland"){
map_data[i,"Percentage_Remain"]=55.78}else{ map_data[i,"Percentage_Remain"]=100*leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]/(leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Remain"]+leave_or_remain[leave_or_remain$HASC_code==map_data$HASC_2[i],"Leave"])
}
}
plot(UK, col = "grey", border = "gray1", xlim = c(-4,-2), ylim = c(50, 59),bg="#A6CAE0")
plot(IRL, col = "lightgrey", border = "gray1",add=TRUE)
plot(FR, col = "lightgrey", border = "gray1",add=TRUE)
choroLayer(spdf = UK,
df = map_data,
var = "Percentage_Remain",
breaks = seq(0,100,10),
col = cols,
legend.pos = "topright",
legend.title.txt = "",
legend.values.rnd = 2,
add = TRUE)
library (maptools)
library (rgdal)
library (PBSmapping)
##
## -----------------------------------------------------------
## PBS Mapping 2.70.4 -- Copyright (C) 2003-2017 Fisheries and Oceans Canada
##
## PBS Mapping comes with ABSOLUTELY NO WARRANTY;
## for details see the file COPYING.
## This is free software, and you are welcome to redistribute
## it under certain conditions, as outlined in the above file.
##
## A complete user guide 'PBSmapping-UG.pdf' is located at
## /home/arthur/R/x86_64-pc-linux-gnu-library/3.3/PBSmapping/doc/PBSmapping-UG.pdf
##
## Packaged on 2017-06-28
## Pacific Biological Station, Nanaimo
##
## All available PBS packages can be found at
## https://github.com/pbs-software
##
## To see demos, type '.PBSfigs()'.
## -----------------------------------------------------------
download.file("http://freakonometrics.free.fr/paris-cartelec.zip","paris-cartelec.zip")
unzip("paris-cartelec.zip")
paris <- readShapeSpatial ( "paris-cartelec.shp" )
## Warning: use rgdal::readOGR or sf::st_read
## Warning: use rgdal::readOGR or sf::st_read
## Warning in Polygon(coords = crds): less than 4 coordinates in polygon
## Warning in Polygon(coords = crds): less than 4 coordinates in polygon
## Warning in Polygon(coords = crds): less than 4 coordinates in polygon
## Warning in Polygon(coords = crds): less than 4 coordinates in polygon
## Warning in Polygon(coords = crds): less than 4 coordinates in polygon
proj4string ( paris ) <- CRS( "+init=epsg:2154" )
paris <- spTransform(paris , CRS( "+proj=longlat +ellps=GRS80" ) )
plot ( paris )
length(paris)
## [1] 850
poly_paris <- SpatialPolygons2PolySet(paris)
point=c(2.33,48.849)
plot(point,xlim=point[1]+c(-.01,.01),ylim=point[2]+c(-.01,.01))
plot(paris,add=TRUE)
plot(point,xlim=point[1]+c(-.01,.01),ylim=point[2]+c(-.01,.01))
plot(paris,add=TRUE)
library(rgeos)
library(geosphere)
idx=unique(poly_paris$PID)
for(i in idx){
ct=centroid(as.matrix(poly_paris[poly_paris$PID==i,c("X","Y")]))
text(ct[1],ct[2],i,col="blue")}
which.poly <- function ( point ) {
idx=unique(poly_paris$PID)
idx_valid=NULL
for ( i in idx ) {
pip=point.in.polygon (point[ 1 ] , point[ 2 ] ,
poly_paris [ poly_paris$PID==i , "X" ] ,
poly_paris [ poly_paris$PID==i , "Y" ] )
if ( pip >0) idx_valid=c ( idx_valid , i ) }
return ( idx_valid ) }
which.poly(point)
## [1] 86
library(ggmap)
## Loading required package: ggplot2
library(geosphere)
#(Hotel <- geocode("8 sinaran drive singapore"))
Hotel = c(103.845,1.320279)
#(Airport <- geocode("singapore airport"))
Airport = c(103.9915,1.36442)
distHaversine(Airport,Hotel, r=6378.137)
## [1] 17
#mapdist(as.numeric(Hotel),as.numeric(Airport))
Loc=data.frame(rbind(Airport, Hotel))
library(ggmap)
library(RgoogleMaps)
CenterOfMap <- (apply(Loc,2,mean))
names(CenterOfMap)=c("lon","lat")
SG <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "terrain", source = "google")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=terrain&language=en-EN&sensor=false
SGMap <- ggmap(SG)
SGMap
SGMap <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "roadmap")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false
ggmap(SGMap)
SGMap <- get_map(c(lon=CenterOfMap["lon"], lat=CenterOfMap["lat"]),zoom = 11, maptype = "satellite")
## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=1.34235,103.91825&zoom=11&size=640x640&scale=2&maptype=satellite&language=en-EN&sensor=false
ggmap(SGMap)