1. Load Required Packages

Load the packages that will be used in this analysis:

#Required to convert jp2 to tif
require(gdalUtils)

#Need to read tif as raster
require(raster)

#Required for doing transformation and reading KML
require(sf)

2. Converting Sentinel Images from JP2 to TIF

Due to raster package in R not supporting jp2000, the downloaded sentinel images were converted to tif.

Band 4 (red) and band 8 (NIR) were converted since they’re ones that will be used to calculate the NDVI.

The downloaded images were for 25-05-2019

#Convert red image
gdal_translate("S2A_MSIL1C_20190525T112121_N0207_R037_T29SND_20190525T132011.SAFE/GRANULE/L1C_T29SND_A020483_20190525T112609/IMG_DATA/T29SND_20190525T112121_B04.jp2", "B04.tif")
## NULL
#Convert NIR image
gdal_translate("S2A_MSIL1C_20190525T112121_N0207_R037_T29SND_20190525T132011.SAFE/GRANULE/L1C_T29SND_A020483_20190525T112609/IMG_DATA/T29SND_20190525T112121_B08.jp2", "B08.tif")
## NULL

3. Read TIF as Raster

To manipulate and extract information, the tif files are read as raster objects in R:

#Read red image
imgred4<-raster("B04.tif")

#Read NIR image
imgnir8<-raster("B08.tif")

In order to visualise what we’re dealing with, we’ll plot both images:

plot(imgred4, main="Red Band - band 4")

plot(imgnir8, main="NIR Band - band 8")

4. Compute and Create NDVI Image:

As per definition, the NDVI is:

\[NDVI=\frac{NIR-Red}{NIR+Red}\]

Now we’ll create an NDVI image and plot it:

#Calculate NDVI and plot it
ndvi<-(imgnir8-imgred4)/(imgnir8+imgred4)

plot(ndvi, main="NDVI image")

5. Import KML File

To pin point the areas in question, the KML file will be used to read the coordinates:

#Read KML file
kmlpoints<-st_read("pontos.kml")
## Reading layer `Pontos_Limpa_e_Aduba' from data source `/Users/omarnooreddin/Documents/GitHub/Fire Fly challenge/pontos.kml' using driver `KML'
## Simple feature collection with 5 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: -8.352288 ymin: 39.49594 xmax: -8.282165 ymax: 39.51637
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
#Show points in quesion
kmlpoints$geometry
## Geometry set for 5 features 
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: -8.352288 ymin: 39.49594 xmax: -8.282165 ymax: 39.51637
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## POINT Z (-8.288045 39.51637 -4.235562e-05)
## POINT Z (-8.352288 39.49939 -4.232582e-05)
## POINT Z (-8.289838 39.50874 -4.234072e-05)
## POINT Z (-8.28952 39.50381 -4.233327e-05)
## POINT Z (-8.282165 39.49594 -4.231837e-05)

6. Convert Projections into UTM

If we inspect the KML file, we’ll see that it uses a “latlong” projection, as per above’s output.

While the Sentinel 2 images use a UTM projection, as per below:

ndvi
## class      : RasterLayer 
## dimensions : 10980, 10980, 120560400  (nrow, ncol, ncell)
## resolution : 10, 10  (x, y)
## extent     : 499980, 609780, 4290240, 4400040  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /private/var/folders/jc/h9v1lxtj4z3b50lm5tgbfp140000gn/T/RtmpSGD88X/raster/r_tmp_2019-05-28_011900_12215_08760.grd 
## names      : layer 
## values     : -0.7224607, 0.9994402  (min, max)

As such to extract the pixels from the correct coordinates we’re going to:

#Extract coordinates from KML and create a dataframe
xy<-data.frame(x=c(-8.288045, -8.352288, -8.289838, -8.28952, -8.282165), 
               y=c( 39.51637,  39.49939,  39.50874,  39.50381, 39.49594 ))

#Convert projection of these points from latlong to UTM
coordinates(xy) <- c("x", "y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84")
utmcoord <- spTransform(xy, CRS("+proj=utm +zone=29 ellps=WGS84"))
utmcoord
## class       : SpatialPoints 
## features    : 5 
## extent      : 555690.4, 561722.6, 4372060, 4374323  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=29 ellps=WGS84 +ellps=WGS84

And just to confirm that these points are indeed present in the image we downloaded, we’ll plot the NDVI image with the points superimposed on them:

#Plot NDVI and add points of interest
plot(ndvi)
points(utmcoord)

7. Extract NDVI Values for Points of Interest

Using the converted coordinates, we’ll extract the NDVI values of these points:

ndvi25MAY<-extract(ndvi, utmcoord)
ndvi25MAY
## [1] 0.3311726 0.3968167 0.2025846 0.2826562 0.3899876

The above shows NDVI values for each provided coordinates, in the same order. i.e 1st element is point 1, 2nd element is point 2 etc…

8. Calculating the Change

We’ll repeat the above steps but for images taken on the 15th of May, and compare it with images taken on the 25th of May (they seem to have skipped the 20th!):

#Convert red image
gdal_translate("S2A_MSIL1C_20190515T112121_N0207_R037_T29SND_20190515T132120.SAFE/GRANULE/L1C_T29SND_A020340_20190515T113007/IMG_DATA/T29SND_20190515T112121_B04.jp2", "B04-1505.tif")
## NULL
#Convert NIR image
gdal_translate("S2A_MSIL1C_20190515T112121_N0207_R037_T29SND_20190515T132120.SAFE/GRANULE/L1C_T29SND_A020340_20190515T113007/IMG_DATA/T29SND_20190515T112121_B08.jp2", "B08-1505.tif")
## NULL
#Read red image
imgred4_15may<-raster("B04-1505.tif")

#Read NIR image
imgnir8_15may<-raster("B08-1505.tif")

#Calculate NDVI
ndvi_15may<-(imgnir8_15may - imgred4_15may)/(imgnir8_15may + imgred4_15may)

#Extract
ndvi15MAY<-extract(ndvi_15may, utmcoord)
ndvi15MAY
## [1] 0.3178066 0.4292224 0.1969309 0.3111055 0.4237710

Having calculated the NDVI 10 days earlier on the 15th of May, we now can compare with the 25th of May and see whether the index has increased or decreased for each point.

((ndvi25MAY-ndvi15MAY)/ndvi15MAY)*100
## [1]  4.205716 -7.549877  2.870864 -9.144600 -7.972076

From the above, we can see if the vegetation has decreased (-ve change) or increased (+ve change) for each coordinate:

Area Change
1 +ve
2 -ve
3 +ve
4 -ve
5 -ve