Repository navigation
Expand file tree
/
Copy pathPoints in grid.R
More file actions
130 lines (89 loc) · 4.68 KB
/
Copy pathPoints in grid.R
File metadata and controls
130 lines (89 loc) · 4.68 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
# Original author: Jim Regetz
# Modified for Mosquito cases by: Jose Luis Delgado Davara - @jldelda
#
#********************************************************************************
#*******************************PART 1*******************************************
#********************************************************************************
# POINTS IN POLYGON
# Script to know in which cell of a grid, each point belongs to.
#
# Input Data
# * CSV table of all the cases of Zika, Dengue and Chikungunya in Río de Janeiro.
# Each row is a patient with its residencial geo-referenciated with lon and lat.
# lat-lon coordinate pair. These are derived from occurrence data
#
# * Polygon shapefile containing a grid of 1200 square metres.
# It contains 2651 cells starting from 0. The is of each cell is the
# number of it. Cell #0 is the one in the upper left corner, and cell #2651
# is the one in the right down corner. This file was created with Qgis.
# Qgis > Vector > Research Tools > Vector Grid > Set same coordenates with CRS.
#
# Workflow
# 1. Read and clean the Mosquito Borne Disease dataset and tell R to treat
# it as a set of spatial points.
# 2. Read in grid polygon of Río de Janeiro.
# 3. Identify which points lie within the cells.
# 4. For each point, get the name of the containing cell (if any), and
# add it to the Mosquito Borne Disease data table.
# 5. Write results to file, in both CSV and ESRI Shapefile formats,
# and draw a map of the Mosquitos sightings and cells.
#
# Output
# * CSV table similar to the input dataset, but with an additional
# column specifying the cell (if any) in which the patient lives.
# * Point shapefile identical to the CSV, but in a format more amenable
# to direct manipulation in a GIS.
# * Map visualization.
require(sp)
require(rgdal)
require(maps)
# read in Mosquito data
Mosquito <- read.csv("Mosquito.Borne.Disease.csv")
#Clean missing values
Mosquito.clean <- Mosquito[!is.na(Mosquito$latitude), ]
Mosquito.clean <- Mosquito.clean[!is.na(Mosquito.clean$longitude), ]
Mosquito.clean <- as.data.frame(Mosquito.clean)
# Turn it into a SpatialPointsDataFrame
coordinates(Mosquito.clean) <- c("longitude", "latitude")
# read in Grid polygons (For this you need the shp, dbf, prj, qpj and shx)
Grid <- readOGR(".", "Grid1")
# tell R that Mosquito coordinates are in the same lat/lon reference system
# as the parks data -- BUT ONLY BECAUSE WE KNOW THIS IS THE CASE!
proj4string(Mosquito.clean) <- proj4string(Grid)
# combine is.na() with over() to do the containment test; note that we
# need to "demote" Grid to a SpatialPolygons object first
inside.Grid <- !is.na(over(Mosquito.clean, as(Grid, "SpatialPolygons")))
# use 'over' again, this time with Grid as a SpatialPolygonsDataFrame
# object, to determine which Cell (if any) contains each sighting, and
# store the Cell id as an attribute of the Mosquito data
Mosquito.clean$Grid <- over(Mosquito.clean, Grid)$id
# Now your Mosquito.clean dataset has a new column with the id of the cell.
#********************************************************************************
#*******************************PART 2*******************************************
#********************************************************************************
# Quick analysis of the dataset
# Count number of point inside and outside the Grid
#Outside the grid
dim(Mosquito.clean[!inside.Grid, ])
# 1305 points
dim(Mosquito.clean[inside.Grid, ])
# 108851
# what fraction of sightings were inside the whole Grid?
mean(inside.Grid)
# [1] 0.9881532. 98% of the points are inside the Grid.
# 2% is from outside Río de Janeiro.
#Our points of interes are those inside the grid
Mosquito <- Mosquito.clean[inside.Grid, ]
# write the augmented Mosquito dataset to CSV
write.csv(Mosquito, "Mosquito-by-cells.csv", row.names=FALSE)
#********************************************************************************
#*******************************PART 3*******************************************
#********************************************************************************
# CREATE DATASET FOR ML ALGORITHM
# With this script we are going to
# 1.- create an empty cell dataset with the time periods.
# 2.- aggregate and count the number of occurrences each cell per time period.
# 3.- merge both dataset to have the number of occurrences in each cell.
# 1.- create an empty cell dataset with the time periods.
# 2.- aggregate and count the number of occurrences each cell per time period.
library(dplyr)