Mercurial > repos > ecology > interpolation_run_idw_interpolation

comparison run_idw_interpolation.R @ 0:d07fcc660f3c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/interpolation commit 450e4496f243d6e94d5238358873bbc014fe2f08
author ecology
date Mon, 08 Jan 2024 10:32:25 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:d07fcc660f3c
1 library("getopt")
2 library("sf")
3 library("tmap")
4 library("RColorBrewer")
5 library("raster")
6 library("gstat")
7
8 args = commandArgs(trailingOnly=TRUE)
9 option_specification = matrix(c(
10 'observationsCsv', 'i1', 1, 'character',
11 'latitudeColumn', 'i2', 1, 'double',
12 'longitudeColumn', 'i3', 1, 'double',
13 'observationsColumn', 'i4', 1, 'double',
14 'studyArea', 'i5', 1, 'character',
15 'idwPower', 'i6', 1, 'double',
16 'samplePoints', 'i7', 1, 'double',
17 'sampleType', 'i8', 1, 'character',
18 'legendLabel', 'i9', 1, 'character',
19 'legendPosition', 'i10', 1, 'character',
20 'numberClasses', 'i11', 1, 'double',
21 'dotSize', 'i12', 1, 'double',
22 'colorType', 'i13', 1, 'character',
23 'testCase', 'i14', 1, 'character',
24 'outputData', 'o', 2, 'character'
25 ), byrow=TRUE, ncol=4);
26 options = getopt(option_specification);
27
28 obsData <- read.csv(file=options$observationsCsv, sep = ',', header = TRUE)
29 latitudeColumn <- options$latitudeColumn
30 longitudeColumn <- options$longitudeColumn
31 observationsColumn <- options$observationsColumn
32 studyArea <- options$studyArea
33 idwPower <- options$idwPower
34 samplePoints <- options$samplePoints
35 sampleType <- options$sampleType
36 legendLabel <- options$legendLabel
37 legendPosition <- options$legendPosition
38 numberClasses <- options$numberClasses
39 dotSize <- options$dotSize
40 colorType <- options$colorType
41 testCase <- options$testCase
42
43 #cat("\n observationsCsv", options$observationsCsv)
44 cat("\n latitudeColumn", latitudeColumn)
45 cat("\n longitudeColumn", longitudeColumn)
46 cat("\n observationsColumn", observationsColumn)
47 #cat("\n studyArea", studyArea)
48 cat("\n idwPower", idwPower)
49 cat("\n samplePoints", samplePoints)
50 cat("\n sampleType", sampleType)
51 cat("\n legendLabel", legendLabel)
52 cat("\n legendposition", legendPosition)
53 cat("\n numberClasses", numberClasses)
54 cat("\n dotSize", dotSize)
55 cat("\n colorType", colorType)
56 cat("\n testCase", testCase)
57 #cat("\n outputData: ", options$outputData)
58
59 coordinates(obsData) <- c(colnames(obsData)[longitudeColumn], colnames(obsData)[latitudeColumn])
60 sf_obsData <- as_Spatial(st_as_sf(obsData))
61
62 polygon <- as_Spatial(st_read(studyArea))
63 sf_obsData@bbox<-polygon@bbox
64
65 runInterpolation <- function(points, values, interpolation_power, sample_points, sample_type){
66 if (testCase == "true") {
67 cat("\n set seed!")
68 set.seed(123)
69 }
70 grd <- as.data.frame(spsample(points, sample_type, n=sample_points))
71 names(grd) <- c("X", "Y")
72 coordinates(grd) <- c("X", "Y")
73 gridded(grd) <- TRUE
74 fullgrid(grd) <- TRUE
75
76 proj4string(points) <- proj4string(points)
77 proj4string(grd) <- proj4string(points)
78 return(gstat::idw(values ~ 1, points, newdata=grd, idp=interpolation_power))
79 }
80
81 plotInterpolationMap <- function(raster, points, legend_label){
82 plot <- tm_shape(raster) +
83 tm_raster(n=numberClasses,palette = rev(brewer.pal(7, colorType)), auto.palette.mapping = FALSE,
84 title=legend_label) +
85 tm_shape(points) + tm_dots(size=dotSize) +
86 tm_legend(legend.outside=legendPosition)
87 return(plot)
88 }
89
90 sf_obsData.idw <- runInterpolation(sf_obsData, obsData$measurement, idwPower, samplePoints, sampleType)
91
92 raster_object <- raster(sf_obsData.idw)
93 raster_object.mask <- mask(raster_object, polygon)
94
95 idw <- plotInterpolationMap(raster_object.mask, sf_obsData, legendLabel)
96 idw
97
98 png(options$outputData)
99 idw
100 dev.off()