forked from FarmMaps/FarmMapsApiClient
179 lines
8.9 KiB
R
179 lines
8.9 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 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
|