From 334d890f1353a8ce3fffa4ec58c1b34e309cae97 Mon Sep 17 00:00:00 2001 From: Pepijn van Oort Date: Fri, 11 Mar 2022 15:10:33 +0100 Subject: [PATCH] update DataDownloadApplication --- FarmmapsApi/Constants.cs | 2 +- FarmmapsBulkSatDownload/ShowGeotiff.r | 178 ------------------ .../DataDownloadApplication.cs | 125 +++++++++--- FarmmapsDataDownload/DataDownloadInput.json | 32 ++-- .../FarmmapsDataDownload.csproj | 2 +- .../Models/DataDownloadInput.cs | 3 +- FarmmapsDataDownload/ShowGeotiff.r | 114 +++++++++++ 7 files changed, 232 insertions(+), 224 deletions(-) delete mode 100644 FarmmapsBulkSatDownload/ShowGeotiff.r create mode 100644 FarmmapsDataDownload/ShowGeotiff.r diff --git a/FarmmapsApi/Constants.cs b/FarmmapsApi/Constants.cs index 1546edb..05ac234 100644 --- a/FarmmapsApi/Constants.cs +++ b/FarmmapsApi/Constants.cs @@ -19,7 +19,7 @@ namespace FarmmapsApiSamples public const string VRAHAULMKILLING_TASK = "vnd.farmmaps.task.vrahaulmkilling"; public const string VRAPLANTING_TASK = "vnd.farmmaps.task.vrapoten"; public const string VRAZONERING_TASK = "vnd.farmmaps.task.vrazonering"; - public const string SATELLITE_TASK = "vnd.farmmaps.task.satellite"; + public const string SATELLITE_TASK = "vnd.farmmaps.task.sentinelhub"; //"vnd.farmmaps.task.satellite"; public const string VANDERSAT_TASK = "vnd.farmmaps.task.vandersat"; public const string TASKMAP_TASK = "vnd.farmmaps.task.taskmap"; public const string WORKFLOW_TASK = "vnd.farmmaps.task.workflow"; diff --git a/FarmmapsBulkSatDownload/ShowGeotiff.r b/FarmmapsBulkSatDownload/ShowGeotiff.r deleted file mode 100644 index bd28e92..0000000 --- a/FarmmapsBulkSatDownload/ShowGeotiff.r +++ /dev/null @@ -1,178 +0,0 @@ -# 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 diff --git a/FarmmapsDataDownload/DataDownloadApplication.cs b/FarmmapsDataDownload/DataDownloadApplication.cs index e4534d0..a1d73e9 100644 --- a/FarmmapsDataDownload/DataDownloadApplication.cs +++ b/FarmmapsDataDownload/DataDownloadApplication.cs @@ -38,7 +38,7 @@ namespace FarmmapsDataDownload public async Task RunAsync() { - var fieldsInputJson = File.ReadAllText("DataDownloadInput.json"); + string fieldsInputJson = File.ReadAllText("DataDownloadInput.json"); List fieldsInputs = JsonConvert.DeserializeObject>(fieldsInputJson); @@ -61,6 +61,28 @@ namespace FarmmapsDataDownload private async Task Process(List roots, DataDownloadInput input) { + //PO20220311: first time a call is made to download satellite images or statistics, an empty list is returned + //If we wait a bit longer, e.g. 10 secs, then e.g. a list of 3 images may be returned + //If we wait still longer, maybe 4 images. + //The solution implemented below is to fire calls as long as the number of images returned keeps increasing + //While in between each call, sleep for sleepSecs + //Continue this until the number no longer increases or the maximum number of calls has been reached + //Out of politeness, don't be too impatient. Don't set sleepSecs to 5 or 10 or 30 secs. Just accept this may take a while, have a coffee, we suggest sleepSecs = 60; + int sleepSecs = 60; + int callCntMax = 4; + //For example we may set: "sleepSecs = 10;" and "callCntMax = 24;" and following result: + //Call no: 1. Giving FarmMaps 10 seconds to get SatelliteItems... + //Call no: 1: Received 2 images + //Call no: 2. Giving FarmMaps 10 seconds to get SatelliteItems... + //Call no: 2: Received 7 images + //Call no: 3. Giving FarmMaps 10 seconds to get SatelliteItems... + //Call no: 3: Received 7 images + //And the firing of calls would stop because the number of images returned is no longer increasing + //In the worst case, this could would lead to a total sleeping period of "sleepSecsSum = sleepSecs * callCntMax" seconds. After that we give up + + //This is an ugly fix. Neater would if FarmMaps would just take a bit longer and then do always deliver all satellite images on first call. + //Once this has been fixed on the side of FarmMaps we can set callCntMax = 0 and the code below will work smoothly without any sleeping + string downloadFolder = input.DownloadFolder; if (string.IsNullOrEmpty(downloadFolder)) { downloadFolder = "Downloads"; @@ -74,7 +96,8 @@ namespace FarmmapsDataDownload var fieldName = input.fieldName; bool storeSatelliteStatistics = input.StoreSatelliteStatisticsSingleImage; bool storeSatelliteStatisticsCropYear = input.StoreSatelliteStatisticsCropYear; - List SatelliteBands = new List(1) { input.SatelliteBand }; + //List SatelliteBands = new List(1) { input.SatelliteBand }; + List satelliteBands = input.SatelliteBands; string headerLineStats = $"FieldName,satelliteDate,satelliteBand,max,min,mean,mode,median,stddev,minPlus,curtosis,maxMinus,skewness,variance,populationCount,variationCoefficient,confidenceIntervalLow, confidenceIntervalHigh,confidenceIntervalErrorMargin" + Environment.NewLine; @@ -166,12 +189,58 @@ namespace FarmmapsDataDownload SaveSettings(settingsfile); } - - // Select all satellite items + //Call first time + int callCnt = 1; + int sleepSecsSum = 0; + //if callCntMax == 0 then don't sleep + //if callCntMax = 1 then sleep first 1x + if (callCntMax > 0) + { + _logger.LogInformation($"Call no: {callCnt}. Giving FarmMaps {sleepSecs} seconds to get SatelliteItems..."); + System.Threading.Thread.Sleep(1000 * sleepSecs); + sleepSecsSum = sleepSecsSum + sleepSecs; + } List satelliteItemsCropYear = await _generalService.FindSatelliteItems(cropfieldItem, _settings.SatelliteTaskCode); + int satelliteItemsCropYearCntPrev = satelliteItemsCropYear.Count; + _logger.LogInformation($"Call no: {callCnt}. Received {satelliteItemsCropYearCntPrev} images"); + callCnt++; + int satelliteItemsCropYearCnt = satelliteItemsCropYearCntPrev; + //if callCntMax > 1 then sleep untill (1) no more increase in number of images received OR (2) maximum number of calls reached + if (callCntMax > 1) + { + //Call second time + _logger.LogInformation($"Call no: {callCnt}. Giving FarmMaps another {sleepSecs} seconds to get SatelliteItems..."); + System.Threading.Thread.Sleep(1000 * sleepSecs); + satelliteItemsCropYear = await _generalService.FindSatelliteItems(cropfieldItem, _settings.SatelliteTaskCode); + satelliteItemsCropYearCnt = satelliteItemsCropYear.Count; + _logger.LogInformation($"Call no: {callCnt}. Received {satelliteItemsCropYearCnt} images"); + sleepSecsSum = sleepSecsSum + sleepSecs; + //As long as there is progress, keep calling + callCnt++; + while (callCnt <= callCntMax && (satelliteItemsCropYearCnt == 0 || satelliteItemsCropYearCnt > satelliteItemsCropYearCntPrev)) + { + _logger.LogInformation($"Surprise! The longer we wait, the more images we get. Sleep and call once more"); + satelliteItemsCropYearCntPrev = satelliteItemsCropYearCnt; + _logger.LogInformation($"Call no: {callCnt} (max: {callCntMax}). Giving FarmMaps another {sleepSecs} seconds to get SatelliteItems..."); + System.Threading.Thread.Sleep(1000 * sleepSecs); + satelliteItemsCropYear = await _generalService.FindSatelliteItems(cropfieldItem, _settings.SatelliteTaskCode); + satelliteItemsCropYearCnt = satelliteItemsCropYear.Count; + _logger.LogInformation($"Call no: {callCnt}. Received {satelliteItemsCropYearCnt} images"); + callCnt++; + sleepSecsSum = sleepSecsSum + sleepSecs; + } + } + + if (satelliteItemsCropYearCnt == 0) + { + _logger.LogWarning($"DataDownloadApplication.cs: after calling one or more times and " + + $"sleeping in total {sleepSecsSum} seconds, still no images found. " + + $"Please check your settings for parameters callCntMax and sleepSecs in DataDownloadApplication.cs or contact FarmMaps"); + } + satelliteItemsCropYear = satelliteItemsCropYear.OrderBy(x => x.DataDate).ToList(); - if (input.StoreSatelliteStatisticsSingleImage == true) { + if (input.StoreSatelliteStatisticsSingleImage == true && satelliteItemsCropYearCnt > 0) { _logger.LogInformation("Available satellite images:"); var count = 0; TimeSpan.FromSeconds(0.5); @@ -188,43 +257,49 @@ namespace FarmmapsDataDownload var SatelliteDate = selectedSatelliteItem.DataDate.Value.ToString("yyyyMMdd"); string fileName = string.Format($"satelliteGeotiff_{fieldName}_{SatelliteDate}"); // no need to add satelliteBand in the name because the tif contains all bands - string fileNameZip = string.Format($"{fileName}.zip"); - string fileNameGeotiff = string.Format($"{fileName}.tif"); - await _farmmapsApiService.DownloadItemAsync(selectedSatelliteItem.Code, Path.Combine(downloadFolder, fileNameZip)); + string fileNameZip = Path.Combine(downloadFolder, string.Format($"{fileName}.zip")); + await _farmmapsApiService.DownloadItemAsync(selectedSatelliteItem.Code, fileNameZip); // Download a csv file with stats - List selectedSatalliteItems = new List(1) { selectedSatelliteItem }; + List selectedSatelliteItems = new List(1) { selectedSatelliteItem }; string fileNameStats = Path.Combine(downloadFolder, string.Format($"satelliteStats_{fieldName}_{SatelliteDate}.csv")); - string downloadedStats = await _generalService.DownloadSatelliteStats(selectedSatalliteItems, fieldName, SatelliteBands, downloadFolder); + _logger.LogInformation($"First call to get DownloadSatelliteStats for selected image..."); + string downloadedStats = await _generalService.DownloadSatelliteStats(selectedSatelliteItems, fieldName, satelliteBands, downloadFolder); + //rename the csv file with stats + //if the targe file already exists, delete it File.Delete(fileNameStats); + //rename File.Move(downloadedStats, fileNameStats); - // wenr.tif. Contains 5 layers: (1) ndvi, (2) wdvi, (3) Red, (4) Green and (5) Blue - // download the geotiffs. Returns a zip file with always these three files: - // data.dat.aux.xml + + // name the tif file + string fileNameTifzipped = Path.Combine(downloadFolder, string.Format($"sentinelhub_{SatelliteDate}.tif")); + string fileNameGeotiff = Path.Combine(downloadFolder, string.Format($"sentinelhub_{fieldName}_{SatelliteDate}.tif")); + // download the geotiffs. Returns a zip file with always these two files: // thumbnail.jpg - // wenr.tif. Contains 5 layers: (1) ndvi, (2) wdvi, (3) Red, (4) Green and (5) Blue + // sentinelhub_yyyyMMdd.tif. Contains 4 layers: (1) ndvi, (2) wdvi, (3) ci-red and (4) natural. Natural has 3 layers inside: redBand, blueBand and greenBand if (true) { - // Extract the file "wenr.tif" from zip, rename it to fileNameGeotiff - ZipFile.ExtractToDirectory(Path.Combine(downloadFolder, fileNameZip), downloadFolder, true); - File.Delete(Path.Combine(downloadFolder, fileNameGeotiff)); // Delete the fileNameGeotiff file if exists - File.Move(Path.Combine(downloadFolder, "wenr.tif"), Path.Combine(downloadFolder, fileNameGeotiff)); // Rename the oldFileName into newFileName + // Extract the file fileNameTifzipped from zip, rename it to fileNameGeotiff + ZipFile.ExtractToDirectory(fileNameZip, downloadFolder, true); + //if the targe file already exists, delete it + File.Delete(fileNameGeotiff); + //rename + File.Move(fileNameTifzipped, fileNameGeotiff); // Cleanup - string[] filesToDelete = new string[] { fileNameZip, "wenr.tif", "thumbnail.jpg", "data.dat.aux.xml" }; - foreach (string f in filesToDelete) - { - File.Delete(Path.Combine(downloadFolder, f)); - } + File.Delete(fileNameZip); + File.Delete(Path.Combine(downloadFolder, "thumbnail.jpg")); } - _logger.LogInformation($"Downloaded files {fileNameGeotiff} and {fileNameStats} to {downloadFolder}"); + //_logger.LogInformation($"Downloaded files {fileNameGeotiff} and {fileNameStats} to {downloadFolder}"); + _logger.LogInformation($"Downloaded files to {downloadFolder}"); } if (input.StoreSatelliteStatisticsCropYear == true) { string fileNameStats = Path.Combine(downloadFolder, string.Format($"satelliteStats_{fieldName}_{cropYear}.csv")); File.Delete(fileNameStats); - string downloadedStats = await _generalService.DownloadSatelliteStats(satelliteItemsCropYear, fieldName, SatelliteBands, downloadFolder); + _logger.LogInformation($"First call to get DownloadSatelliteStats for whole cropYear..."); + string downloadedStats = await _generalService.DownloadSatelliteStats(satelliteItemsCropYear, fieldName, satelliteBands, downloadFolder); File.Move(downloadedStats, fileNameStats); _logger.LogInformation($"Downloaded file {fileNameStats} with stats for field '{fieldName}', cropyear {cropYear}"); } diff --git a/FarmmapsDataDownload/DataDownloadInput.json b/FarmmapsDataDownload/DataDownloadInput.json index a675294..e5b124e 100644 --- a/FarmmapsDataDownload/DataDownloadInput.json +++ b/FarmmapsDataDownload/DataDownloadInput.json @@ -1,34 +1,30 @@ [ { - "UseCreatedCropfield": true, - "outputFileName": "TestData", - "fieldName": "TestField", + "UseCreatedCropfield": false, // if false, program will make new CropfieldItemCode and SatelliteTaskCode; if true, the program will read CropfieldItemCode and SatelliteTaskCode from a file called "..\FarmmapsDataDownload\bin\Debug\netcoreapp3.1\Settings_{fieldName}.json", which will be faster + "outputFileName": "test_BvdTFieldlabG92", + "fieldName": "test_BvdTFieldlabG92", "DownloadFolder": "Downloads", //"C:\\workdir\\groenmonitor\\", // "Downloads", -> if you just put "Downloads" the program will download to somewhere in ..\FarmMapsApiClient_WURtest\FarmmapsDataDownload\bin\Debug\netcoreapp3.1\Downloads\ - "GetCropRecordings": true, + "GetCropRecordings": false, "CrprecItem": "...", //item code of de crop recording parrent - can be found by opening the crop recording page of a field. "GetShadowData": false, - "GetSatelliteData": false, - "SatelliteBand": "wdvi", // "natural", "ndvi" or "wdvi" - "StoreSatelliteStatisticsSingleImage": false, - "StoreSatelliteStatisticsCropYear": false, + "GetSatelliteData": true, + "SatelliteBands": [ "ndvi", "wdvi", "ci-red" ], // ["ndvi"] or ["wdvi"] or ["ci-red"] or multiple: [ "wdvi", "ndvi" ] + "StoreSatelliteStatisticsSingleImage": true, + "StoreSatelliteStatisticsCropYear": true, "GetVanDerSatData": false, "StoreVanDerSatStatistics": false, - "CropYear": 2020, + "CropYear": 2022, "geometryJson": { "type": "Polygon", "coordinates": [ [ - [ 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 ] + [ 5.563472073408009, 52.547554398144172 ], + [ 5.567425915520115, 52.547725375100377 ], + [ 5.567917474269188, 52.540608459298582 ], + [ 5.563878143678981, 52.54048022658143 ], + [ 5.563472073408009, 52.547554398144172 ] ] ] } } - - ] \ No newline at end of file diff --git a/FarmmapsDataDownload/FarmmapsDataDownload.csproj b/FarmmapsDataDownload/FarmmapsDataDownload.csproj index 22b003c..173a795 100644 --- a/FarmmapsDataDownload/FarmmapsDataDownload.csproj +++ b/FarmmapsDataDownload/FarmmapsDataDownload.csproj @@ -2,7 +2,7 @@ Exe - netcoreapp3.0 + netcoreapp3.1 diff --git a/FarmmapsDataDownload/Models/DataDownloadInput.cs b/FarmmapsDataDownload/Models/DataDownloadInput.cs index 19f1bd1..1912c1a 100644 --- a/FarmmapsDataDownload/Models/DataDownloadInput.cs +++ b/FarmmapsDataDownload/Models/DataDownloadInput.cs @@ -1,4 +1,5 @@ using System; +using System.Collections.Generic; using Newtonsoft.Json.Linq; namespace FarmmapsDataDownload.Models @@ -16,7 +17,7 @@ namespace FarmmapsDataDownload.Models public string fieldName { get; set; } public bool GetSatelliteData { get; set; } public bool GetVanDerSatData { get; set; } - public string SatelliteBand { get; set; } + public List SatelliteBands { get; set; } public bool StoreSatelliteStatisticsSingleImage { get; set; } public bool StoreSatelliteStatisticsCropYear { get; set; } public bool StoreVanDerSatStatistics { get; set; } diff --git a/FarmmapsDataDownload/ShowGeotiff.r b/FarmmapsDataDownload/ShowGeotiff.r new file mode 100644 index 0000000..1ce3d8d --- /dev/null +++ b/FarmmapsDataDownload/ShowGeotiff.r @@ -0,0 +1,114 @@ +# 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 +