forked from FarmMaps/FarmMapsApiClient
153 lines
7.2 KiB
R
153 lines
7.2 KiB
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 image 8 for date '2020-06-14'
|
|
|
|
library(raster)
|
|
library(sf)
|
|
library(rgdal)
|
|
setwd("C:/workdir/groenmonitor/")
|
|
# FarmmapsDataDownload generates two files:
|
|
fileGeotiff <- "wenr.tif"
|
|
# fileJpg <- "thumbnail.jpg"
|
|
# FarmmapsBulkSatDownload generates many files. Here is what I tried when inputing the same field for a number of years using BulkSatDownloadInput.json
|
|
# 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') #59.91667
|
|
|
|
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.1654627
|
|
# Same as in downloaded file SatelliteDataStatistics_test_satData_wdvi_2020-06-14.csv
|
|
# "mean": 0.16564258790046743
|
|
# So FindSatelliteItem in C:\git\FarmMapsApiClient_WURtest\FarmmapsApi\Services\GeneralService.cs returns a a stacked geotiff in
|
|
# i.e. the cellstats are calculated from the rectangle
|
|
|
|
# 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((
|
|
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') #0.1654627 # Stats from rectangle, RD
|
|
cellStats(r.wenr.wgs84.wdvi,'mean') #0.1656399 # Stats from rectangle, WGS84
|
|
cellStats(r.wenr.rd.wdvi.pol,'mean') #0.1658702 # Stats from raster clipped by polygon, WGS84
|
|
cellStats(r.wenr.wgs84.wdvi.pol,'mean') #0.1658611 # Stats from raster clipped by polygon, WGS84
|
|
# Conclusion: in this example little difference, but can be quite different!
|
|
|
|
# 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
|