forked from FarmMaps/FarmMapsApiClient
Merge remote-tracking branch 'upstream/master'
And updated ShowGeotiff.r using new 'r' 'terra' packages
And added comment for the DataDownloadApplication and for the BulkSatDownloadApplication:
//To run this app, first go to farmmaps datastore at https://farmmaps.eu/en/editor/plan (or on test)
//goto 'Apps and Data', goto 'Data', buy (or get for free?): 'SATELLITE'
This commit is contained in:
82
FarmmapsDataDownload/ShowGeotiff.r
Normal file
82
FarmmapsDataDownload/ShowGeotiff.r
Normal file
@@ -0,0 +1,82 @@
|
||||
# ShowGeotiff.r
|
||||
# Have a look at a downloaded satellite image and check if stats are correctly calculated
|
||||
# downloaded and calculated the stats for the polygon defined in C:\git\FarmMapsApiClient\FarmmapsDataDownload\DataDownloadInput.json
|
||||
# and in which in the console I requested the image for date '2022-08-25'
|
||||
|
||||
library(terra)
|
||||
setwd("C:/git/FarmMapsApiClient/FarmmapsDataDownload/bin/Debug/netcoreapp3.1/Downloads")
|
||||
|
||||
# FarmmapsDataDownload
|
||||
fileGeotiff <- "sentinelhub_test_BvdTFieldlabG92_20220825.tif"
|
||||
lenfilename <- nchar(fileGeotiff)
|
||||
year <- substr(fileGeotiff,lenfilename-11,lenfilename-8)
|
||||
imgdate <- substr(fileGeotiff,lenfilename-11,lenfilename-4)
|
||||
|
||||
r.sentinelhub <- rast(x=fileGeotiff)
|
||||
# plot(r.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(r.sentinelhub) <- c("ndvi","wdvired","ci-red","red","green","blue")
|
||||
# crs(r.sentinelhub)
|
||||
# CRS arguments: +proj=longlat +datum=WGS84 +no_defs
|
||||
r.sentinelhub.rd <- project(r.sentinelhub,'+init=EPSG:28992')
|
||||
# crs(r.sentinelhub.rd)
|
||||
|
||||
r.sentinelhub.rd.wdvi <- subset(r.sentinelhub.rd,2)
|
||||
# dev.off()
|
||||
plot(r.sentinelhub.rd.wdvi,main=paste("wdvi",imgdate),xlab="RDX",ylab="RDY")
|
||||
global(r.sentinelhub.rd.wdvi, fun="mean")
|
||||
|
||||
# Convert the .rd.wdvi raster to WGS84
|
||||
r.sentinelhub.wgs84.wdvi <- project(r.sentinelhub.rd.wdvi, '+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 <- vect(pol$wkt,crs='+init=EPSG:4326')
|
||||
pol.rd <- project(pol.wgs84, "+init=epsg:28992")
|
||||
|
||||
# Have a look at both polygons
|
||||
# wg84
|
||||
plot(r.sentinelhub.wgs84.wdvi,main=paste("wdvi",imgdate),xlab="LON",ylab="LAT")
|
||||
plot(pol.wgs84,add=TRUE, col="transparent",border="red")
|
||||
# RD
|
||||
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)
|
||||
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
|
||||
global(r.sentinelhub.wgs84.wdvi,'mean') # [1] 0.1733634 # Stats from rectangle, WGS84
|
||||
global(r.sentinelhub.rd.wdvi,'mean') # [1] 0.173531 # Stats from rectangle, RD. Almost but not exactly same as above
|
||||
global(r.sentinelhub.wgs84.wdvi.pol,'mean',na.rm=TRUE) # [1] 0.1740224 # Stats from raster clipped by polygon, WGS84
|
||||
global(r.sentinelhub.rd.wdvi.pol,'mean',na.rm=TRUE) # [1] 0.1738386 # Stats from raster clipped by polygon, RD. Almost but not exactly same as above
|
||||
# file satelliteStats_test_BvdTFieldlabG92_20220825.csv
|
||||
# "wdvi" "mean": 0.173905644768987
|
||||
# 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.
|
||||
global(r.sentinelhub.wgs84.wdvi.pol,'median', na.rm=TRUE) # Error in fun(values(x[[i]]), ...) : could not find function "fun"
|
||||
r.sentinelhub.wgs84.wdvi.vals <- values(r.sentinelhub.wgs84.wdvi)
|
||||
median(r.sentinelhub.wgs84.wdvi.vals) # [1] 0.1459627
|
||||
median(r.sentinelhub.wgs84.wdvi.vals,na.rm=TRUE) # [1] 0.1459627
|
||||
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.1453323
|
||||
Reference in New Issue
Block a user