Library
library(terra) # loads the terra library
## terra 1.7.78
library(rpart) # loasd the recursive partitioning and regression tree library
Upload TIF File
rawImage <- rast("C:/Users/154455/Dropbox/2.dAILYr/R terra/LC08_016028_20210824.TIF")
Data Structure
rawImage
## class : SpatRaster
## dimensions : 2169, 3545, 19 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 369645, 475995, 4982775, 5047845 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source : LC08_016028_20210824.TIF
## names : SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, ...
Plot TIF File
plot(rawImage)
Plot grey TIF File
plot(rawImage$SR_B5, col=gray(0:255 / 255))
Plot RGB
plotRGB(c(rawImage$SR_B4, rawImage$SR_B3, rawImage$SR_B2), stretch="lin")
Calculation
NDVI <- ((rawImage$SR_B5 - rawImage$SR_B4)/ (rawImage$SR_B5 + rawImage$SR_B4))
NDVI
## class : SpatRaster
## dimensions : 2169, 3545, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 369645, 475995, 4982775, 5047845 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## varname : LC08_016028_20210824
## name : SR_B5
## min value : -0.1436602
## max value : 0.9495935
Histogram
hist(NDVI, xlab="NDVI", ylab="Frequency")
## Warning: [hist] a sample of 13% of the cells was used (of which 20% was NA)
Plot NDVI
plot(NDVI, col=rev(terrain.colors(10)), main="NDVI")
Raster
rast_nbrhds <- rast(NDVI)
rast_nbrhds
## class : SpatRaster
## dimensions : 2169, 3545, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 369645, 475995, 4982775, 5047845 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
CRS
crs(NDVI)
## [1] "PROJCRS[\"WGS 84 / UTM zone 18N\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 18N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-75,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n AREA[\"Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamaica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.\"],\n BBOX[0,-78,84,-72]],\n ID[\"EPSG\",32618]]"
Vector
nbrhds <- vect("C:/Users/154455/Dropbox/2.dAILYr/R terra/NeighbourhoodNDVI.gpkg")
NBRHDS
nbrhds_proj <- terra::project(nbrhds, "EPSG:32618") # reproject Neighbourhoods to UTM 18N
nbrhdIDrast <- rasterize(nbrhds_proj, rast_nbrhds, "ONS_ID")
Plot NBRHDS
plot(nbrhdIDrast)
lines(nbrhds_proj)
Zonal
NDVIbyNbrhd <- zonal(NDVI, nbrhdIDrast, "mean", na.rm=TRUE)
zonal(NDVI, nbrhdIDrast, "mean", na.rm=TRUE)
## ONS_ID SR_B5
## 1 3001 0.2718333
## 2 3002 0.2549597
## 3 3003 0.1818396
## 4 3004 0.2005813
## 5 3005 0.2509934
## 6 3006 0.3060857
## 7 3007 0.2427334
## 8 3008 0.1789586
## 9 3009 0.2088654
## 10 3010 0.2776771
## 11 3011 0.2790439
## 12 3012 0.2314434
## 13 3013 0.2249226
## 14 3014 0.2197906
## 15 3015 0.2414667
## 16 3016 0.2422117
## 17 3017 0.3026289
## 18 3018 0.2158740
## 19 3019 0.1583075
## 20 3020 0.1855722
## 21 3021 0.3573304
## 22 3022 0.2200641
## 23 3023 0.2223734
## 24 3024 0.1176931
## 25 3025 0.2633298
## 26 3026 0.2779497
## 27 3027 0.1877151
## 28 3028 0.1984221
## 29 3029 0.1898260
## 30 3030 0.1498353
## 31 3031 0.1942585
## 32 3032 0.2432974
## 33 3033 0.3161335
## 34 3034 0.2319828
## 35 3035 0.2305140
## 36 3036 0.2744971
## 37 3037 0.3933624
## 38 3038 0.3104095
## 39 3039 0.3747193
## 40 3040 0.2396934
## 41 3041 0.2756470
## 42 3042 0.3246351
## 43 3043 0.2117383
## 44 3044 0.2881309
## 45 3045 0.2030530
## 46 3046 0.3420093
## 47 3047 0.1776597
## 48 3048 0.2187456
## 49 3049 0.3208577
## 50 3050 0.3539877
## 51 3051 0.3185859
## 52 3052 0.2385591
## 53 3053 0.2198501
## 54 3054 0.2020800
## 55 3055 0.1505994
## 56 3056 0.2066737
## 57 3057 0.1994933
## 58 3058 0.2072209
## 59 3059 0.2034846
## 60 3060 0.2510204
## 61 3061 0.2229429
## 62 3062 0.3766720
## 63 3063 0.1987518
## 64 3064 0.1552076
## 65 3065 0.2096238
## 66 3066 0.2268546
## 67 3067 0.1461444
## 68 3068 0.1250520
## 69 3069 0.2382274
## 70 3070 0.3486916
## 71 3071 0.3640045
## 72 3072 0.3247830
## 73 3073 0.3913914
## 74 3074 0.3748339
## 75 3075 0.3263166
## 76 3076 0.4126848
## 77 3077 0.2253755
## 78 3078 0.3917863
## 79 3079 0.1727064
## 80 3080 0.2459583
## 81 3081 0.2741774
## 82 3082 0.1721655
## 83 3083 0.2062937
## 84 3084 0.1282868
## 85 3085 0.2276173
## 86 3086 0.3894057
## 87 3087 0.1762890
## 88 3088 0.1913196
## 89 3089 0.2119673
## 90 3090 0.2432336
## 91 3091 0.2947103
## 92 3092 0.2309628
## 93 3093 0.2190038
## 94 3094 0.2669745
## 95 3095 0.2179608
## 96 3097 0.2084639
## 97 3098 0.2339770
## 98 3099 0.2794043
## 99 3100 0.3032717
## 100 3101 0.1792289
## 101 3102 0.2833546
## 102 3103 0.2429410
## 103 3104 0.1529106
## 104 3105 0.1674837
## 105 3106 0.2355038
## 106 3107 0.2232746
## 107 3108 0.2241174
## 108 3109 0.2225854
## 109 3110 0.2465480
## 110 3111 0.1778035
## 111 3112 0.1507260
## 112 3113 0.4059962
## 113 3114 0.2657269
## 114 3115 0.1236391
## 115 3116 0.2027260
## 116 3117 0.2162043
Data
names(NDVIbyNbrhd)[2] <- "NDVI"
names(NDVI) <- "NDVI"
newNbrhds <- merge(nbrhds_proj, NDVIbyNbrhd)
newNbrhds[newNbrhds$NDVI == max(newNbrhds$NDVI)]$ONS_Name
## [1] "NAVAN - SARSFIELD"
Upload vector
lcpolys <- vect("D:/ICT Backup_154455/Downloads/SpatialR_SMexamples.gpkg")
Plot RGB
plotRGB(c(rawImage$SR_B4, rawImage$SR_B3, rawImage$SR_B2), stretch="lin")
polys(lcpolys)
text(lcpolys, lcpolys$LandCover)
Subset
LS8Image <- subset(rawImage, 1:9) # use just the first 9 bands in the source file
LS8Image
## class : SpatRaster
## dimensions : 2169, 3545, 9 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 369645, 475995, 4982775, 5047845 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source : LC08_016028_20210824.TIF
## names : SR_B1, SR_B2, SR_B3, SR_B4, SR_B5, SR_B6, ...
Sample data plot
set.seed(1) # create a seed for random number generation
trainingSamples<-spatSample(lcpolys, 600, method="random") # Sample 600 points randomly from the polygon data
plot(trainingSamples)
table(trainingSamples$LandCover)
##
## Builtup Crops Parkland Residential Water Wetlands
## 40 133 208 28 128 63
sampledImage <- extract(LS8Image, trainingSamples, ID=FALSE)
head(sampledImage)
## SR_B1 SR_B2 SR_B3 SR_B4 SR_B5 SR_B6 SR_B7 SR_QA_AEROSOL ST_B10
## 1 8738 9005 9821 9783 15483 13011 10980 96 48595
## 2 8245 8227 8899 8214 19546 13068 9574 96 45097
## 3 7294 7432 7687 7567 7600 7552 7516 164 44357
## 4 8426 8564 9427 8815 23238 14026 10162 96 46170
## 5 8366 8447 9122 8754 20449 14078 10395 96 46425
## 6 7328 7469 7766 7656 7647 7591 7563 164 44290
sampData <- data.frame(class = trainingSamples$LandCover, sampledImage)
head(sampData)
## class SR_B1 SR_B2 SR_B3 SR_B4 SR_B5 SR_B6 SR_B7 SR_QA_AEROSOL ST_B10
## 1 Residential 8738 9005 9821 9783 15483 13011 10980 96 48595
## 2 Parkland 8245 8227 8899 8214 19546 13068 9574 96 45097
## 3 Water 7294 7432 7687 7567 7600 7552 7516 164 44357
## 4 Crops 8426 8564 9427 8815 23238 14026 10162 96 46170
## 5 Crops 8366 8447 9122 8754 20449 14078 10395 96 46425
## 6 Water 7328 7469 7766 7656 7647 7591 7563 164 44290
rpart
cartmodel <- rpart(as.factor(class)~., data=sampData, method='class', minsplit=5)
print(cartmodel)
## n= 600
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 600 392 Parkland (0.067 0.22 0.35 0.047 0.21 0.1)
## 2) SR_B3>=8281 472 264 Parkland (0.085 0.28 0.44 0.059 0 0.13)
## 4) ST_B10< 45285 202 6 Parkland (0 0.03 0.97 0 0 0)
## 8) SR_B5>=26179.5 4 0 Crops (0 1 0 0 0 0) *
## 9) SR_B5< 26179.5 198 2 Parkland (0 0.01 0.99 0 0 0) *
## 5) ST_B10>=45285 270 143 Crops (0.15 0.47 0.044 0.1 0 0.23)
## 10) SR_B5>=18695.5 119 6 Crops (0.017 0.95 0.017 0.017 0 0) *
## 11) SR_B5< 18695.5 151 88 Wetlands (0.25 0.093 0.066 0.17 0 0.42)
## 22) SR_B2>=8319.5 74 36 Builtup (0.51 0.12 0.041 0.32 0 0)
## 44) ST_B10>=48752 37 0 Builtup (1 0 0 0 0 0) *
## 45) ST_B10< 48752 37 13 Residential (0.027 0.24 0.081 0.65 0 0)
## 90) SR_B6>=14142 10 1 Crops (0 0.9 0 0.1 0 0) *
## 91) SR_B6< 14142 27 4 Residential (0.037 0 0.11 0.85 0 0) *
## 23) SR_B2< 8319.5 77 14 Wetlands (0 0.065 0.091 0.026 0 0.82)
## 46) ST_B10< 46442.5 16 9 Parkland (0 0.31 0.44 0.063 0 0.19)
## 92) ST_B10< 45602.5 7 0 Parkland (0 0 1 0 0 0) *
## 93) ST_B10>=45602.5 9 4 Crops (0 0.56 0 0.11 0 0.33) *
## 47) ST_B10>=46442.5 61 1 Wetlands (0 0 0 0.016 0 0.98) *
## 3) SR_B3< 8281 128 0 Water (0 0 0 0 1 0) *
plot(cartmodel, uniform=TRUE, main="Classification Tree")
text(cartmodel, cex = 1)
Plot Classify
classify <- predict(LS8Image, cartmodel, na.rm=TRUE) # will take some time
## |---------|---------|---------|---------|=========================================
plot(classify)
Landcover
landcover <- which.max(classify)
landcover
## class : SpatRaster
## dimensions : 2169, 3545, 1 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 369645, 475995, 4982775, 5047845 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : which.max
## min value : 1
## max value : 6
Plot
myClasses <- c("Built-up", "Crops", "Parkland", "Residential", "Water", "Wetlands")
levels(landcover) <- data.frame(id = 1:6, class=myClasses)
mycolors <- c("red", "lightgreen","darkgreen","yellow","blue","brown")
plot(landcover, col=mycolors)
For value
k <- 10 # number of times to cross-validate (k-fold cross-validation)
j <- sample(rep(1:k, each = round(nrow(sampData))/k))
resultList <- list() # create a new empty list to hold each validation result
for (x in 1:k) {
train <- sampData[j != x, ]
validation <- sampData[j == x, ]
# Create a CART model as above but with new training subset
cartmodel <- rpart(as.factor(class)~., train, method = 'class', minsplit = 5)
pclass <- predict(cartmodel, validation, na.rm = TRUE)
# assign class to maximum probablity
pc <- apply(pclass, 1, which.max)
# use labels instead of numbers
pc <- colnames(pclass)[pc]
# create a data.frame using the reference and prediction
resultList[[x]] <- cbind(validation$class, pc)
}
# create a data.frame of results based on the lists we just created
resultDF <- do.call(rbind, resultList)
resultDF <- data.frame(resultDF)
colnames(resultDF) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(resultDF)
print(conmat)
## predicted
## observed Builtup Crops Parkland Residential Water Wetlands
## Builtup 35 3 0 1 0 1
## Crops 0 121 7 1 0 4
## Parkland 0 4 202 2 0 0
## Residential 1 3 0 21 0 3
## Water 0 0 0 0 128 0
## Wetlands 0 4 0 4 0 55
Calculate overall accuracy
# Calculate overall accuracy
totalSamples <- sum(conmat)
diagonal <- diag(conmat)
overallAcc <- sum(diagonal) / totalSamples
Summary
summary(zonal(NDVI, classify, "mean", na.rm=TRUE))
## Builtup Crops Parkland Residential
## Min. :0.0000 Min. :0.000000 Min. :0.00000 Min. :0.000000
## 1st Qu.:0.0000 1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.000000
## Median :0.0000 Median :0.005051 Median :0.00000 Median :0.008197
## Mean :0.1054 Mean :0.341524 Mean :0.21178 Mean :0.109616
## 3rd Qu.:0.0126 3rd Qu.:0.813889 3rd Qu.:0.08754 3rd Qu.:0.079202
## Max. :1.0000 Max. :1.000000 Max. :1.00000 Max. :0.851852
## Water Wetlands NDVI
## Min. :0.0 Min. :0.0000 Min. :0.02458
## 1st Qu.:0.0 1st Qu.:0.0000 1st Qu.:0.25814
## Median :0.0 Median :0.0000 Median :0.32840
## Mean :0.1 Mean :0.1317 Mean :0.30390
## 3rd Qu.:0.0 3rd Qu.:0.0000 3rd Qu.:0.37446
## Max. :1.0 Max. :0.9836 Max. :0.52672