FarmMapsApiClient_KB34_MAST/FarmmapsBulkSatDownload/ShowGeotiff.r

# ShowGeotiff.r
# I downloaded and calculated the stats for the polygon defined in C:\git\FarmMapsApiClient_WURtest\FarmmapsDataDownload\DataDownloadInput.json
# in which I set "SatelliteBand": "wdvi" and in which in the console I requested the image for date '2020-09-22'
# FarmmapsBulkSatDownload generates many files. Here is what I tried when inputing the same field for a number of years using BulkSatDownloadInput.json
# see list below

library(raster)
library(sf)
library(rgdal)
setwd("C:/workdir/groenmonitor/DataDownload/")

# FarmmapsDataDownload and BulkSatDownload can be used to download zip files with inside two files:
fileGeotiff <- "wenr.tif"
fileJpg <- "thumbnail.jpg"
# Here is what I tried when inputing the same field for a number of years:
# fileGeotiff <- "wheat_fld5641_20210224.tif"	# 2 files for year 2021. This file is a nice example having in upperleft corner no data
# fileGeotiff <- "wheat_fld5641_20210331.tif"	# 2 files for year 2021
# fileGeotiff <- "wheat_fld5641_20200321.tif"	# 14 files for year 2020, earliest
# fileGeotiff <- "wheat_fld5641_20200922.tif"	# 14 files for year 2020, latest
# fileGeotiff <- "wheat_fld5641_20190121.tif"	# 9 files for year 2019, earliest
# fileGeotiff <- "wheat_fld5641_20191117.tif"	# 9 files for year 2019, latest
# 1 file for year 2018, with error message 'End of Central Directory record could not be found' and invalid wheat_fld5641_20180630.zip 
# fileGeotiff <- "wheat_fld5641_20170526.tif"	# 1 file for year 2017
# Zero files for 2016
lenfilename <- nchar(fileGeotiff)
year <- substr(fileGeotiff,lenfilename-11,lenfilename-8)
imgdate <- substr(fileGeotiff,lenfilename-11,lenfilename-4)

# The thumbnail has the polygon clipped out, has 1 layer, no crs and the mean value is not the mean wdvi we are looking for
r.thumbnail <- raster(fileJpg)
plot(r.thumbnail)
crs(r.thumbnail)
#CRS arguments: NA 
cellStats(r.thumbnail,'mean')	#87.5128 # nonsense

stk.wenr <- stack(x=fileGeotiff)
# plot(stk.wenr) shows 5 plots (5 bands)
# I think these are:
# 1. ndvi (since it runs from 0 to 1)
# 2. wdvi (since it runs from 0 to 0.5)
# 3-5: RGB (since they run from 0 to 255)
plot(stk.wenr)
# CRS arguments:
# +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
# Or use st_crs(stk.wenr) to get more info, but no EPSG code in there. 
# Likely it is epsg:28992 (Amersfoort)
crs(stk.wenr)
stk.wenr <- projectRaster(stk.wenr, crs = CRS('+init=EPSG:28992'))
crs(stk.wenr)
# Looks the same but strangely, if we don't do projectRaster(stk.wenr, crs = CRS('+init=EPSG:28992')), we find below the bottom left corner of the polygon missing

r.wenr.rd.wdvi <- subset(stk.wenr,2)
dev.off()
plot(r.wenr.rd.wdvi,main=paste("wdvi",imgdate),xlab="RDX",ylab="RDY")
cellStats(r.wenr.rd.wdvi,'mean')	#0.1350561

# Furthermore we can see
# shows coordinates in RD
# returns a rectangle, thus the shape of the polygon submitted is not clipped. 
# The polygon was provided in WGS84. Let's draw it on top. 
# First convert the raster to WGS84
r.wenr.wgs84.wdvi <- projectRaster(r.wenr.rd.wdvi, crs = CRS('+init=EPSG:4326'))

# Draw a polygon on top of the raster
# Example polygon p1
# p1 <- data.frame(id = 1, wkt = 'POLYGON((4.963 52.801, 4.966 52.801, 4.966 52.803, 4.963 52.803, 4.963 52.801))')
# p1 <- st_as_sf(p1, wkt = 'wkt', crs = targetcrs)
# plot(p1,add=TRUE, col="transparent",border="black")
# Draw the polygon on top of the raster
# Polygon p2 from C:\git\FarmMapsApiClient_WURtest\FarmmapsDataDownload\DataDownloadInput.json
p2 <- data.frame(id = 1, wkt = gsub("\n","",'POLYGON((
          4.960707146896585 52.800583669708487,
          4.960645975538824 52.800470217610922,
          4.962140695752897 52.799177147194797,
          4.967523821195745 52.801502400041208,
          4.966336768950911 52.802543735879809,
          4.961711880764330 52.801009996856429,
          4.960707146896585 52.800583669708487))'))
# Polygon p2 from C:\git\FarmMapsApiClient_WURtest\FarmmapsBulkSatDownload\BulkSatDownloadInput.json
p2 <- data.frame(id = 1, wkt = gsub("\n","",'POLYGON((
			3.37837807779104   51.3231095796538,
			3.38065689232502   51.3212527499355,
			3.38022924592256   51.3210683536359,
			3.37980548452565   51.3208801127141,
			3.37959556105776   51.3207540143696,
			3.3793691292654   51.3205959677371,
			3.37822219207335   51.3215667913007,
			3.37816999925795   51.3216109809456,
			3.37646704574705   51.3208025481261,
			3.37646695791282   51.3208025061493,
			3.37608401443192   51.3206231652693,
			3.37607169507628   51.3206173959751,
			3.37606021048754   51.320612017601,
			3.37582728410659   51.3205029306946,
			3.37580409779263   51.3206502985963,
			3.37575872019649   51.3207993094705,
			3.37575476634361   51.3208122883487,
			3.37571181656268   51.3208797459348,
			3.3756624532907   51.3209415238446,
			3.37557609963811   51.3210110142077,
			3.37541089899821   51.3211055871218,
			3.37477516102591   51.3214102985009,
			3.37473173914127   51.3214311108204,
			3.37455904622072   51.3215138815012,
			3.37415098054777   51.3217199232877,
			3.37313700916272   51.3222422862785,
			3.37748824689601   51.3242852920348,
			3.37749760805371   51.3242713084009,
			3.37811903757028   51.3233437635596,
			3.37818758851947   51.3232647797363,
			3.37823803668144   51.3232236798646,
			3.37837807779104   51.3231095796538))'))
p2.wgs84 <- st_as_sf(p2, wkt = 'wkt', crs = CRS('+init=EPSG:4326'))
# Or other way round, in RD Amersfoort. That looks ok
p2.rd <- st_transform(p2.wgs84, "+init=epsg:28992")

# Have a look at both
# wg84
dev.off()
plot(r.wenr.wgs84.wdvi,main=paste("wdvi",imgdate),xlab="LON",ylab="LAT")
plot(p2.wgs84,add=TRUE, col="transparent",border="red")
# RD
dev.off()
plot(r.wenr.rd.wdvi,main=paste("wdvi",imgdate),xlab="RDX",ylab="RDY")
plot(p2.rd,add=TRUE, col="transparent",border="red")

#Let's clip the polygon
r.wenr.rd.wdvi.pol <- mask(r.wenr.rd.wdvi,p2.rd)
r.wenr.wgs84.wdvi.pol <- mask(r.wenr.wgs84.wdvi,p2.wgs84)
dev.off()
plot(r.wenr.wgs84.wdvi.pol,main=paste("wdvi",imgdate),xlab="LON",ylab="LAT")
plot(p2.wgs84,add=TRUE, col="transparent",border="red")
#That's what we want! Now compare the stats
cellStats(r.wenr.rd.wdvi,'mean')	# [1] 0.1350561 # Stats from rectangle, RD
cellStats(r.wenr.wgs84.wdvi,'mean')	# [1] 0.1351411 # Stats from rectangle, WGS84
cellStats(r.wenr.rd.wdvi.pol,'mean')	# [1] 0.05723957 # Stats from raster clipped by polygon, RD
cellStats(r.wenr.wgs84.wdvi.pol,'mean')	# [1] 0.05723607 # Stats from raster clipped by polygon, WGS84
# file SatelliteDataStatistics_test_satData_wdvi_2020-09-22.csv
# "mean": 0.057430520945401985	# SatelliteDataStatistics_test_satData_wdvi_2020.csv returns stats for the clipped raster (.pol). 'mean' almost the same, maybe 
# cellStats cannot return median, just a few stats.
cellStats(r.wenr.wgs84.wdvi.pol,'median')	# Error in .local(x, stat, ...) : invalid 'stat'. Should be sum, min, max, sd, mean, or 'countNA'
r.wenr.wgs84.wdvi.vals <- values(r.wenr.wgs84.wdvi)
median(r.wenr.wgs84.wdvi.vals)	# [1] NA
median(r.wenr.wgs84.wdvi.vals,na.rm=TRUE)	# [1] 0.076
r.wenr.wgs84.wdvi.pol.vals <- values(r.wenr.wgs84.wdvi.pol)
median(r.wenr.wgs84.wdvi.pol.vals)	# [1] NA
median(r.wenr.wgs84.wdvi.pol.vals,na.rm=TRUE)	# [1] 0.048
# "median": 0.04800000041723251	# SatelliteDataStatistics_test_satData_wdvi_2020.csv returns stats for the clipped raster (.pol). 
# An image may contain NA values. Check:
cellStats(r.wenr.wgs84.wdvi,'countNA')	# [1] 22956
ncell(r.wenr.wgs84.wdvi)	# [1] 221696
cellStats(r.wenr.wgs84.wdvi,'countNA') / ncell(r.wenr.wgs84.wdvi)	# [1] 0.1035472	# 10% no data? doesn't show in the plot?
cellStats(r.wenr.wgs84.wdvi.pol,'countNA')	# [1] 147387
summary(r.wenr.wgs84.wdvi.pol.vals)	# shows the same: NA's: 147387
ncell(r.wenr.wgs84.wdvi.pol)	# [1] 221696
cellStats(r.wenr.wgs84.wdvi.pol,'countNA') / ncell(r.wenr.wgs84.wdvi.pol)	# [1] 0.6648158	# 66% no data? doesn't show in the plot?

# The project FarmmapsNbs can generate a wenr.tif file, application.tif, uptake.tif (in rtest1.uptake.zip)and shape.shp (in rtest1.taskmap.zip)
r.application <- raster("C:/git/FarmMapsApiClient_WURtest/FarmmapsNbs/bin/Debug/netcoreapp3.1/Downloads/application.tif")
dev.off()
plot(r.application)
plot(p2.rd,add=TRUE, col="transparent",border="red")
# The application.tif file is a rectangle (polygon not yet clipped), in projection Amersfoort RD New (EPSG:28992)

r.uptake <- raster("C:/git/FarmMapsApiClient_WURtest/FarmmapsNbs/bin/Debug/netcoreapp3.1/Downloads/uptake.tif")
dev.off()
plot(r.uptake)
plot(p2.rd,add=TRUE, col="transparent",border="red")
# The uptake.tif file is a rectangle (polygon not yet clipped), in projection Amersfoort RD New (EPSG:28992)

shp.wgs84 <- readOGR(dsn="C:/git/FarmMapsApiClient_WURtest/FarmmapsNbs/bin/Debug/netcoreapp3.1/Downloads", layer="shape")
crs(shp.wgs84)
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
dev.off()
plot(r.wenr.wgs84.wdvi,main="wdvi",xlab="LON",ylab="LAT")
plot(shp.wgs84,add=TRUE, col="transparent",border="black")
plot(p2.wgs84,add=TRUE, col="transparent",border="red")
# The shape file is in WGS84