Guest User

Untitled

a guest
Jun 25th, 2018
85
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 3.92 KB | None | 0 0
  1. # arc de crise Afrique
  2. ged171 <- read.csv("data/ged171.csv")
  3.  
  4. #selection pays arc de crise Afrique.
  5. nigeria<- ged171[ged171$country=="Nigeria",]
  6. mali <- ged171[ged171$country=="Mali",]
  7. niger <- ged171[ged171$country=="Niger",]
  8. chad <- ged171[ged171$country=="Chad",]
  9. sudan <- ged171[ged171$country=="Sudan",]
  10. south_sudan <- ged171[ged171$country=="South Sudan",]
  11. ethiopia <- ged171[ged171$country=="Ethiopia",]
  12. eritrea<- ged171[ged171$country=="Eritrea",]
  13. libya<- ged171[ged171$country=="Libya",]
  14. algeria<- ged171[ged171$country=="Algeria",]
  15. egypt<- ged171[ged171$country=="Egypt",]
  16. somalia<- ged171[ged171$country=="Somalia",]
  17.  
  18. #fusion des lignes de chaque tableaux pour créer un tableau unique.
  19. arc_de_crise <- rbind(niger,mali,nigeria,chad,sudan,south_sudan,ethiopia,eritrea,
  20. libya,algeria,egypt,somalia)
  21.  
  22.  
  23. # load Packages
  24. library(sf)
  25. library(sp)
  26. library(cartography)
  27. library(rnaturalearth)
  28. library("SpatialPosition")
  29.  
  30.  
  31. #importation du fond de carte
  32. country <- ne_countries(scale = 110, type = "sovereignty", returnclass = "sf")
  33. afr_ucdp <- ne_countries(continent = 'africa',scale = 110, type = "sovereignty", returnclass = "sf")
  34. projrob <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
  35. afr_ucdp <- st_transform(x = afr_ucdp, crs = projrob)
  36. projrob <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
  37. country <- st_transform(x = country, crs = projrob)
  38.  
  39.  
  40. # Tranformation du data.frame nigeria en couche de points (objet sf)
  41.  
  42.  
  43. arc_de_crise <- st_as_sf(arc_de_crise, coords = c("longitude", "latitude"),remove = F,
  44. crs = 4326)
  45.  
  46. # Changer la projection de la couche des conflits wgs84 => robinson
  47.  
  48. arc_de_crise<- st_transform(x = arc_de_crise, crs = projrob)
  49. par(mar = c(0,0,1.2,0))
  50.  
  51. # Affichage des points sur la carte. (carte semi de points)
  52.  
  53. opar <- par(mar = c(0,0,1.2,0))
  54. plot(st_geometry(afr_ucdp),add=F, col = NA, border = NA, bg = "lightblue")
  55. plot(st_geometry(afr_ucdp), col = "grey90", border = NA, add=T)
  56. plot(st_geometry(afr_ucdp), add=T, col ="lightyellow", lwd = 0.5, border = "grey")
  57. plot(st_geometry(arc_de_crise), col="#96000040", cex = .5, pch = 16, add=TRUE)
  58.  
  59.  
  60. # pour travailler sur le nombre d'evenements
  61. arc_de_crise$cpt <- 1
  62.  
  63.  
  64. # Créer la grille de points
  65. mygrid <- CreateGrid(w = as(afr_ucdp,"Spatial"), resolution = 100000)
  66.  
  67. #years <- list(c(1989:1995),c(1996:2002), c(2003:2009), c(2010:2016))
  68.  
  69. # selectionner les evt des 2013 à 2016
  70. evtx <- arc_de_crise[arc_de_crise$year %in% c(1989:2016),]
  71. mymat <- CreateDistMatrix(knownpts = as(evtx,"Spatial"), unknownpts = mygrid)
  72.  
  73. # calcul des potentiels
  74. mystewart <- stewart(knownpts = as(evtx,"Spatial"), unknownpts = mygrid,
  75. matdist = mymat, varname = "best",
  76. typefct = "exponential", span = 75000,
  77. beta = 2, mask = as(afr_ucdp, "Spatial"))
  78.  
  79. # transformation en Raster
  80. ras <- rasterStewart(mystewart)
  81.  
  82. # chouix des bornes de classes
  83.  
  84. bb <- c(1,2,5,10,20,50,100,200,500,1000, 2000, 5000, 10000, 20000, 50000)
  85. #ou bb<- c(25,50,100,200,500,1000, 2000, 5000, 10000, 20000, 50000)
  86.  
  87. # transformation du raster en polygones
  88. contourpoly <- rasterToContourPoly(r = ras,breaks = bb,
  89. mask = as(afr_ucdp, "Spatial"))
  90. plot(contourpoly)
  91.  
  92.  
  93. opar <- par(mar = c(0,0,1.2,0))
  94. plot(afr_ucdp$geometry, add=F, col = NA, border = NA, bg = "lightblue")
  95. plot(country$geometry, col = "grey90", border = NA, add=T)
  96. plot(afr_ucdp$geometry, add=T, col ="lightyellow", lwd = 0.5, border = "grey")
  97. choroLayer(spdf = contourpoly,
  98. col = carto.pal("wine.pal", length(bb)-1),
  99. df = contourpoly@data,
  100. var = "center", legend.pos = c(-1294561, -3027262 ),
  101. breaks = bb, border = NA,
  102. lwd = 0.2,add=T,
  103. legend.title.txt = "Nb de morts\ndans les conflits",
  104. legend.values.rnd = 2)
  105. plot(afr_ucdp$geometry, border = "grey", add=T, lwd = 0.2)
  106. layoutLayer(title = "1989:2016", tabtitle = T, frame = F,scale = NULL)
Add Comment
Please, Sign In to add comment