-
Notifications
You must be signed in to change notification settings - Fork 2
/
My R Cookbook
139 lines (111 loc) · 5.29 KB
/
My R Cookbook
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
131
132
133
134
135
136
137
###This is my personal R-cookBook :) ###
##crop - pacote raster
r <- raster(nrow=45, ncol=90)
r[] <- 1:ncell(r)
e <- extent(-160, 10, 30, 60)
rc <- crop(r, e)
# use row and column numbers:
rc2 <- crop(r, extent(r, 5, 10, 7, 15))
# crop Raster* with Spatial* object
b <- as(extent(6, 6.4, 49.75, 50),'SpatialPolygons')
crs(b) <- crs(r)
rb <- crop(r, b)
# crop a SpatialPolygon* object with another one
if (require(rgdal) & require(rgeos)) {
p <- shapefile(system.file("external/lux.shp", package="raster"))
pb <- crop(p, b)
}
##maxima verossimilhanca para um modelo linear##
#gerando os dados
x = 1:200
y = 1 + 2*x + rnorm(200,0,40)
#construindo o modelo a ser ajustado
LL = function(a,b,mu,sigma){ #parametros do modelo e parametros de residuo
R = y - (a + b*x) #equacao do modelo
R = suppressWarnings(dnorm(R,mu,sigma)) #probabilidade de achar um valor de R para um determinado conjunto de parametros
-sum(log(R)) #soma do log das probs.
}
#algoritmo para estimativa de maxima verossimilhanca
fit = mle(LL, start=list(a=1,b=1,mu=0,sigma=10),fixed=list(mu=0),nobs=length(y)) #obs: mantendo mu=0 fixo, pois pois e do residuo
#tirando a prova
linTest = lm(y~x)
glmTest = glm(y~x)
summary(linTest)
summary(glmTest)
summary(fit)
##usando o pacote raster para calcular slope do terreno a a partir de uma imagem raster de elevacao do terreno (i.e. um DEM).
#tres funcoes estao dispsonivel no pacote raster: a de Evans (1980), Zevenbergen and Tore (1987) e a de Moore et al. (1993).
#respectivamente: evans(), zev.tho(), moore()
#fonte: https://www.r-bloggers.com/terrain-attributes-with-the-raster-package/
#Exemplo:
DTM = raster("caminho do raster")
slope = DEMderiv(DTM,"slope","evans") #(1)raster; (2)objetivo; (3)metodo
##calculando area usando raster e shape
#1)shape file
#get shapefile of Georgia’s borders
geo_border=getData(‘GADM’, country=’GEO’, level=0)
#calculate area [m2] of the polygon
sqm<-areaPolygon(geo_border)
#convert sqm to km2
sqkm<-sqm/1000000
#print area of Georgia according to shapefile, rounded to one digit
print(paste(“Area of Georgia (shapefile):”,round(sqkm, digits=1),”km2″))
#2)raster file
##getting SRTM data of Georgia
geo_raster = getData(‘alt’, country=’GEO’, mask=TRUE)
#get sizes of all cells in raster [km2]
cell_size<-area(geo_raster, na.rm=TRUE, weights=FALSE)
#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]
#compute area [km2] of all cells in geo_raster
raster_area<-length(cell_size)*median(cell_size)
#print area of Georgia according to raster object
print(paste(“Area of Georgia (raster):”,round(raster_area, digits=1),”km2″))
#3)calculando area de uma porcao especifica do raster file
##lowland zone 0-999 m
geo_raster1000 <- geo_raster
geo_raster1000[geo_raster1000 <=-1] <- NA
geo_raster1000[geo_raster1000 >999] <- NA
#calculate area of regions under 0 m asl
#get sizes of all cells under 0 m
cell_size<-area(geo_raster1000, na.rm=TRUE, weights=FALSE)
#delete NAs from vector of all raster cells
##NAs lie outside of the rastered region, can thus be omitted
cell_size<-cell_size[!is.na(cell_size)]
#compute area [km2] of all cells in geo_raster1000
lowland_area<-length(cell_size)*median(cell_size)
#print area of Georgia according to raster object
print(paste(“Area of lowland regions (0-999 m):”,round(lowland_area, digits=1),”km2″))
#plot lowland zone
X11(width=8, height=5)
plot(geo_raster1000,width=15, height=10,main=”Georgia lowland areas”,sub=paste(“0-999 meters asl, area =”,round(lowland_area, digits=1),”km2″))
plot(geo_border,add=TRUE)
## construindo uma equação empírica para modelar dados (i.e. fazer previsoes) ##
sp.data = read.csv("/home/anderson/PosDoc/dados_ocorrencia/PO_unique/Lagostomus maximus.csv",h=T) #abrindo tabela de ocorrencia de sp
var.pres = extract(predictors,sp.data[,2:3]) #OBS: predictors e um stack de rasters de variaveis climaticas no presente
var.pres = var.pres[complete.cases(var.pres),]
var.pass = extract(predictorsProjection,sp.data[,2:3]) #OBS: predictors e um stack de rasters de variaveis climaticas no passado
var.pass = var.pass[complete.cases(var.pass),]
modelo = glm(var.pres[,1]~var.pres[,2]) #esse modelo e so um exemplo (nao significa nada)
plot(var.pres[,1]~var.pres[,2]) #dados
points(fitted(modelo)~var.pres[,2],pch=20) #valores modelados para o presente (i.e. para os dados conhecidos)
points(fitted(modelo)~var.pass[,2],pch=1,col='blue') #valores modelados para o passado (i.e. para os dados desconhecidos)
qqplot(proj,fitted(modelo)) #comparando a distribuicao dos dados para presente e passado
points(proj,proj,type='l',col='red') #linha de referencia (linha de 45 graus)
##deixando dois rasters com mesma extent, origin e resolution##
##30/set/2017
library(raster)
newLayer = projectRaster(rasterProblema, rasterReferencia)
##excluir da amostra aqueles pontos que estiverem fora do poligono do shapefile
##30/set/2017
##fonte: http://robinlovelace.net/r/2014/07/29/clipping-with-r.html
library(rgdal)
stations <- readOGR("data", "lnd-stns") #pontos
zones <- readOGR("data", "london_sport") #poligono
stations <- spTransform(stations, CRS(proj4string(zones))) # transform CRS
plot(zones)
points(stations)
stations_subset <- stations[zones, ] ##proceidmento de cliping with poligon
plot(zones)
points(stations_subset)