Computation of parcels area in Northern corridor for the AVANSE project

QUESTION

Is there a significative difference between computed and field area of parcels? Can we use quantitive characteristic of parcels to classify this exercice?
This research is conducted to answer the following indicators:

1. AVANSE Indicator 1.4 (FtF 4.5.2-2) Number of Hectares under Improved Technologies or Management Practices as a Result of USG Assistance (RiA) (WOG)
2. AVANSE Indicator 1.4.1 (F, FTF 4.5.1-28): Hectares under New or Improved/Rehabilitated Irrigation and Drainage Services as a Result of USG Assistance (RiA) (WOG)
3. AVANSE Indicator 1.5.1 (Custom) Number of Rural Hectares Mapped and Adjudicated (S)
4. AVANSE Indicator 2.1 (F 4.8.1-26) Number of Hectares of Biological Significance and/or Natural Resources under Improved Natural Resource Management as a Result of USG Assistance
5. AVANSE Indicator 2.3 (F 4.8.1-1) Number of Hectares of Biological Significance and/or Natural Resources Showing Improved Physical Conditions as a Result of USG Assistance

INPUT DATA

Data are collected in the field using GPS tracklogs and stored in .gpx format for every parcel. Each parcel is saved with a unique code that refers to the name of the surveyor and the crop. GPS receiver tracks are organized in the office by surveyor and device ID using a file storage system.

FEATURES

We read GPS tracks per file with the original coordinate system if they have coordinates. otherwise, the parcel without features is skipped. We transform set of points from each file into polyline, then into polygon by linking the first and end point. We add attribute data to spatial polygon by retrieving the name of the parcel from the metadata file, the name of file, the device ID, the name of surveyor from the GPS file storage system. We finally write each file into shapefile.

Show R code

# GETTING AND PROCESSING .GPX FILES
setwd("./parcels")
myfiles <- dir(pattern = "\\.gpx", recursive=TRUE)
for (i in 1:length(myfiles)) {   
      tryCatch({
            myfiles[i] %>%
                  readOGR(layer=ogrListLayers(myfiles[1])[3]) %>% 
                  SpatialLines2PolySet() %>% 
                  PolySet2SpatialPolygons() %>% 
                  SpatialPolygonsDataFrame(data=data.frame(surveyor=strsplit(
                        unlist(strsplit(myfiles[i],"_"))[3],"/")[[1]][1], device=unlist(strsplit(
                              myfiles[i],"_"))[2], filename=gsub("Piste_", "", strsplit(unlist(
                                    strsplit(myfiles[i],"/"))[2],"[.]")[[1]][1]), 
                        name=readOGR(myfiles[i], layer="tracks")@data[1,1])) %>% 
                  writeOGR("./shapefile", paste0('plot',i), driver="ESRI Shapefile")
            }, error=function(e){cat("ERROR :",conditionMessage(e), "in",myfiles[i],"\n")})
      }

We compile separate shapefiles into one shapefile. We create a numeric unique ID based on the name value from the metadata file.

Show R code

# COMBINING SHAPEFILES IN ONE SHAPEFILE
# retrieve shapefiles and append them in a list
shps <- dir("./shapefile", "*[0-9].shp")
shps <- gsub('.shp','',shps)
for (shp in shps) assign(shp,  readOGR('./shapefile',layer=shp))
lsts <- lapply(shps,get)
# to have unique ID we can write a function using the spChFIDs function from sp:
UniqueID <- function(SPDF){
      ID <- tolower(SPDF@data[1,4])
      newSPDF <- spChFIDs(SPDF,ID)    
      return(newSPDF)
}
# now we can insert new IDs for all of our parcels
newlst <- lapply(lsts, function(x) UniqueID(x))
parcels <- do.call(rbind, newlst)
writeOGR(parcels,"./shapefile", "plots", driver="ESRI Shapefile", overwrite_layer=TRUE)

In other to better visualize parcels in either static maps or web mapping application, we read parcels, sub-watersheds, road network spatial data and reproject them in the UTM zone 18 WGS 84 coordinate system to calculate the area of parcels in hectares. We also read the banana spreadsheet that will be joined to spatial parcels to compare computed and field parcels area.

Show R code

# reading necessary files and reprojecting
library(rgdal);library(plyr)
banana_data <- read.csv("C:/GIS/parcels/banana_parcel.csv")
plots <- readOGR(dsn = "C:/GIS/parcels/shapefile", "plots",verbose=FALSE)
plots_utm <- spTransform(plots, CRS('+proj=utm +zone=18 
                 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0'))
plots_utm$area_m2 <- sapply(plots_utm@polygons, function(x) x@Polygons[[1]]@area)
plots_utm$area_ha <- plots_utm$area_m2*10^-4
plots_utm@data <- join(plots_utm@data, banana_data, by="name", type="left", match="first")
subwatersheds <- readOGR(dsn = "C:/GIS/parcels/shapefile", "subwatesheds", verbose=FALSE)
subwatersheds_utm <- spTransform(subwatersheds, CRS('+proj=utm +zone=18 
                 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0'))
road <- readOGR(dsn = "C:/GIS/parcels/shapefile", "road", verbose=FALSE)
road_utm <- spTransform(road, CRS('+proj=utm +zone=18 
                 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0'))
plots_utm$subwatershed <- over(plots_utm, subwatersheds_utm)$Name

We also calculate the size and the length of raw files and processed files. Knowing the size of the raw file is a good indicator to set how much memory we will allocate to our Integrated Development Environment(IDE), in our case Rstudio.

Summary of processed parcels GPS tracks

Show R code

library(knitr)
setwd("C:/GIS/parcels")
myfiles <- dir(pattern = "\\.gpx", recursive=TRUE)
size_file <- round(sum(file.info(myfiles)$size)/2^20,2)
size_plot <- round(file.info("./shapefile/plots.shp")$size/2^20,2)
length_files <- length(myfiles)
nrow_plots <- nrow(plots)
ha_GPS <- round(sum(banana_data$field_ha[!is.na(banana_data$field_ha)]),2)
ha_plot <- round(sum(plots_utm$area_ha),2)
size <- c(size_file, size_plot)
length <- c(length_files,nrow_plots)
ha <- c(ha_GPS, ha_plot)
DT<-data.frame(cbind(size,length,ha),row.names=c("GPS files","Shapefiles"))
colnames(DT) <- c('Size in MB','Number of parcels',"Number of hectares")
kable(DT, align='c')

	Size in MB	Number of parcels	Number of hectares
GPS files	7.44	1645	727.62
Shapefiles	0.90	1446	646.94

The number of files decrease after processing from 1645 to 1446. It is mainly due to missing features from GPS tracks. The number of size decrease from 7.44 to 0.9. It is one of the objective of spatial data processing. GPS files is converted in Shapefile because it is a handy format when dealing with spatial data.

Exploratory analysis

The goal of this task is to understand the basic relationships observed between surveyors in the dataset and prepare to build a model. We can organize surveyor by their frequency, which refers to the frequency of the parcels collected.

Show R code

library(ggplot2)
# linear relationships between computed and field area
sort(summary(plots_utm$surveyor),decreasing=TRUE)

##     Ronald      Velia   Rochelin PierreJoel      Nixon      Gonel 
##        161        148        139        127        101        100 
##     Jacquy       Bony       Rony    Suzette    William      Sonel 
##         84         68         68         49         43         41 
##    Henriot   Geffrard    Tardier    Anthony     Sandro   Dorceant 
##         40         38         36         27         27         24 
##     Folero    Otilien    Semeran   Nicodeme      Billy       Remy 
##         24         24         23         20         18         12 
##      Lesly 
##          4

We can also aggregate computed area by surveyor. The basic idea is to explore in a bar plot the area covered by each surveyor.

Show R code

aggr_plots <- aggregate(area_ha ~ surveyor, data=plots_utm@data, sum)
qplot(surveyor,area_ha,data=aggr_plots, ymax=area_ha,ymin=0,label=" ", geom=c("pointrange", "text"), vjust=-0.5,color=I("blue")) + ggtitle("Number of computed area covered per surveyors") + geom_text(aes(label=round(area_ha,1),size=10)) + coord_flip() + theme(legend.position ="none") + theme(axis.text.y=element_text(size = 10, hjust=1,color="black"))

We can also make a representation on the cloud of the surveyors showing their relative size based on the number of parcels collected or the total computed parcels area. The two graphs shows respectively the frequency of parcels per susrveyor and the number of parcels area computed per surveyor. We can see that Ronald and Velia are the top surveyors.

Show R code

# Plot the 40 most frequently occurring words in the corpus
library(wordcloud)
par(mfrow=c(1,2))
set.seed(142)   
dark2 <- brewer.pal(8, "Dark2")
aggr_name <- table(plots_utm$surveyor)
wordcloud(names(aggr_name), as.numeric(aggr_name), max.words=25, rot.per=0.2, colors=dark2) 
wordcloud(aggr_plots$surveyor, aggr_plots$area_ha, max.words=20, rot.per=0.2, colors=dark2) 

We can see that Ronald has collected more parcels than his buddies but his efficiency in area covered is less than Gonel and Velia

Can we compare the efficiency of surveyors.

First of all let’s explore the relationship between the two source of parcels area disregrading surveyors.

Show R code

library(ggplot2)
# linear relationships between computed and field area
qplot(field_ha,area_ha,data=plots_utm@data, geom=c("point","smooth"), method="lm", main="Computed area versus field area",xlab="GPS lecture", ylab="Computation")

There’s a positive relationships and it seems to be strong for all surveyors

Show R code

qplot(field_ha,area_ha,data=plots_utm@data[plots_utm$surveyor==c("Ronald","Velia","Rochelin", "PierreJoel"),],color=surveyor, geom=c("point","smooth"), method="lm", facets=.~surveyor, main="Computation area versus field area",xlab="GPS lecture",ylab="Computation")

If we split this relationship by the top four(4) frequent surveyors to have a closer look at the relationship, we see it is positive and strong for all of them, but much more for Velia.

Regression models

Let’s build a model with computed area as outcome and only surveyor as predictor. The basic idea in this exercice is that we want to understand the difference in computed area between surveyors. The efficiency of a surveyor can be mostly explained by the total computed area, because number of hecatares are key terms in our five (5) indicators above mentioned.

Show R code

fit1<-lm(area_ha ~ surveyor, data=plots_utm@data)
summary(fit1)$coef

##                       Estimate Std. Error    t value     Pr(>|t|)
## (Intercept)         0.38709625 0.08627239  4.4869079 7.809141e-06
## surveyorBilly       0.26395647 0.13640862  1.9350425 5.318331e-02
## surveyorBony        0.21699358 0.10197158  2.1279809 3.351038e-02
## surveyorDorceant   -0.18115750 0.12576253 -1.4404728 1.499540e-01
## surveyorFolero     -0.05595666 0.12576253 -0.4449390 6.564316e-01
## surveyorGeffrard    0.06762164 0.11283314  0.5993065 5.490640e-01
## surveyorGonel       0.32111415 0.09722404  3.3028266 9.809488e-04
## surveyorHenriot    -0.02168810 0.11165527 -0.1942416 8.460145e-01
## surveyorJacquy      0.27378341 0.09917305  2.7606635 5.842578e-03
## surveyorLesly      -0.11339899 0.24017216 -0.4721571 6.368872e-01
## surveyorNicodeme   -0.02845405 0.13225306 -0.2151485 8.296824e-01
## surveyorNixon      -0.05461095 0.09712166 -0.5622942 5.740043e-01
## surveyorOtilien     0.12663408 0.12576253  1.0069301 3.141398e-01
## surveyorPierreJoel -0.06339348 0.09500147 -0.6672895 5.046956e-01
## surveyorRemy       -0.09680108 0.15552975 -0.6223959 5.337814e-01
## surveyorRochelin   -0.01147140 0.09427976 -0.1216741 9.031744e-01
## surveyorRonald     -0.02643322 0.09322616 -0.2835387 7.768052e-01
## surveyorRony        0.20450040 0.10197158  2.0054646 4.510271e-02
## surveyorSandro      0.31438707 0.12200758  2.5767831 1.007257e-02
## surveyorSemeran    -0.20711291 0.12720169 -1.6282245 1.036990e-01
## surveyorSonel       0.07237510 0.11110519  0.6514106 5.148868e-01
## surveyorSuzette    -0.03882504 0.10744360 -0.3613527 7.178895e-01
## surveyorTardier    -0.06247662 0.11412764 -0.5474276 5.841709e-01
## surveyorVelia       0.06050020 0.09381234  0.6449067 5.190918e-01
## surveyorWilliam     0.22634095 0.11007447  2.0562530 3.994043e-02

The estimate of billy is 0.26, but the P-value is 5.31e-02 > 0.05. The test statistic is significative. in other words, the t-test for \(H_0: \beta_{surveyor} = 0\) versus \(H_a: \beta_{surveyor} \neq 0\) with P-value equal to 0.0531 is significant. Notice that we have, Anthony who is missing in the table. It is because he is the reference category. The number 0.26 is the estimated increase in computed area comparing billy parcels and Anthony parcels. We reject the null hypothesis and claim that there is a significant difference in the mean computed area between anthony parcels and that of billy parcels.

Show R code

fit2<-lm(area_ha ~ surveyor-1, data=plots_utm@data)
summary(fit2)$coef

##                     Estimate Std. Error   t value     Pr(>|t|)
## surveyorAnthony    0.3870963 0.08627239  4.486908 7.809141e-06
## surveyorBilly      0.6510527 0.10566166  6.161674 9.355895e-10
## surveyorBony       0.6040898 0.05436248 11.112257 1.440576e-27
## surveyorDorceant   0.2059387 0.09150568  2.250557 2.456600e-02
## surveyorFolero     0.3311396 0.09150568  3.618787 3.063666e-04
## surveyorGeffrard   0.4547179 0.07272134  6.252881 5.318441e-10
## surveyorGonel      0.7082104 0.04482845 15.798236 6.408757e-52
## surveyorHenriot    0.3654082 0.07088000  5.155307 2.889094e-07
## surveyorJacquy     0.6608797 0.04891185 13.511649 3.156538e-39
## surveyorLesly      0.2736973 0.22414223  1.221087 2.222555e-01
## surveyorNicodeme   0.3586422 0.10023945  3.577855 3.580614e-04
## surveyorNixon      0.3324853 0.04460597  7.453829 1.568936e-13
## surveyorOtilien    0.5137303 0.09150568  5.614190 2.371586e-08
## surveyorPierreJoel 0.3237028 0.03977881  8.137567 8.718479e-16
## surveyorRemy       0.2902952 0.12940858  2.243245 2.503457e-02
## surveyorRochelin   0.3756248 0.03802299  9.878887 2.638428e-22
## surveyorRonald     0.3606630 0.03532977 10.208475 1.164087e-23
## surveyorRony       0.5915967 0.05436248 10.882445 1.507519e-26
## surveyorSandro     0.7014833 0.08627239  8.131030 9.179785e-16
## surveyorSemeran    0.1799833 0.09347377  1.925496 5.436674e-02
## surveyorSonel      0.4594714 0.07001027  6.562913 7.374373e-11
## surveyorSuzette    0.3482712 0.06404064  5.438285 6.325297e-08
## surveyorTardier    0.3246196 0.07471408  4.344825 1.492872e-05
## surveyorVelia      0.4475965 0.03684876 12.146854 2.255475e-32
## surveyorWilliam    0.6134372 0.06836274  8.973268 8.881315e-19

If we omit the intercept, then the model includes now all the surveyors. Now Billy is not a linear combination of Anthony. there are 25 means in the dataset, the total number of surveyors. The expected value of the outcome should be the mean for one or another surveyor. As we can see in the table billy is about 0.65 and anthony is about 0.38 and it is clearly illustrated in the boxplot below.

Show R code

plots_utm$surveyor <- as.character(plots_utm$surveyor) 
boxplot(area_ha ~ surveyor, data = plots_utm@data, subset=surveyor==c("Anthony","Billy"),
        xlab = "Name of surveyor", ylab = "Computed area in ha",
        main = "Comparison of surveyors efficiency", col = "lightgray")

If I were to substract 0.38 to 0.65, we would get exactly 0.27. In other words, the model is perfectly consistent. If you are going to fix the mean of the billy and anthony, any way you do this, the model will alway be consistent.

Show R code

alpha <- 0.05
n<-nrow(mtcars)
pe <- coef(summary(fit1))["surveyorBilly", "Estimate"]
se <- coef(summary(fit1))["surveyorBilly", "Std. Error"]
tstat <- qt(1 - alpha/2, n - 2)  # n - 2 for model with intercept and slope
pe + c(-1, 1) * (se * tstat)

## [1] -0.01462709  0.54254004

If we were willing to choose the model 1 as our best model, then the confidence interval for the 0.263 computed area difference between Billy and Anthony would be -0.01462709 and 0.54254004. The graph below shows residuals versus fitted computed area value.

Show R code

rownames(plots_utm@data) <-plots_utm$name
plot(lm(area_ha ~ surveyor, data=plots_utm@data), which=c(1))

Conclusion 1

We have checked the difference between computed parcels area and written parcels area from surveyors. The difference is not signicant. We correctly accept the null hypothesis and claim that there is not a significant difference in computed parcels area versus written parcels area from surveyors.

STATIC MAPS

At this step, we visualize the difference in computed area and field area in a spatial context. Field area are read from the GPS device by the surveyor and transcripted manually in a paper survey form. They are then entered in a spreadsheet table by data entry clerks. Given How parcels are scattered over the watersheds, they are subsetted at the level of subwatershed for each map to have a closer look when comparing the two sources of area.

Show R code

# plotting parcels area in Pont Parois
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(798000,2178500), scale = 1000, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(798000,2178400), "0",which=2, cex = .7)
text2 = list("sp.text", c(799000,2178400), "1000 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(799000,2178500), scale = 700,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_parois <- subset(plots_utm, subwatershed == "Pont Parois")
# plotting with scale bar, north arrow, and Title
spplot(plots_parois, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Pont Parois sub-watershed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.90, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Tabary
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(800000,2176500), scale = 1000, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(800000,2176400), "0",which=2, cex = .7)
text2 = list("sp.text", c(801000,2176400), "1000 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(800000,2176000), scale = 700,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_tabary <- subset(plots_utm, subwatershed == "Tabary")
# plotting with scale bar, north arrow, and Title
spplot(plots_tabary, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Tabary sub-watershed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.90, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Bois-de-Lance
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(802000,2171000), scale = 500, fill=c("transparent","black"), which=2)
text1 = list("sp.text", c(802000,2170900), "0",which=2, cex = .7)
text2 = list("sp.text", c(802500,2170900), "500 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(802200,2171000), scale = 500, which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_lance <- subset(plots_utm, subwatershed == "Bois de Lance")
# plotting with scale bar, north arrow, and Title
spplot(plots_lance, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Bois de Lance sub-watershed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.90, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Fourier
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(799500,2172000), scale = 1000, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(799500,2171900), "0",which=2, cex = .7)
text2 = list("sp.text", c(800500,2171900), "1000 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(800500,2172000), scale = 600,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_fourier <- subset(plots_utm, subwatershed == "Fourier")
# plotting with scale bar, north arrow, and Title
spplot(plots_fourier, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Fourier sub-watershed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.85, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Sanson
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(791000,2172000), scale = 1000, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(791000,2171900), "0",which=2, cex = .7)
text2 = list("sp.text", c(792000,2171900), "1000 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(794000,2174800), scale = 600,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_sanson <- subset(plots_utm, subwatershed == "Sanson")
# plotting with scale bar, north arrow, and Title
spplot(plots_sanson, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Sanson sub-wateshed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.90, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Lamatry
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(837000,2166500), scale = 2000, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(837000,2166300), "0",which=2, cex = .7)
text2 = list("sp.text", c(839000,2166300), "2000 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(835000,2170000), scale = 700,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_lamatry <- subset(plots_utm, subwatershed == "Lamatry")
# plotting with scale bar, north arrow, and Title
spplot(plots_lamatry, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Lamatry sub-wateshed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(.915, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Show R code

# plotting parcels area in Colette
library(sp);library(lattice);library(grid)
trellis.par.set(add.line = list(col = "transparent"))
scale = list("SpatialPolygonsRescale", layout.scale.bar(), 
             offset = c(831000,2169500), scale = 500, fill=c("transparent","black"),which=2)
text1 = list("sp.text", c(831000,2169400), "0",which=2, cex = .7)
text2 = list("sp.text", c(831500,2169400), "500 m",which=2, cex = .7)
arrow = list("SpatialPolygonsRescale", layout.north.arrow(type=1), 
             offset = c(831500,2169500), scale = 600,which=1)
rd = list("sp.lines",road_utm, col='gray')
plots_colette <- subset(plots_utm, subwatershed == "Colette")
# plotting with scale bar, north arrow, and Title
spplot(plots_colette, c("area_ha", "field_ha"), names.attr = c("computed","written"),
       colorkey=list(space="right", align="bottom"), scales = list(draw = TRUE), main = 
       "Computed versus field parcels area in Colette sub-wateshed", col.regions = 
       colorRampPalette(c("yellow","orange","purple","blue","dark blue"))(16), 
       ylab = "Northing (m)", xlab = "Easting (m)",
       sp.layout = list(scale,text1,text2,arrow,rd), as.table = TRUE)
grid.text("area (ha)", x=unit(0.86, "npc"), y=unit(0.50, "npc"), rot=-90, check.overlap = TRUE)

Conclusion 2

We have compared with a panel visualization parcels area from both computation and GPS lecture method. we have aggregated this comparison at the subwatershed level so as to have legible information. The difference between the two methods is not significant. We correctly accept the null hypothesis and claim that there is not a significant difference in computed parcels area versus GPS lecture parcels area.