pvanoort/FarmMapsApiClient_WUR
1
0
Fork 0
forked from FarmMaps/FarmMapsApiClient
Code Pull Requests Activity

update DataDownloadApplication

Browse Source
This commit is contained in:
Pepijn van Oort
2022-03-11 15:10:33 +01:00
parent f8a412be64
commit 334d890f13
7 changed files with 232 additions and 224 deletions
Download Patch File Download Diff File Expand all files Collapse all files

125
FarmmapsDataDownload/DataDownloadApplication.cs
View File

@@ -38,7 +38,7 @@ namespace FarmmapsDataDownload
public async Task RunAsync()
{
var fieldsInputJson = File.ReadAllText("DataDownloadInput.json");
string fieldsInputJson = File.ReadAllText("DataDownloadInput.json");
List<DataDownloadInput> fieldsInputs = JsonConvert.DeserializeObject<List<DataDownloadInput>>(fieldsInputJson);
@@ -61,6 +61,28 @@ namespace FarmmapsDataDownload
private async Task Process(List<UserRoot> 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<string> SatelliteBands = new List<string>(1) { input.SatelliteBand };
//List<string> SatelliteBands = new List<string>(1) { input.SatelliteBand };
List<string> 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<Item> 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<Item> selectedSatalliteItems = new List<Item>(1) { selectedSatelliteItem };
List<Item> selectedSatelliteItems = new List<Item>(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}");
}

32
FarmmapsDataDownload/DataDownloadInput.json
View File

@@ -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 ]
]
]
}
}
]

2
FarmmapsDataDownload/FarmmapsDataDownload.csproj
View File

@@ -2,7 +2,7 @@
<PropertyGroup>
<OutputType>Exe</OutputType>
<TargetFramework>netcoreapp3.0</TargetFramework>
<TargetFramework>netcoreapp3.1</TargetFramework>
</PropertyGroup>
<ItemGroup>

3
FarmmapsDataDownload/Models/DataDownloadInput.cs
View File

@@ -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<string> SatelliteBands { get; set; }
public bool StoreSatelliteStatisticsSingleImage { get; set; }
public bool StoreSatelliteStatisticsCropYear { get; set; }
public bool StoreVanDerSatStatistics { get; set; }

114
FarmmapsDataDownload/ShowGeotiff.r Normal file
View File

@@ -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