Calculating NDVI index in Warsaw
The goal is to calculate NDVI index for Warsaw city for each cell separately and save each observation into a data frame, along with the cell’s coordinates.
Step 1 - downloading Sentinel-2 data
Log in on the Scihub copernicus website and mark your area of interest (in my case, Warsaw city):
Use Advanced Search to filter observations with 0% of cloud cover:
Click “search” and download observations you’re interested in (in my case, one observation for lower-left tile, one for upper left tile and one for upper-right tile. Put together, they cover the whole Warsaw’s area)
Open downloaded folder, and then GRANULE->(really long and senseless name of the folder)->IMG_DATA. There you will find some files with extensions .jp2. Copy those files to your working directory.
Step 2 - reading data into R and extracting the image of Warsaw from a tile (in my case, UDC tile was the first one)
library(devtools)
library(sen2r)
library(raster)
library(RStoolbox)
library(tidyr)
library(dplyr)
library(rlang)
library(ggplot2)
library(caret)
library(spdep)
library(rgdal)
library(maptools)
library(sp)
library(classInt)
library(GISTools)
library(maps)
library(rgeos)
library(readr)
library(pryr)
library(tabularaster)
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat/Bands0")
# Reading Sentinel-2 data into a raster object
s=stack("T34UDC_20211009T095029_B02.jp2"," T34UDC_20211009T095029_B03.jp2","T34UDC_20211009T095029_B04.jp2"," T34UDC_20211009T095029_B08.jp2")
# Reading & extracting Warsaw area
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat")
pow.sf <- st_read("powiaty.shp")## Reading layer `powiaty' from data source
## `C:\Users\ewado\OneDrive\Dokumenty\Studia\R-mat\powiaty.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 380 features and 29 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 171677.6 ymin: 133223.7 xmax: 861895.7 ymax: 774923.7
## Projected CRS: ETRS89 / Poland CS92pow.sf <- st_transform(pow.sf, crs = "+proj=longlat +datum=NAD83")
waw.pow.sf<-pow.sf[pow.sf$jpt_nazwa_=='powiat Warszawa',]
waw.pow.sf2<-st_transform(waw.pow.sf, crs(s)@projargs)
waw.sat.s<-crop(s, extent(waw.pow.sf2))
waw.sat_przyciety<-mask(waw.sat.s, mask=waw.pow.sf2)
was<-waw.sat_przyciety
ndvi<-(was[[4]]-was[[3]])/(was[[4]]+was[[3]])
plot(waw.sat_przyciety[[1]]) 
plot(ndvi) 
Step 3 - creating a data frame
Create a spatial object using xyFromCell function (the function that returns coordinates of the center of the raster cell). I used xyFromCell(name_of_the_raster, spatial=TRUE, all_cells) and change it’s projection. Then, save coordinates of this spatial object using a function coordinates():
ndvi.wsp<-spTransform(xyFromCell(ndvi, spatial=TRUE, 1:5591900), CRS("+proj=longlat +datum=NAD83"))
crds<-coordinates(ndvi.wsp)
class(crds)## [1] "matrix" "array"dim(crds)## [1] 5591900 2head(crds)## x y
## [1,] 20.85120 52.35033
## [2,] 20.85135 52.35033
## [3,] 20.85150 52.35033
## [4,] 20.85164 52.35034
## [5,] 20.85179 52.35034
## [6,] 20.85194 52.35034It looks nice, but it’s a matrix, and we want a data frame:
crds.df<-data.frame(ID=c(1:5591900), X=crds[,1], Y=crds[,2])
class(crds.df)## [1] "data.frame"dim(crds.df)## [1] 5591900 3head(crds.df)## ID X Y
## 1 1 20.85120 52.35033
## 2 2 20.85135 52.35033
## 3 3 20.85150 52.35033
## 4 4 20.85164 52.35034
## 5 5 20.85179 52.35034
## 6 6 20.85194 52.35034The data frame of the coordinates is ready, now we need to add corresponding values of NDVI index:
ndvi.wart<-as.data.frame(ndvi)
crds.df$NDVI<-ndvi.wart[,1]
head(crds.df)## ID X Y NDVI
## 1 1 20.85120 52.35033 NA
## 2 2 20.85135 52.35033 NA
## 3 3 20.85150 52.35033 NA
## 4 4 20.85164 52.35034 NA
## 5 5 20.85179 52.35034 NA
## 6 6 20.85194 52.35034 NAsummary(crds.df)## ID X Y NDVI
## Min. : 1 Min. :20.85 Min. :52.10 Min. :-1.0
## 1st Qu.:1397976 1st Qu.:20.92 1st Qu.:52.16 1st Qu.: 0.2
## Median :2795950 Median :21.00 Median :52.22 Median : 0.4
## Mean :2795950 Mean :21.00 Mean :52.22 Mean : 0.4
## 3rd Qu.:4193925 3rd Qu.:21.07 3rd Qu.:52.29 3rd Qu.: 0.6
## Max. :5591900 Max. :21.14 Max. :52.35 Max. : 1.0
## NA's :1567343As you can wee, it worked :)
crds.df[100000:100010,]## ID X Y NDVI
## 100000 100000 20.92447 52.34591 0.6436494
## 100001 100001 20.92462 52.34591 0.6038692
## 100002 100002 20.92477 52.34591 0.5959860
## 100003 100003 20.92491 52.34591 0.5638566
## 100004 100004 20.92506 52.34591 0.5771626
## 100005 100005 20.92521 52.34591 0.6607187
## 100006 100006 20.92535 52.34591 0.6931373
## 100007 100007 20.92550 52.34591 0.6766121
## 100008 100008 20.92565 52.34591 0.6623377
## 100009 100009 20.92579 52.34591 0.6694659
## 100010 100010 20.92594 52.34591 0.6800000library(ggplot2)
ggplot(crds.df, aes(x=X, y=Y, col=NDVI))+
geom_point()
You can also remove observations with NA as ndvi index:
dim(crds.df)## [1] 5591900 4crds.df<-crds.df[!is.na(crds.df$NDVI),]
dim(crds.df)## [1] 4024557 4head(crds.df)## ID X Y NDVI
## 439 439 20.91551 52.3504 -0.1923536
## 440 440 20.91566 52.3504 -0.2149068
## 441 441 20.91580 52.3504 -0.2272175
## 442 442 20.91595 52.3504 -0.2250923
## 443 443 20.91610 52.3504 -0.2322738
## 444 444 20.91624 52.3504 -0.2307692ggplot(crds.df, aes(x=X, y=Y, col=NDVI))+
geom_point()
Step 4 - repeat steps 2&3 for the rest of the tiles and join data frames.
UEC tile
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat/Bands02")
s2=stack("T34UEC_20211029T095029_B02.jp2"," T34UEC_20211029T095029_B03.jp2"," T34UEC_20211029T095029_B04.jp2"," T34UEC_20211029T095029_B08.jp2")
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat")
waw.pow.sf2<-st_transform(waw.pow.sf, crs(s2)@projargs)
waw.sat.s<-crop(s2, extent(waw.pow.sf2))
waw.sat_przyciety<-mask(waw.sat.s, mask=waw.pow.sf2)
was<-waw.sat_przyciety
ndvi2<-(was[[4]]-was[[3]])/(was[[4]]+was[[3]])
plot(waw.sat_przyciety[[1]])
plot(ndvi2)
ndvi2.wsp<-spTransform(xyFromCell(ndvi2, spatial=TRUE, 1:5215360), CRS("+proj=longlat +datum=NAD83"))
crds2<-coordinates(ndvi2.wsp)
crds2.df<-data.frame(ID=c(1:5215360), X=crds2[,1], Y=crds2[,2])
ndvi2.wart<-as.data.frame(ndvi2)
crds2.df$NDVI<-ndvi2.wart[,1]
head(crds2.df)## ID X Y NDVI
## 1 1 20.99978 52.35043 0.5420726
## 2 2 20.99993 52.35043 0.5314183
## 3 3 21.00007 52.35043 0.5178908
## 4 4 21.00022 52.35043 0.5089035
## 5 5 21.00037 52.35043 0.4528302
## 6 6 21.00051 52.35043 0.4827769crds2.df<-crds2.df[!is.na(crds2.df$NDVI),]
ggplot(crds2.df, aes(x=X, y=Y, col=NDVI))+
geom_point()
UDD tile
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat/Bands03")
s3=stack("T34UDD_20210909T095029_B02.jp2","T34UDD_20210909T095029_B03.jp2","T34UDD_20210909T095029_B04.jp2","T34UDD_20210909T095029_B08.jp2")
setwd("C:/Users/ewado/OneDrive/Dokumenty/Studia/R-mat")
waw.pow.sf2<-st_transform(waw.pow.sf, crs(s3)@projargs)
waw.sat.s<-crop(s3, extent(waw.pow.sf2))
waw.sat_przyciety<-mask(waw.sat.s, mask=waw.pow.sf2)
was<-waw.sat_przyciety
ndvi3<-(was[[4]]-was[[3]])/(was[[4]]+was[[3]])
plot(waw.sat_przyciety[[1]])
plot(ndvi3)
ndvi3.wsp<-spTransform(xyFromCell(ndvi3, spatial=TRUE, 1:2338250), CRS("+proj=longlat +datum=NAD83"))
crds3<-coordinates(ndvi3.wsp)
crds3.df<-data.frame(ID=c(1:2338250), X=crds3[,1], Y=crds3[,2])
ndvi3.wart<-as.data.frame(ndvi3)
crds3.df$NDVI<-ndvi3.wart[,1]
head(crds3.df)## ID X Y NDVI
## 1 1 20.85114 52.36805 NA
## 2 2 20.85129 52.36805 NA
## 3 3 20.85144 52.36805 NA
## 4 4 20.85158 52.36805 NA
## 5 5 20.85173 52.36805 NA
## 6 6 20.85188 52.36805 NAcrds3.df[1000000:1000010,]## ID X Y NDVI
## 1000000 1000000 21.00081 52.32301 0.1560510
## 1000001 1000001 21.00095 52.32301 0.2024480
## 1000002 1000002 21.00110 52.32301 0.2729913
## 1000003 1000003 21.00125 52.32301 0.2508193
## 1000004 1000004 21.00139 52.32301 0.2122213
## 1000005 1000005 21.00154 52.32301 0.2384089
## 1000006 1000006 21.00169 52.32301 0.3286040
## 1000007 1000007 21.00183 52.32301 0.3052569
## 1000008 1000008 21.00198 52.32301 0.2257840
## 1000009 1000009 21.00213 52.32301 0.2270270
## 1000010 1000010 21.00227 52.32301 0.2874570crds3.df<-crds3.df[!is.na(crds3.df$NDVI),]
ggplot(crds3.df, aes(x=X, y=Y, col=NDVI))+
geom_point()
Merging data frames:
df.ca<-rbind(crds.df, crds2.df, crds3.df)
dim(df.ca)## [1] 8823978 4head(df.ca)## ID X Y NDVI
## 439 439 20.91551 52.3504 -0.1923536
## 440 440 20.91566 52.3504 -0.2149068
## 441 441 20.91580 52.3504 -0.2272175
## 442 442 20.91595 52.3504 -0.2250923
## 443 443 20.91610 52.3504 -0.2322738
## 444 444 20.91624 52.3504 -0.2307692ggplot(df.ca, aes(x=X, y=Y, col=NDVI))+
geom_point()