require("gstat")
require('sp')
require("spatstat")
require("ncf")
require("spdep")
require("raster")
require("ggplot2")
require("ggmap")


birds <- readRDS("F:\\Flash Drive\\ESCI 597\\Homework\\Rmarkdown\\HW7\\birdRichnessMexico.rds")
mexMap <- get_googlemap(center = c(lon = -102, lat = 23), zoom = 5, maptype = 'terrain')

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=23,-102&zoom=5&size=640x640&scale=2&maptype=terrain&sensor=false

map1 <- ggmap(mexMap) 

map1 <- map1 + geom_point(data=data.frame(birds), aes(x=long, y=lat,color=nSpecies),size=6)

map1 <- map1 + labs(x="Longtitude",y="Latitude",title="Bird Species Richness")
map1

Here is a different look at the same map using the method from the muese map.

birds <- readRDS("birdRichnessMexico.rds")
data(birds)
#coordinates(birds) <- ~x+y
proj4string(birds) <- CRS("+init=epsg:28992")
bubble(birds, xcol="lat",ycol="long",zcol="nSpecies", maxsize = 1.75, main = "Bird Species Richness")

The following is a series of lagged scatter plots with decreasing sizes of distances. The distances of bins so to speak are 200 km, 150 km and 55 km from the top down. In the 55 km lagged scatter plot we see that that correlation becomes very strong. With out much further analysis it appears that many of the data points are much further than 55 km apart. Further, the correlations decay to less than an r of 0.1 by approximately 1000 km.

hscat(nSpecies~1,data=birds,breaks = seq(0,2000000,by=200000))

hscat(nSpecies~1,data=birds,breaks = seq(0,2000000,by=150000))

hscat(nSpecies~1,data=birds,breaks = seq(0,2000000,by=55000))

Here, I extract the coordinates from the birds data set and transform them into a neighbor object using the tri2nb function so it will be usable by the moran.test function.

coords <- coordinates(birds)
birds.nb <- tri2nb(coords, row.names=NULL)
## 
##      PLEASE NOTE:  The components "delsgs" and "summary" of the
##  object returned by deldir() are now DATA FRAMES rather than
##  matrices (as they were prior to release 0.0-18).
##  See help("deldir").
##  
##      PLEASE NOTE: The process that deldir() uses for determining
##  duplicated points has changed from that used in version
##  0.0-9 of this package (and previously). See help("deldir").
moran.test(birds$nSpecies, nb2listw(birds.nb, style="W"))
## 
##  Moran I test under randomisation
## 
## data:  birds$nSpecies  
## weights: nb2listw(birds.nb, style = "W")  
## 
## Moran I statistic standard deviate = 17.166, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.703649266      -0.005025126       0.001704339

I also wanted to explore the spline.correlate that is described in the homework. We see in the plot below that the correlations decay to zero at about a distance equivalent to 10\(^\circ\) which is about 1000 km as it turns out. This trend clearly matches that which is seen in the previous lagged scatter plots.

The ‘spline.correlog’ function is calculating Morans I at increasing distances and then plotting its results gives a nice correlellogram. The summary tells us that there is strong auto correlation with a relatively low variation (0.001704339) and somewhat large distances. The strength of the correlations in bird richness is high until about 500km if you translate the degrees of latitude or longitude to km. Given that birds are migratory and cover large territories and seek far and wide feeding grounds we expect these correlations to stay strong for relatively long distances. These analysis of corrections thus far not telling about spatial direction or the direction of strenght. Perhaps this could be analogous to the direction of a magnetic field. Its there and you can measure it but you can’t directly see it.

fit2 <- spline.correlog(x=birds$long, y=birds$lat, z=birds$nSpecies, resamp = 100, quiet = TRUE)
plot(fit2)

summary.spline.correlog(fit2)
## $call
## [1] "spline.correlog(x = birds$long, y = birds$lat, z = birds$nSpecies, "
## [2] "    resamp = 100, quiet = TRUE)"                                    
## 
## $estimate
##                 x        e         y
## estimate 9.920275 2.968102 0.8989139
## 
## $quantiles
##               x        e         y
## 0      6.757321 1.227406 0.5160518
## 0.025  7.161961 1.580178 0.5302334
## 0.25   9.587548 2.185836 0.7112004
## 0.5    9.928654 3.149392 0.8809220
## 0.75  10.354809 3.809466 0.9987259
## 0.975 11.521012 5.069922 1.3083744
## 1     11.671988 5.675113 1.5816945

The following is a variogram with the distance cut off at 3500 km.The sill is not very apparent and we don’t see a strong asymptotic behavior. No direction is considered in this function. Direction is considered and contrasted in the next section

require("gstat")
  birdVar <- variogram(birds$nSpecies~1, birds, cloud = FALSE, cutoff=3500000)
  plot(birdVar,pch=20,cex=1.5,col="red",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Bird Richness")

The following is a series of variograms exploring different angels of direction to get a sense of the anisotropy that is seemingly present and also tinkering with the distance cut off. When the direction 160\(^\circ\) counterclockwise from North is used. the semi variance doesn’t level off into a clear sill. As we increase the distances that are used we begin to see some of what is likely edge effects. The max distance used is 4000 km. The variograms are improved slightly when the direction 160\(^\circ\) counterclockwise from North is used. Deviations from that typically deviate from the expected shape of the variagram. When applying the direction 135\(^\circ\) counterclockwise from North is used we see something more of an apparent sill but again it is not a preferably asymptotic sill.

This

require("gstat")
  birdVar <- variogram(birds$nSpecies~1, birds, cloud = FALSE, alpha=c(0,45,90,160), cutoff=2500000)
  plot(birdVar,pch=20,cex=1.5,col="red",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Bird Richness")

      require("gstat")
  birdVar <- variogram(birds$nSpecies~1, birds, cloud = FALSE, alpha=c(0,45,90,160), cutoff=4000000)
  plot(birdVar,pch=20,cex=1.5,col="red",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Bird Richness")

  require("gstat")
  birdVar <- variogram(birds$nSpecies~1, birds, cloud = FALSE, alpha=c(0,45,90,135), cutoff=2000000)
  plot(birdVar,pch=20,cex=1.5,col="red",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Bird Richness")

  require("gstat")
  birdVar <- variogram(birds$nSpecies~1, birds, cloud = FALSE, alpha=c(0,45,90,120), cutoff=1000000)
  plot(birdVar,pch=20,cex=1.5,col="red",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)", main = "Bird Richness")