Patrón de puntos
Sebastian franco
8/2/2022
set.seed(1019093771)
seqlat=seq(from=-73.3,to=-73.25,by=0.001)
seqlong=seq(from=5.54,to=5.58,by=0.001)
latitude=sample(seqlat,size = 100,replace=TRUE)
longitud=sample(seqlong,size = 100,replace=TRUE)
xy=data.frame(x=longitud,y=latitude)
plot(xy$x,xy$y,xlab="x",ylab="y",main = "Nube de puntos")
SMI=sort.int(runif(100,0.7,0.95),partial=10)
NDVI=sort.int(rnorm(100,0.45,0.06),partial=10)
LST=sort.int(26*rbeta(100,shape1=0.87,shape2=0.91),partial=10)
crop=factor(ifelse(rgamma(n=100,rate=0.8,shape=0.5)<0.5,0,1))
df1=data.frame(xy,SMI,NDVI,LST,crop)- Eliminar puntos duplicados
- Convertir el data.frame xy en un objeto ppp
library(spatstat)
dataxy=ppp(xy$x,xy$y,xrange=c(min(xy$x),max(xy$x)),yrange=c(min(xy$y),max(xy$y)))
any(duplicated(dataxy)) #El conjunto presenta puntos duplicados.## [1] TRUEy=unique(dataxy);y #se borran los puntos duplicados.## Planar point pattern: 98 points
## window: rectangle = [5.54, 5.58] x [-73.3, -73.25] unitsany(duplicated(y)) #se comprueba que el conjunto y no tiene duplicados.## [1] FALSE- El documento en la página 40 muestra la forma de colocar marcas a los puntos. Cree un conjunto de datos llamado Y1 usando como marca la variable “crop”. No olvide usar df1 y eliminar duplicados.
y1=unique(ppp(xy$x,xy$y,xrange=c(min(xy$x),max(xy$x)),yrange=c(min(xy$y),
max(xy$y)),marks = df1$crop))
any(duplicated(y1)) #Mostramos que los datos y1 no tienen duplicados.## [1] FALSE- En la misma página 40 con la forma de colocar marcas a los puntos. Cree un conjunto de datos llamado Y2 usando como marca la variable “NDVI”. No olvide usar df1 y eliminar duplicados.
#Se crea y2 sin datos duplicados.
y2=unique(ppp(xy$x,xy$y,xrange=c(min(xy$x),max(xy$x)),yrange=c(min(xy$y),
max(xy$y)),marks = df1$NDVI))
any(duplicated(y2)) #Mostramos que los datos y2 no tienen duplicados.## [1] FALSE- Grafique con plot de R tanto Y, como Y1 y Y2. Use el argumento size de plot para ajustar el tamaño de los puntos [plot(Y1, size=0.002), por ejemplo] de modo que se note cierta variabilidad. Cambie los colores de los puntos.Poner título al gráfico.
plot(y,size=0.8,cols = "orange",main="Gráfico de dispersión de puntos de y") #Gráfica de y
plot(y1,size=0.8,cols = "violet",main="Grafico de dispersión de puntos de y1") #Gráfica de y1
plot(y2,size=0.0008,cols="blue",main="Grafico de dispersión de puntos de y2") #Grádica de y2
6.) El código siguiente permite graficar las diferentes marcas. Grafique las marcas faltantes todo en un solo gráfico.
dataxy3=ppp(xy$x,xy$y,xrange=c(min(xy$x),max(xy$x)),yrange=c(min(xy$y),max(xy$y)),
marks = df1[,3:6])
y3=unique(dataxy3)
par(mfrow=c(1,2))
#Gráfica con marca NDVI.
plot(y3,which.marks = "NDVI",size=0.0009,main = "Gráfica con marca NDVI")
#Gráfica con la marca crop.
plot(y3,which.marks="crop",size=0.7,main = "Gráfica con marca crop") 
par(mfrow=c(1,1))
#Gráfica con todas las marcas.
plot(y3,size=0.002,main = "Gráfica con todas las marcas")
- En la página 46 se muestra cómo crear la ventana de observación del patrón de puntos. Utilice la función owin para crear la ventana de observación del patrón de puntos.
Creemos una ventana de observación
ventana1=owin(c(5.546,5.558),c(-73.289,-73.277));ventana1## window: rectangle = [5.546, 5.558] x [-73.289, -73.277] units- Use la funciones previas para obtener los atributos asociados en el conjunto Y3.
npoints(y3) #Número de puntos.## [1] 100marks(y3) #enseña las marcas. ## SMI NDVI LST crop
## 1 0.7003165 0.3739830 0.1265497 1
## 2 0.7003747 0.3798445 0.3797227 0
## 3 0.7067153 0.3711108 0.4365587 0
## 4 0.7128075 0.3142128 0.3815548 0
## 5 0.7047073 0.3855915 0.6773283 1
## 6 0.7140928 0.3234469 0.7076408 0
## 7 0.7153972 0.3842781 0.9111761 1
## 8 0.7020227 0.3287184 1.0375211 0
## 9 0.7077328 0.3748065 1.1126307 0
## 10 0.7242838 0.3865961 1.3627091 0
## 11 0.7296341 0.4434096 3.3562513 0
## 12 0.7319490 0.4314491 1.5478080 0
## 13 0.7457754 0.4017037 2.3855098 1
## 14 0.7320049 0.4263004 3.3612009 0
## 15 0.7344072 0.4238689 4.8809750 1
## 16 0.7703627 0.4482110 3.9325893 0
## 17 0.7649898 0.4466739 3.7992952 0
## 18 0.7492195 0.4353484 2.1885546 1
## 19 0.7761838 0.4047428 1.4671126 1
## 20 0.7573370 0.4452199 7.2625367 0
## 21 0.7375233 0.4485198 12.5144810 0
## 22 0.7341695 0.4166358 17.6239062 0
## 23 0.7778402 0.4097973 9.6418889 0
## 24 0.7321527 0.4235036 7.1406967 1
## 25 0.7370721 0.4364198 16.0279909 1
## 26 0.7684460 0.4361302 9.7191305 1
## 27 0.7795580 0.4530981 16.8807168 1
## 28 0.7788749 0.4522561 9.3645360 0
## 29 0.8105033 0.3953779 19.4101567 0
## 30 0.7947396 0.4168373 6.4009503 0
## 31 0.7818783 0.3977001 9.7446449 0
## 32 0.7857513 0.4316305 6.4535195 0
## 33 0.7850290 0.4432195 12.5788287 0
## 34 0.7853501 0.4276869 12.3586475 1
## 35 0.8015467 0.4282852 12.1561370 0
## 36 0.8397781 0.4132762 18.9662134 0
## 37 0.8151492 0.4456396 12.0619919 0
## 38 0.8172903 0.4083254 8.6285218 0
## 39 0.8187871 0.4539371 13.9224548 0
## 40 0.8126998 0.4325954 13.6287986 1
## 41 0.8327393 0.4049181 12.5069824 0
## 42 0.8368972 0.4470381 9.0011605 0
## 43 0.8331060 0.3937042 14.6805212 1
## 44 0.8388723 0.4086306 19.3387994 0
## 45 0.8318273 0.4462569 13.0219701 0
## 46 0.8124585 0.4625120 11.2681518 0
## 47 0.8214975 0.4629162 10.9375906 0
## 48 0.8175047 0.4668809 15.2516592 1
## 49 0.8170245 0.4576802 5.2752427 0
## 50 0.8119136 0.4751325 18.8192274 0
## 51 0.8276453 0.4613167 15.9436019 0
## 52 0.8172499 0.4592499 13.4006411 0
## 53 0.8242347 0.4811090 5.6641643 0
## 54 0.9312632 0.4651838 6.2405634 0
## 55 0.8615634 0.4771018 17.9894842 0
## 56 0.8806234 0.4749217 16.3789806 1
## 57 0.9181670 0.4639923 9.5357516 0
## 58 0.9365731 0.4613304 19.2412623 0
## 59 0.8433337 0.4554698 12.2923053 0
## 60 0.8830226 0.4746682 13.0110663 0
## 61 0.8918438 0.4756508 18.3385562 0
## 62 0.9176971 0.4791754 14.4945889 1
## 63 0.9125944 0.4687186 9.8029344 0
## 64 0.9228547 0.4619359 6.4780725 1
## 65 0.8710444 0.4684945 5.6562794 1
## 66 0.9038192 0.4540035 13.7735573 1
## 67 0.9279473 0.4624590 9.4460530 0
## 68 0.9017287 0.4564109 13.9631210 0
## 69 0.8761597 0.4554806 14.9960455 1
## 70 0.8483555 0.4843592 6.5934764 0
## 71 0.9315989 0.4886965 19.4723050 0
## 72 0.8801244 0.4962105 20.7642691 0
## 73 0.9437523 0.5015097 21.6593215 0
## 74 0.8483180 0.5193189 21.9290657 1
## 75 0.8737233 0.5188195 20.8287891 1
## 76 0.8736008 0.5228016 22.3774975 0
## 77 0.8487480 0.5182990 20.5725925 0
## 78 0.8430211 0.5048396 19.5778783 1
## 79 0.8967701 0.4983458 21.0523718 0
## 80 0.8550617 0.5033632 20.9269331 1
## 81 0.8529713 0.4849543 20.6841693 0
## 82 0.9172109 0.5435368 20.6217952 1
## 83 0.9152510 0.5452659 21.3425125 1
## 84 0.9053281 0.4988086 21.0304511 0
## 85 0.9254509 0.5642902 22.5060250 0
## 86 0.9445492 0.4883584 20.2820008 1
## 87 0.9217566 0.5207169 19.4905041 0
## 88 0.8984132 0.4874660 23.6274196 0
## 89 0.8826963 0.5401072 24.5748269 1
## 90 0.8783752 0.5942982 25.6120607 0
## 91 0.9206606 0.5527431 25.4225967 1
## 92 0.8910429 0.5518161 24.8395426 1
## 93 0.8399429 0.4999890 25.1998956 0
## 94 0.9169917 0.5252800 24.9916013 0
## 95 0.9260182 0.4978491 25.1972704 0
## 96 0.8430462 0.5511001 23.6371497 0
## 97 0.9039299 0.4848599 24.1861358 0
## 98 0.8979095 0.4843733 24.1001887 0
## 99 0.8785747 0.5038417 24.4403340 1
## 100 0.9393632 0.4973090 22.5154271 0coords(y3) #enseña las coordenadas.## x y
## 1 5.559 -73.255
## 2 5.550 -73.293
## 3 5.548 -73.269
## 4 5.548 -73.251
## 5 5.572 -73.269
## 6 5.567 -73.299
## 7 5.547 -73.254
## 8 5.554 -73.291
## 9 5.554 -73.278
## 10 5.540 -73.297
## 11 5.569 -73.267
## 12 5.568 -73.269
## 13 5.576 -73.290
## 14 5.550 -73.250
## 15 5.572 -73.273
## 16 5.579 -73.273
## 17 5.552 -73.256
## 18 5.553 -73.300
## 19 5.546 -73.297
## 20 5.548 -73.254
## 21 5.580 -73.289
## 22 5.574 -73.287
## 23 5.551 -73.287
## 24 5.561 -73.261
## 25 5.578 -73.278
## 26 5.559 -73.270
## 27 5.568 -73.257
## 28 5.545 -73.265
## 29 5.554 -73.253
## 30 5.557 -73.282
## 31 5.570 -73.298
## 32 5.573 -73.293
## 33 5.573 -73.272
## 34 5.575 -73.296
## 35 5.566 -73.263
## 36 5.551 -73.256
## 37 5.551 -73.264
## 38 5.577 -73.253
## 39 5.551 -73.271
## 40 5.576 -73.274
## 41 5.543 -73.272
## 42 5.568 -73.288
## 43 5.552 -73.260
## 44 5.572 -73.290
## 45 5.571 -73.290
## 46 5.571 -73.291
## 47 5.561 -73.264
## 48 5.545 -73.293
## 49 5.556 -73.275
## 50 5.540 -73.269
## 51 5.573 -73.284
## 52 5.565 -73.300
## 53 5.561 -73.253
## 54 5.572 -73.293
## 55 5.540 -73.258
## 56 5.567 -73.258
## 57 5.549 -73.269
## 58 5.577 -73.292
## 59 5.580 -73.295
## 60 5.559 -73.272
## 61 5.572 -73.263
## 62 5.556 -73.297
## 63 5.543 -73.268
## 64 5.579 -73.262
## 65 5.574 -73.299
## 66 5.568 -73.262
## 67 5.542 -73.293
## 68 5.541 -73.288
## 69 5.570 -73.275
## 70 5.574 -73.288
## 71 5.563 -73.290
## 72 5.545 -73.255
## 73 5.552 -73.252
## 74 5.542 -73.275
## 75 5.566 -73.280
## 76 5.580 -73.263
## 77 5.568 -73.263
## 78 5.556 -73.292
## 79 5.558 -73.281
## 80 5.561 -73.261
## 81 5.540 -73.281
## 82 5.541 -73.270
## 83 5.568 -73.256
## 84 5.575 -73.251
## 85 5.552 -73.285
## 86 5.546 -73.282
## 87 5.547 -73.297
## 88 5.564 -73.256
## 89 5.553 -73.279
## 90 5.554 -73.289
## 91 5.554 -73.292
## 92 5.546 -73.288
## 93 5.555 -73.260
## 94 5.557 -73.297
## 95 5.574 -73.253
## 96 5.577 -73.280
## 97 5.541 -73.296
## 98 5.551 -73.274
## 99 5.550 -73.279
## 100 5.556 -73.297as.owin(y3) #enseña la ventana definida en la creación.## window: rectangle = [5.54, 5.58] x [-73.3, -73.25] unitsas.data.frame(y3) #enseña las coordenadas y las marcas creadas.## x y SMI NDVI LST crop
## 1 5.559 -73.255 0.7003165 0.3739830 0.1265497 1
## 2 5.550 -73.293 0.7003747 0.3798445 0.3797227 0
## 3 5.548 -73.269 0.7067153 0.3711108 0.4365587 0
## 4 5.548 -73.251 0.7128075 0.3142128 0.3815548 0
## 5 5.572 -73.269 0.7047073 0.3855915 0.6773283 1
## 6 5.567 -73.299 0.7140928 0.3234469 0.7076408 0
## 7 5.547 -73.254 0.7153972 0.3842781 0.9111761 1
## 8 5.554 -73.291 0.7020227 0.3287184 1.0375211 0
## 9 5.554 -73.278 0.7077328 0.3748065 1.1126307 0
## 10 5.540 -73.297 0.7242838 0.3865961 1.3627091 0
## 11 5.569 -73.267 0.7296341 0.4434096 3.3562513 0
## 12 5.568 -73.269 0.7319490 0.4314491 1.5478080 0
## 13 5.576 -73.290 0.7457754 0.4017037 2.3855098 1
## 14 5.550 -73.250 0.7320049 0.4263004 3.3612009 0
## 15 5.572 -73.273 0.7344072 0.4238689 4.8809750 1
## 16 5.579 -73.273 0.7703627 0.4482110 3.9325893 0
## 17 5.552 -73.256 0.7649898 0.4466739 3.7992952 0
## 18 5.553 -73.300 0.7492195 0.4353484 2.1885546 1
## 19 5.546 -73.297 0.7761838 0.4047428 1.4671126 1
## 20 5.548 -73.254 0.7573370 0.4452199 7.2625367 0
## 21 5.580 -73.289 0.7375233 0.4485198 12.5144810 0
## 22 5.574 -73.287 0.7341695 0.4166358 17.6239062 0
## 23 5.551 -73.287 0.7778402 0.4097973 9.6418889 0
## 24 5.561 -73.261 0.7321527 0.4235036 7.1406967 1
## 25 5.578 -73.278 0.7370721 0.4364198 16.0279909 1
## 26 5.559 -73.270 0.7684460 0.4361302 9.7191305 1
## 27 5.568 -73.257 0.7795580 0.4530981 16.8807168 1
## 28 5.545 -73.265 0.7788749 0.4522561 9.3645360 0
## 29 5.554 -73.253 0.8105033 0.3953779 19.4101567 0
## 30 5.557 -73.282 0.7947396 0.4168373 6.4009503 0
## 31 5.570 -73.298 0.7818783 0.3977001 9.7446449 0
## 32 5.573 -73.293 0.7857513 0.4316305 6.4535195 0
## 33 5.573 -73.272 0.7850290 0.4432195 12.5788287 0
## 34 5.575 -73.296 0.7853501 0.4276869 12.3586475 1
## 35 5.566 -73.263 0.8015467 0.4282852 12.1561370 0
## 36 5.551 -73.256 0.8397781 0.4132762 18.9662134 0
## 37 5.551 -73.264 0.8151492 0.4456396 12.0619919 0
## 38 5.577 -73.253 0.8172903 0.4083254 8.6285218 0
## 39 5.551 -73.271 0.8187871 0.4539371 13.9224548 0
## 40 5.576 -73.274 0.8126998 0.4325954 13.6287986 1
## 41 5.543 -73.272 0.8327393 0.4049181 12.5069824 0
## 42 5.568 -73.288 0.8368972 0.4470381 9.0011605 0
## 43 5.552 -73.260 0.8331060 0.3937042 14.6805212 1
## 44 5.572 -73.290 0.8388723 0.4086306 19.3387994 0
## 45 5.571 -73.290 0.8318273 0.4462569 13.0219701 0
## 46 5.571 -73.291 0.8124585 0.4625120 11.2681518 0
## 47 5.561 -73.264 0.8214975 0.4629162 10.9375906 0
## 48 5.545 -73.293 0.8175047 0.4668809 15.2516592 1
## 49 5.556 -73.275 0.8170245 0.4576802 5.2752427 0
## 50 5.540 -73.269 0.8119136 0.4751325 18.8192274 0
## 51 5.573 -73.284 0.8276453 0.4613167 15.9436019 0
## 52 5.565 -73.300 0.8172499 0.4592499 13.4006411 0
## 53 5.561 -73.253 0.8242347 0.4811090 5.6641643 0
## 54 5.572 -73.293 0.9312632 0.4651838 6.2405634 0
## 55 5.540 -73.258 0.8615634 0.4771018 17.9894842 0
## 56 5.567 -73.258 0.8806234 0.4749217 16.3789806 1
## 57 5.549 -73.269 0.9181670 0.4639923 9.5357516 0
## 58 5.577 -73.292 0.9365731 0.4613304 19.2412623 0
## 59 5.580 -73.295 0.8433337 0.4554698 12.2923053 0
## 60 5.559 -73.272 0.8830226 0.4746682 13.0110663 0
## 61 5.572 -73.263 0.8918438 0.4756508 18.3385562 0
## 62 5.556 -73.297 0.9176971 0.4791754 14.4945889 1
## 63 5.543 -73.268 0.9125944 0.4687186 9.8029344 0
## 64 5.579 -73.262 0.9228547 0.4619359 6.4780725 1
## 65 5.574 -73.299 0.8710444 0.4684945 5.6562794 1
## 66 5.568 -73.262 0.9038192 0.4540035 13.7735573 1
## 67 5.542 -73.293 0.9279473 0.4624590 9.4460530 0
## 68 5.541 -73.288 0.9017287 0.4564109 13.9631210 0
## 69 5.570 -73.275 0.8761597 0.4554806 14.9960455 1
## 70 5.574 -73.288 0.8483555 0.4843592 6.5934764 0
## 71 5.563 -73.290 0.9315989 0.4886965 19.4723050 0
## 72 5.545 -73.255 0.8801244 0.4962105 20.7642691 0
## 73 5.552 -73.252 0.9437523 0.5015097 21.6593215 0
## 74 5.542 -73.275 0.8483180 0.5193189 21.9290657 1
## 75 5.566 -73.280 0.8737233 0.5188195 20.8287891 1
## 76 5.580 -73.263 0.8736008 0.5228016 22.3774975 0
## 77 5.568 -73.263 0.8487480 0.5182990 20.5725925 0
## 78 5.556 -73.292 0.8430211 0.5048396 19.5778783 1
## 79 5.558 -73.281 0.8967701 0.4983458 21.0523718 0
## 80 5.561 -73.261 0.8550617 0.5033632 20.9269331 1
## 81 5.540 -73.281 0.8529713 0.4849543 20.6841693 0
## 82 5.541 -73.270 0.9172109 0.5435368 20.6217952 1
## 83 5.568 -73.256 0.9152510 0.5452659 21.3425125 1
## 84 5.575 -73.251 0.9053281 0.4988086 21.0304511 0
## 85 5.552 -73.285 0.9254509 0.5642902 22.5060250 0
## 86 5.546 -73.282 0.9445492 0.4883584 20.2820008 1
## 87 5.547 -73.297 0.9217566 0.5207169 19.4905041 0
## 88 5.564 -73.256 0.8984132 0.4874660 23.6274196 0
## 89 5.553 -73.279 0.8826963 0.5401072 24.5748269 1
## 90 5.554 -73.289 0.8783752 0.5942982 25.6120607 0
## 91 5.554 -73.292 0.9206606 0.5527431 25.4225967 1
## 92 5.546 -73.288 0.8910429 0.5518161 24.8395426 1
## 93 5.555 -73.260 0.8399429 0.4999890 25.1998956 0
## 94 5.557 -73.297 0.9169917 0.5252800 24.9916013 0
## 95 5.574 -73.253 0.9260182 0.4978491 25.1972704 0
## 96 5.577 -73.280 0.8430462 0.5511001 23.6371497 0
## 97 5.541 -73.296 0.9039299 0.4848599 24.1861358 0
## 98 5.551 -73.274 0.8979095 0.4843733 24.1001887 0
## 99 5.550 -73.279 0.8785747 0.5038417 24.4403340 1
## 100 5.556 -73.297 0.9393632 0.4973090 22.5154271 0marks(y3)=seq(1,100) #cambio de marcas por una definida arbitrariamente
#cambio de coordenadas por otras arbitrariamente escogidas
coords(y3)=matrix(seq(5.54,5.58,0.001),seq(-73.3,-73.25,0.001),ncol = 2,nrow = 100) ## Warning in matrix(seq(5.54, 5.58, 0.001), seq(-73.3, -73.25, 0.001), ncol =
## 2, : la longitud de los datos [41] no es un submúltiplo o múltiplo del número de
## filas [100] en la matriz## Warning: 100 points were rejected as lying outside the specified window#Creación de una ventana
y3[ventana1] #Usamos la ventana creada en el punto 7.## Marked planar point pattern: 0 points
## marks are numeric, of storage type 'integer'
## window: rectangle = [5.546, 5.558] x [-73.289, -73.277] units- Con la información de la página 58 y usando solo las coordenadas x e y originales eliminando los puntos duplicados, obtenga el gráfico de densidad y el de contornos para el patrón de puntos.
Debido a que ya se eliminaron los puntos duplicados del punto 1 y 2, se procede a obtener el gráfico de densidad y el de contornos para el patrón de puntos.
#Gráfico de densidad de las coordenadas x y y originales.
plot(density(y))
#Gráfico de contorno de las coordenadas x y y originales.
contour(density(y))
- Con la página 59 hacer el histograma de puntos para cada coordenada para ver si existe alguna concentración mayor de puntos en alguna latitud o longitud.
hist(xy$x,xlab = "x",ylab="Frecuencia",main = "Histograma de frecuencia de x",
col = "cadetblue3") # Histograma de la longitud.
hist(xy$y,xlab = "y",ylab="Frecuencia",main = "Histograma de frecuencia de y",
col = "burlywood1") #Histograma de la latitud.
Si analizamos la longitud no se ve gran diferencia con respecto a la frecuencia de estos puntos, el punto entre 5.56 y 5.565 es el que menos se repite.
Respecto a la latitud no notamos una concentración grande de puntos en alguna coordenada vemos que las barras difieren muy poco en altura.
- Con la página 59 hacer el gráfico de densidad de modo que se perciba en él el patrón de puntos original sin marcas y sin duplicados.
plot(density(y))
plot(y,add=TRUE)
- Construya un patrón de puntos solo para aquellos sitios donde estaba presente un cultivo (llamarlo Y4) [pruebe la función split para crop].
#se usa la función split para dividir en dos tablas el df1, según el si hay
#o no ausencia de cultivo.
dividir=split(df1,crop,drop = TRUE)
#Se toman los puntos donde hay ausencia de cultivos.
ausencia=dividir$"1";ausencia## x y SMI NDVI LST crop
## 1 5.559 -73.255 0.7003165 0.3739830 0.1265497 1
## 5 5.572 -73.269 0.7047073 0.3855915 0.6773283 1
## 7 5.547 -73.254 0.7153972 0.3842781 0.9111761 1
## 13 5.576 -73.290 0.7457754 0.4017037 2.3855098 1
## 15 5.572 -73.273 0.7344072 0.4238689 4.8809750 1
## 18 5.553 -73.300 0.7492195 0.4353484 2.1885546 1
## 19 5.546 -73.297 0.7761838 0.4047428 1.4671126 1
## 24 5.561 -73.261 0.7321527 0.4235036 7.1406967 1
## 25 5.578 -73.278 0.7370721 0.4364198 16.0279909 1
## 26 5.559 -73.270 0.7684460 0.4361302 9.7191305 1
## 27 5.568 -73.257 0.7795580 0.4530981 16.8807168 1
## 34 5.575 -73.296 0.7853501 0.4276869 12.3586475 1
## 40 5.576 -73.274 0.8126998 0.4325954 13.6287986 1
## 43 5.552 -73.260 0.8331060 0.3937042 14.6805212 1
## 48 5.545 -73.293 0.8175047 0.4668809 15.2516592 1
## 56 5.567 -73.258 0.8806234 0.4749217 16.3789806 1
## 62 5.556 -73.297 0.9176971 0.4791754 14.4945889 1
## 64 5.579 -73.262 0.9228547 0.4619359 6.4780725 1
## 65 5.574 -73.299 0.8710444 0.4684945 5.6562794 1
## 66 5.568 -73.262 0.9038192 0.4540035 13.7735573 1
## 69 5.570 -73.275 0.8761597 0.4554806 14.9960455 1
## 74 5.542 -73.275 0.8483180 0.5193189 21.9290657 1
## 75 5.566 -73.280 0.8737233 0.5188195 20.8287891 1
## 78 5.556 -73.292 0.8430211 0.5048396 19.5778783 1
## 80 5.561 -73.261 0.8550617 0.5033632 20.9269331 1
## 82 5.541 -73.270 0.9172109 0.5435368 20.6217952 1
## 83 5.568 -73.256 0.9152510 0.5452659 21.3425125 1
## 86 5.546 -73.282 0.9445492 0.4883584 20.2820008 1
## 89 5.553 -73.279 0.8826963 0.5401072 24.5748269 1
## 91 5.554 -73.292 0.9206606 0.5527431 25.4225967 1
## 92 5.546 -73.288 0.8910429 0.5518161 24.8395426 1
## 99 5.550 -73.279 0.8785747 0.5038417 24.4403340 1#se crea el objeto ppp con la marca crop.
y4=unique(ppp(ausencia$x,ausencia$y,xrange=c(min(ausencia$x),max(ausencia$x)),yrange=c(min(ausencia$y),
max(ausencia$y)),marks = ausencia[,6]))
#no hay duplicados
any(duplicated(y4)) ## [1] FALSE#Patrón de puntos
plot(y4,cols = "purple",main = "Patrón de puntos donde hay ausencia de cultivos")
- Grafique la densidad y los contornos para Y4.
#Gráfico de densidad de las coordenadas x y y con presencia de cultivos.
plot(density(y4))
#Gráfico de contorno de las coordenadas x y y con presencia de cultivos.
contour(density(y4))
14) busque la función para un cono 3D y construya una imagen con esta función. Haga lo mismo para un paraboloide cualquiera.
#cono
x=seq(-3,3,length=100) #dominio donde se hace la gráfica.
y=seq(-3,3,length=100) #rango donde se hace la gráfica.
cono=function(x,y) sqrt(x^2+y^2)
#valores de z dependiendo del cambio de x y y.
z=outer(x, y, cono)
#Gráfica 3d de datos tomados.
persp(x,y,z) 
#Paraboloide.
x=seq(-3,3,length=100) #dominio donde se hace la gráfica.
y=seq(-3,3,length=100) #rango donde se hace la gráfica.
parabola=function(x,y) x^2+y^2
#valores de z dependiendo del cambio de x y y.
z=outer(x, y, parabola) #valores de z dependiente del cambio.
#Gráfica 3d de datos tomados.
persp(x,y,z)
- Con el código anterior convierta en imagen el patrón pero grafique no la primera marca sino la que identifica los cultivos.
#Gráfica de la marca que identifica los cultivos.
a=as.im(y4) #Convertir en imagen el patrón
plot(a) #
plot(y4,add = TRUE)
- Aplique la función summary sobre Y3.
summary(y3)## Marked planar point pattern: 100 points
## Average intensity 50000 points per square unit
##
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
##
## Mark variables: SMI, NDVI, LST, crop
## Summary:
## SMI NDVI LST crop
## Min. :0.7003 Min. :0.3142 Min. : 0.1265 0:68
## 1st Qu.:0.7774 1st Qu.:0.4273 1st Qu.: 6.4719 1:32
## Median :0.8350 Median :0.4585 Median :13.5147
## Mean :0.8302 Mean :0.4576 Mean :13.3334
## 3rd Qu.:0.8931 3rd Qu.:0.4884 3rd Qu.:20.5849
## Max. :0.9445 Max. :0.5943 Max. :25.6121
##
## Window: rectangle = [5.54, 5.58] x [-73.3, -73.25] units
## (0.04 x 0.05 units)
## Window area = 0.002 square unitsNuevo conjunto de datos.
set.seed(1019093771)
seqlat=seq(from=0,to=2000,by=2) #cm
seqlong=seq(from=0,to=500,by=2)#cm
lat=sample(seqlat,size = 150,replace = TRUE)
lon=sample(seqlong,size=150,replace=TRUE)
xyn=data.frame(x=lon,y=lat)
#se crean otros valores para las variables SMI,NDVI,LST y crop, esta vez con 150 datos.
SMI=sort.int(runif(150,0.7,0.95),partial=10)
NDVI=sort.int(rnorm(150,0.45,0.06),partial=10)
LST=sort.int(26*rbeta(150,shape1=0.87,shape2=0.91),partial=10)
crop=factor(ifelse(rgamma(n=150,rate=0.8,shape=0.5)<0.5,0,1))
df2=data.frame(xyn,SMI,NDVI,LST,crop)
library(spatstat)
dataxyn=ppp(xyn$x,xyn$y,xrange=c(min(xyn$x),max(xyn$x)),yrange=(c(min(xyn$y),max(xyn$y))))
any(duplicated(dataxyn))## [1] FALSEyn=unique(dataxyn)
plot(yn,size=0.5)
summary(yn)## Planar point pattern: 150 points
## Average intensity 0.0001506024 points per square unit
##
## Coordinates are integers
## i.e. rounded to the nearest unit
##
## Window: rectangle = [0, 500] x [2, 1994] units
## (500 x 1992 units)
## Window area = 996000 square units- Del nuevo conjunto de datos Yn y usando la página 78 obtenga la densidad de puntos por cm2 y conviértala a m2.
#La densidad de metros por centimetros cuadrados estimada es 0.00015.
summary(yn)$intensity## [1] 0.0001506024Como se tiene la medida en \(cm^2\) para pasarla a metros debemos multiplicar la intensidad estimada por \(1000^2\) ya que 1 metro son 1000 centímetros.
#calculando la densidad de puntos por metros cuadrados estimada.
lambda=summary(yn)$intensity*(1000)^2;lambda ## [1] 150.6024- Use la página 80 y construya una cuadricula de 4 filas y 4 columnas para ver la distribución del conteo de puntos en cada celda.
#cuenta los puntos en cada cuadrícula.
quadratcount(yn, nx = 4, ny = 4) ## x
## y [0,125) [125,250) [250,375) [375,500]
## [1.5e+03,1.99e+03] 14 6 11 10
## [998,1.5e+03) 13 11 7 3
## [500,998) 10 7 9 10
## [2,500) 10 14 6 9#guardar los conteos como una matriz en el objeto Q.
Q <- quadratcount(yn, nx = 4, ny = 4)
#Gráfica de puntos superpuesta por las cuadrículas.
plot(yn,size=0.7)
plot(Q, add = TRUE)
- En la misma página 80 se muestra la construcción de la densidad del patrón y el gráfico 3d de esta densidad. Aplique sobre Yn.
densidad=density(yn) #densidad el patrón.
#gráfica de la densidad sobrepuesta por el patrón de puntos.
plot(densidad)
plot(yn,add=TRUE)
persp(densidad) #Gráfico 3D de la densidad.
- Use la pag 81 y construya el gráfico de contornos de Yn
#Gráfico de contorno
contour(densidad)
- Usando el código inicial de la pág 82 muestre el patrón de puntos y haga al lado los gráficos de las covariables NDVI, SMI y LST.
No se entiende este punto porque bei.extra es una matriz de 101x201 mientras que las covariables creadas por los datos son vectores, claramente no se puede seguir bien el ejemplo del libro.
#Covariable índice de vegetación de diferencia normalizado.
par(mfrow = c(1, 2))
plot(yn)
plot(NDVI)
par(mfrow = c(1, 1))
#Temperatura en la superficie del suelo.
par(mfrow = c(1, 2))
plot(yn)
plot(LST)
par(mfrow = c(1, 1))
#Existencia de vegetación.
par(mfrow = c(1, 2))
plot(yn)
plot(crop)
par(mfrow = c(1, 1))
#índice de humedad del suelo.
par(mfrow = c(1, 2))
plot(yn)
plot(SMI)
par(mfrow = c(1, 1))- En las pag 89 y 90 aparece una prueba de conteos de cuadrantes para completa aleatoriedad espacial del patrón. Use los cuadrantes antes construidos para probar la hipótesis de CSR (CAE en español)
Se seguirá la prueba de hipótesis \(\chi^2\) para comprobar aleatoriedad espacial completa donde.
\(H_0:\) Los cuadrantes son aleatorios (Distribuyen Poisson). \(H_0:\) Los cuadrantes no son aleatorios (no es un proceso Poisson).
#Prueba estadística
M=quadrat.test(yn, nx = 4, ny = 4);M #no se rechaza H0##
## Chi-squared test of CSR using quadrat counts
##
## data: yn
## X2 = 14.693, df = 15, p-value = 0.9474
## alternative hypothesis: two.sided
##
## Quadrats: 4 by 4 grid of tilesplot(yn,size=0.01)
plot(M, add = TRUE, cex = 0.9)
Como se puede ver el p-valor=0.9474 se rechaza la hipótesis nula, es decir que el patrón de puntos tiene completa aleatorización espacial.
- De la pág 91 use la prueba de Kolmogorov-Smirnov para CSR usando en primer lugar la dirección latitudinal y en segundo la longitudinal.
#En este código no funciona lo que dice en el libro, se intento crear una función
#como la dice el libro pero genera error
#t=function(x,y)x
#ks.test(yn,"x")
ks.test(yn$x,yn$x)## Warning in ks.test(yn$x, yn$x): p-value will be approximate in the presence of
## ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: yn$x and yn$x
## D = 0, p-value = 1
## alternative hypothesis: two-sidedSe hace la prueba de normalidad usando la dirección latitudinal y la dirección longitudinal, según el test de kolmogorov se tiene indicios de normalidad.
- ¿Tendrá el patrón de puntos alguna relación con alguna covariable? Use el código de las pág 92 y 93 y verifique si se puede ignorar la independencia y concluir que el patrón no es uniforme en su intensidad. Use una a una las covariables.
#No quedó clara la explicación de covariables para comparar por lo mencionado en
#el punto 21.
#b <- quantile(NDVI, probs = (0:4)/4)
#NDVIcut <- cut(NDVI, breaks = b, labels = 1:4)
#V <- tess(image = NDVIcut)
#quadrat.test(yn, tess = V)- Con el patrón Yn ajuste un modelo de Poisson Homogéneo (pag 97) . Escriba el modelo ajustado.
Se ajusta el modelo homogéneo de Poisson.
Mp=ppm(yn, ~1);Mp## Stationary Poisson process
## Intensity: 0.0001506024
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) -8.800867 0.08164966 -8.960898 -8.640837 *** -107.7882- Ajuste un modelo que dependa solo de la latitud (pag 97). Escriba el modelo.
Se calcula el modelo que depende sólo de la latitud.
My=ppm(yn~y) ;My## Nonstationary Poisson process
##
## Log intensity: ~y
##
## Fitted trend coefficients:
## (Intercept) y
## -8.788626e+00 -1.229333e-05
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -8.788626e+00 0.1630588042 -9.1082153673 -8.4690366001 ***
## y -1.229333e-05 0.0001420336 -0.0002906742 0.0002660875
## Zval
## (Intercept) -53.89850628
## y -0.08655221este modelo es \(\lambda_{\theta}((y))=\exp \left(-8.7886-0.0000123y\right)\).
- Ajuste un modelo que dependa solo de la longitud (pag 97). Escriba el modelo.
Se calcula el modelo que sólo depende de la longitud.
Mx=ppm(yn~x) ;Mx## Nonstationary Poisson process
##
## Log intensity: ~x
##
## Fitted trend coefficients:
## (Intercept) x
## -8.528357164 -0.001144169
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -8.528357164 0.1527885325 -8.82781719 -8.228897e+00 ***
## x -0.001144169 0.0005705162 -0.00226236 -2.597743e-05 *
## Zval
## (Intercept) -55.818045
## x -2.005497este modelo es \(\lambda_{\theta}((x))=\exp \left(-8.528-0.001144x\right)\).
- Ajuste un modelo en función de x+y . Escriba el modelo.
Se establece ajustar un modelo Poisson no homogéneo con una intensidad logarítica, el modelo es \(\lambda_{\theta}((x, y))=\exp \left(\theta_{0}+\theta_{1} x+\theta_{2} y\right)\).
Mb=ppm(yn, ~x + y);Mb## Nonstationary Poisson process
##
## Log intensity: ~x + y
##
## Fitted trend coefficients:
## (Intercept) x y
## -8.516282e+00 -1.144143e-03 -1.213304e-05
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -8.516282e+00 0.2079552724 -8.9238664955 -8.108697e+00 ***
## x -1.144143e-03 0.0005705170 -0.0022623357 -2.595024e-05 *
## y -1.213304e-05 0.0001420293 -0.0002905053 2.662393e-04
## Zval
## (Intercept) -40.9524681
## x -2.0054495
## y -0.0854263El modelo queda como \(\lambda_{\theta}((x, y))=\exp \left(-8.516-0.00114x+0.0000121y\right)\).
- Ajuste un modelo polinomial cuadrático. Escriba el modelo ajustado.
Para ajustar un modelo con intensidades constantes pero desiguales a cada lado de la línea vertical x = 300
Mpol=ppm(yn, ~polynom(x, y, 2));Mpol## Nonstationary Poisson process
##
## Log intensity: ~x + y + I(x^2) + I(x * y) + I(y^2)
##
## Fitted trend coefficients:
## (Intercept) x y I(x^2) I(x * y)
## -8.290324e+00 -1.673031e-03 -6.436028e-04 1.630735e-06 -2.575337e-07
## I(y^2)
## 3.461424e-07
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -8.290324e+00 3.791356e-01 -9.033416e+00 -7.547232e+00 ***
## x -1.673031e-03 2.387238e-03 -6.351932e-03 3.005870e-03
## y -6.436028e-04 5.944891e-04 -1.808780e-03 5.215745e-04
## I(x^2) 1.630735e-06 4.377496e-06 -6.948999e-06 1.021047e-05
## I(x * y) -2.575337e-07 9.360442e-07 -2.092147e-06 1.577079e-06
## I(y^2) 3.461424e-07 2.684892e-07 -1.800867e-07 8.723714e-07
## Zval
## (Intercept) -21.8663817
## x -0.7008227
## y -1.0826149
## I(x^2) 0.3725269
## I(x * y) -0.2751298
## I(y^2) 1.2892229side <- function(z) factor(ifelse(z < 300, "izquierda","derecha"))
ppm(bei, ~side(x))## Nonstationary Poisson process
##
## Log intensity: ~side(x)
##
## Fitted trend coefficients:
## (Intercept) side(x)izquierda
## -5.1222055 0.5283049
##
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## (Intercept) -5.1222055 0.02188965 -5.1651084 -5.0793026 *** -234.00128
## side(x)izquierda 0.5283049 0.03373948 0.4621768 0.5944331 *** 15.65836Se considera como primer nivel lo que se encuentra a izquiera de 300 y segundo lo que se encuentra a derecha así, el ajuste queda.
\[\lambda_{\theta}((x, y))= \begin{cases}\exp (-5.1222) & \text { if } x<300 \\ \exp (-5.1222055+0.5283049) & \text { if } x \geq 300\end{cases}\]
- Use el anova para comparar el modelo con solo covariable latitud y el modelo con x+y, Señale el modelo que muestra mejor ajuste.
anova(My) #análisis de la desviación para y.## Analysis of Deviance Table
## Terms added sequentially (first to last)
##
## Df Deviance Npar
## NULL 1
## y 1 0.0074914 2anova(Mb) #análisis de desviación para x+y.## Analysis of Deviance Table
## Terms added sequentially (first to last)
##
## Df Deviance Npar
## NULL 1
## x 1 4.0548 2
## y 1 0.0073 3Note que el modelo con sólo la covarialbe longitud tiene mejor ajuste porque su desviasión es menor a la del modelo x+y.
- Use el anova para comparar el modelo con x+y y el modelo cuadrático. Señale el modelo que muestra mejor ajuste.
anova(Mb)## Analysis of Deviance Table
## Terms added sequentially (first to last)
##
## Df Deviance Npar
## NULL 1
## x 1 4.0548 2
## y 1 0.0073 3anova(Mpol)## Analysis of Deviance Table
## Terms added sequentially (first to last)
##
## Df Deviance Npar
## NULL 1
## x 1 4.0548 2
## y 1 0.0073 3
## I(x^2) 1 0.1404 4
## I(x * y) 1 0.0816 5
## I(y^2) 1 1.6323 6Note que el modelo cuadrático tiene más desviasión, por lo tanto es el que tiene menor ajuste, para este caso el mejor modelo es x+y