forked from FarmMaps/FarmMapsApiClient
115 lines
5.3 KiB
R
115 lines
5.3 KiB
R
|
# ShowGeotiff.r
|
||
|
|
# Have a look at a downloaded satellite image and check if stats are correctly calculated
|
||
|
|
# 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 '2022-03-08'
|
||
|
|
|
||
|
|
library(raster)
|
||
|
|
library(sf)
|
||
|
|
library(rgdal)
|
||
|
|
setwd("C:/git/FarmMapsApiClient_WURtest/FarmmapsDataDownload/bin/Debug/netcoreapp3.1/Downloads")
|
||
|
|
|
||
|
|
# FarmmapsDataDownload
|
||
|
|
fileGeotiff <- "sentinelhub_test_BvdTFieldlabG92_20220308.tif"
|
||
|
|
lenfilename <- nchar(fileGeotiff)
|
||
|
|
year <- substr(fileGeotiff,lenfilename-11,lenfilename-8)
|
||
|
|
imgdate <- substr(fileGeotiff,lenfilename-11,lenfilename-4)
|
||
|
|
|
||
|
|
stk.sentinelhub <- stack(x=fileGeotiff)
|
||
|
|
# plot(stk.sentinelhub) shows 6 plots (6 bands)
|
||
|
|
# 1. ndvi
|
||
|
|
# 2. wdvi Note wdvi-red
|
||
|
|
# 3. ci-red
|
||
|
|
# 4. natural: red
|
||
|
|
# 5. natural: green
|
||
|
|
# 6. natural: blue
|
||
|
|
names(stk.sentinelhub) <- c("ndvi","wdvired","ci-red","red","green","blue")
|
||
|
|
plot(stk.sentinelhub)
|
||
|
|
crs(stk.sentinelhub)
|
||
|
|
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
|
||
|
|
stk.sentinelhub.rd <- projectRaster(stk.sentinelhub, crs = CRS('+init=EPSG:28992'))
|
||
|
|
crs(stk.sentinelhub)
|
||
|
|
|
||
|
|
r.sentinelhub.rd.wdvi <- subset(stk.sentinelhub.rd,2)
|
||
|
|
dev.off()
|
||
|
|
plot(r.sentinelhub.rd.wdvi,main=paste("wdvi",imgdate),xlab="RDX",ylab="RDY")
|
||
|
|
cellStats(r.sentinelhub.rd.wdvi,'mean') # 0.2252725
|
||
|
|
|
||
|
|
# Convert the .rd.wdvi raster to WGS84
|
||
|
|
r.sentinelhub.wgs84.wdvi <- projectRaster(r.sentinelhub.rd.wdvi, crs = CRS('+init=EPSG:4326'))
|
||
|
|
|
||
|
|
# Draw a polygon on top of the raster
|
||
|
|
# Polygon pol from C:\git\FarmMapsApiClient_WURtest\FarmmapsDataDownload\DataDownloadInput.json
|
||
|
|
pol <- data.frame(id = 1, wkt = gsub("\n","",'POLYGON((
|
||
|
|
5.563472073408009 52.547554398144172,
|
||
|
|
5.567425915520115 52.547725375100377,
|
||
|
|
5.567917474269188 52.540608459298582,
|
||
|
|
5.563878143678981 52.54048022658143,
|
||
|
|
5.563472073408009 52.547554398144172
|
||
|
|
))'))
|
||
|
|
|
||
|
|
|
||
|
|
pol.wgs84 <- st_as_sf(pol, wkt = 'wkt', crs = CRS('+init=EPSG:4326'))
|
||
|
|
pol.rd <- st_transform(pol.wgs84, "+init=epsg:28992")
|
||
|
|
|
||
|
|
#Calculate approximate middle of polygon
|
||
|
|
res <- as.data.frame(do.call("rbind", lapply(st_geometry(pol.wgs84), st_bbox)))
|
||
|
|
res$latmid <- (res$ymax+res$ymin)/2.0
|
||
|
|
res$lonmid <- (res$xmax+res$xmin)/2.0
|
||
|
|
res
|
||
|
|
# xmin ymin xmax ymax latmid lonmid
|
||
|
|
# 1 5.563472 52.54048 5.567917 52.54773 52.5441 5.565695
|
||
|
|
|
||
|
|
# Have a look at both polygons
|
||
|
|
# wg84
|
||
|
|
dev.off()
|
||
|
|
plot(r.sentinelhub.wgs84.wdvi,main=paste("wdvi",imgdate),xlab="LON",ylab="LAT")
|
||
|
|
plot(pol.wgs84,add=TRUE, col="transparent",border="red")
|
||
|
|
# RD
|
||
|
|
dev.off()
|
||
|
|
plot(r.sentinelhub.rd.wdvi,main=paste("wdvi",imgdate),xlab="RDX",ylab="RDY")
|
||
|
|
plot(pol.rd,add=TRUE, col="transparent",border="red")
|
||
|
|
|
||
|
|
# Clip the polygon from the full rectangle figure
|
||
|
|
r.sentinelhub.rd.wdvi.pol <- mask(r.sentinelhub.rd.wdvi,pol.rd)
|
||
|
|
r.sentinelhub.wgs84.wdvi.pol <- mask(r.sentinelhub.wgs84.wdvi,pol.wgs84)
|
||
|
|
dev.off()
|
||
|
|
plot(r.sentinelhub.wgs84.wdvi.pol,main=paste("wdvi",imgdate),xlab="LON",ylab="LAT")
|
||
|
|
plot(pol.wgs84,add=TRUE, col="transparent",border="red")
|
||
|
|
#That's what we want!
|
||
|
|
|
||
|
|
# Now compare the stats
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi,'mean') # [1] 0.2250987 # Stats from rectangle, WGS84
|
||
|
|
cellStats(r.sentinelhub.rd.wdvi,'mean') # [1] 0.2252725 # Stats from rectangle, RD. Almost but not exactly same as above
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi.pol,'mean') # [1] 0.2275067 # Stats from raster clipped by polygon, WGS84
|
||
|
|
cellStats(r.sentinelhub.rd.wdvi.pol,'mean') # [1] 0.2275073 # Stats from raster clipped by polygon, RD. Almost but not exactly same as above
|
||
|
|
# file satelliteStats_test_BvdTFieldlabG92_20220308.csv
|
||
|
|
# "wdvi" "mean": 0.22744397204465
|
||
|
|
# Mean in csv corresponds with cellStats calculated from clipped tif!
|
||
|
|
# So while the tif returned is a non-clipped image, the downloaded statistics are from the clipped image
|
||
|
|
# Exactly as we wanted.
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi.pol,'median') # Error in .local(x, stat, ...) : invalid 'stat'. Should be sum, min, max, sd, mean, or 'countNA'
|
||
|
|
r.sentinelhub.wgs84.wdvi.vals <- values(r.sentinelhub.wgs84.wdvi)
|
||
|
|
median(r.sentinelhub.wgs84.wdvi.vals) # [1] NA
|
||
|
|
median(r.sentinelhub.wgs84.wdvi.vals,na.rm=TRUE) # [1] 0.2318
|
||
|
|
r.sentinelhub.wgs84.wdvi.pol.vals <- values(r.sentinelhub.wgs84.wdvi.pol)
|
||
|
|
median(r.sentinelhub.wgs84.wdvi.pol.vals) # [1] NA
|
||
|
|
median(r.sentinelhub.wgs84.wdvi.pol.vals,na.rm=TRUE) # [1] 0.2338
|
||
|
|
# file satelliteStats_test_BvdTFieldlabG92_20220308.csv
|
||
|
|
# "wdvi" "mean": 0.233799993991851
|
||
|
|
# Median is same as for median(r.sentinelhub.wgs84.wdvi.pol.vals,na.rm=TRUE)
|
||
|
|
# in csv corresponds with cellStats calculated from clipped tif!
|
||
|
|
# So while the tif returned is a non-clipped image, the downloaded statistics are from the clipped image
|
||
|
|
# Exactly as we wanted.
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi,'countNA') # [1] 27896
|
||
|
|
ncell(r.sentinelhub.wgs84.wdvi) # [1] 272718
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi,'countNA') / ncell(r.sentinelhub.wgs84.wdvi) # [1] 0.1022888 # 10% no data? doesn't show in the plot?
|
||
|
|
cellStats(r.sentinelhub.wgs84.wdvi.pol,'countNA') # [1] 57625
|
||
|
|
summary(r.sentinelhub.wgs84.wdvi.pol.vals) # shows the same: NA's: 57625
|
||
|
|
ncell(r.sentinelhub.wgs84.wdvi.pol) # [1] 272718
|
||
|
|
populationCount = ncell(r.sentinelhub.wgs84.wdvi.pol) - cellStats(r.sentinelhub.wgs84.wdvi.pol,'countNA')
|
||
|
|
populationCount # [1] 215093
|
||
|
|
# file satelliteStats_test_BvdTFieldlabG92_20220308.csv
|
||
|
|
# "wdvi" "populationCount": 214688
|
||
|
|
# similar but not same
|
||
|
|
|