library(spatstat)
## Loading required package: nlme
##
## spatstat 1.45-0 (nickname: 'One Lakh')
## For an introduction to spatstat, type 'beginner'
##
## Note: spatstat version 1.45-0 is out of date by more than 4 months; we recommend upgrading to the latest version.
library(RColorBrewer)
library(sp)
library(raster)
##
## Attaching package: 'raster'
## The following objects are masked from 'package:spatstat':
##
## area, rotate, shift
## The following object is masked from 'package:nlme':
##
## getData
library(rgdal)
## rgdal: version: 1.1-8, (SVN revision 616)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/david.maupin/Desktop/Packages/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/david.maupin/Desktop/Packages/rgdal/proj
## Linking to sp version: 1.2-3
library(rgeos)
## rgeos version: 0.3-19, (SVN revision 524)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-3
## Polygon checking: TRUE
library(maptools)
## Checking rgeos availability: TRUE
library(openxlsx)
library(maps)
##
## # maps v3.1: updated 'world': all lakes moved to separate new #
## # 'lakes' database. Type '?world' or 'news(package="maps")'. #
library(knitr)
library(rmarkdown)
#The purpose of the file is to convert a data file into a PPP object using the Spatstat in order
#to analyze the spatial points. Once the points are analyzed the file can then be converted
#into a shapefile for further use in other programs (Spotfire, Tableau, ArcGis, etc.)
setwd("J:\\2016-07-02 ContourLinesPPP to dataframe")
#Create custom color ramp
ColorRamp= ColorRamp <- colorRampPalette(c("red","yellow",
"green","blue"), bias=1)
ColorRamp=ColorRamp(400)
#Force numbers to numeric
options(scipen=10)
#Read Text File
Example <- read.xlsx("Examplev2.xlsx",sheet=1, colNames = TRUE)
#Creat Spatial Points Dataframe
coordinates(Example) = ~Surface.Longitude+Surface.Latitude
Amountppp <- as(Example["Amount"], "ppp")
plot(Amountppp)
#Planar point pattern
as(as(Amountppp, "SpatialPoints"), "ppp")
## Planar point pattern: 654 points
## window: rectangle = [-103, -101] x [31, 33] units
str(Amountppp)
## List of 6
## $ window :List of 4
## ..$ type : chr "rectangle"
## ..$ xrange: num [1:2] -103 -101
## ..$ yrange: num [1:2] 31 33
## ..$ units :List of 3
## .. ..$ singular : chr "unit"
## .. ..$ plural : chr "units"
## .. ..$ multiplier: num 1
## .. ..- attr(*, "class")= chr "units"
## ..- attr(*, "class")= chr "owin"
## $ n : int 654
## $ x : num [1:654] -101 -103 -102 -103 -102 ...
## $ y : num [1:654] 31.9 31.2 31.6 31.5 31.6 ...
## $ markformat: chr "vector"
## $ marks : num [1:654] 101126 226170 159220 149334 156266 ...
## - attr(*, "class")= chr "ppp"
# Create an owin object from the county outline(s)
WestTx <- map('county', regions = c('texas,midland','texas,glasscock',
'texas,reagan', 'texas,upton') , fill=TRUE, col="transparent", plot=T, border='GREEN')
summary(WestTx)
## Length Class Mode
## x 29 -none- numeric
## y 29 -none- numeric
## range 4 -none- numeric
## names 4 -none- character
str(WestTx)
## List of 4
## $ x : num [1:29] -102 -102 -102 -101 -101 ...
## $ y : num [1:29] 32.1 32.1 32.1 32.1 31.7 ...
## $ range: num [1:4] -102.3 -101.3 31.1 32.1
## $ names: chr [1:4] "texas,glasscock" "texas,midland" "texas,reagan" "texas,upton"
## - attr(*, "class")= chr "map"
WestTxOwin<- map(WestTx, fill=TRUE,plot=T)
countypoly <- map2SpatialPolygons(WestTxOwin, IDs=WestTxOwin$names,
proj4string=CRS("+proj=longlat +datum=WGS84"))
countyowin <- as.owin.SpatialPolygons(countypoly)
Amountppp <- ppp(x=Amountppp$x, y=Amountppp$y, marks=Amountppp$marks, window=countyowin)
## Warning in ppp(x = Amountppp$x, y = Amountppp$y, marks = Amountppp$marks, :
## 472 points were rejected as lying outside the specified window
plot(Amountppp)
## Warning: Interpretation of arguments maxsize and markscale has changed (in
## spatstat version 1.37-0 and later). Size of a circle is now measured by its
## diameter.
## Warning in plot.ppp(Amountppp): 472 illegal points also plotted
#Remove points outside of window#
Amountppp=as.ppp(Amountppp)
plot(Amountppp)
#Create Map outline
map('county', regions = c('texas,midland','texas,glasscock',
'texas,reagan','texas,upton','texas,martin','texas,howard','texas, andrews',
'texas,ector', 'texas,crane','texas,irion','texas,sterling','texas,mitchell'),add=F)
# Plot a smoothed surface using the weights
plot(Smooth.ppp(Amountppp,sigma=.07),add=T,col = ColorRamp)
#Add points
plot(Amountppp, add=T, pch=16, axes=T, cex=.3)
#Contour Density
con=contour(Smooth(Amountppp,.07), axes=T,add=T)
#Add County Line
map('county', regions = c('texas,midland','texas,glasscock',
'texas,reagan','texas,upton','texas,martin','texas,howard','texas, andrews',
'texas,ector', 'texas,crane','texas,irion','texas,sterling','texas,mitchell'),add=T)
AmountSmooth= Smooth.ppp(Amountppp,.07, axes=F)
#Create Spatial Grid Dataframe
SGD= as(AmountSmooth,"SpatialGridDataFrame")
#Create image for gridded data
ImageSGD= as.image.SpatialGridDataFrame(SGD)
#Create contour lines from image
CtrImageSGD= contourLines(ImageSGD)
#Convert to Spatial lines dataframe
SLDF= ContourLines2SLDF(CtrImageSGD)
plot(SLDF)
#Convert to Raster
SmoothRaster= raster(AmountSmooth)
#Convert raster layer to contour lines
R2SLDF=rasterToContour(SmoothRaster)
#Write Shapefile
writeLinesShape(R2SLDF,"CONTOUR")