Geographic Information Systems Asked by daves_ma on May 6, 2021
I want to interpolate some daily point data. The ultimate goal is to receive something like a heat map/IDW map with some interpolation statistics like kriging/IDW.
My data:
Nr. : coordinates : Sensor_Name : X04.Apr : X15.Apr : X01.Mai : X16.Mai : X01.Jun : X16.Jun : X01.Jul
1 (16.623, 48.21022) Wst_cz3_024 22 21 26 33 30 19 18
2 (16.623, 48.21022) Wst_cz2_041 24 23 34 38 35 21 19
3 (16.623, 48.21022) Wst_cz1_019 27 26 30 40 36 23 21
4 (16.62026, 48.21169) Wst_at3_022 30 28 32 38 36 24 21
5 (16.62026, 48.21169) Wst_at2_008 28 26 30 36 34 22 20
6 (16.62026, 48.21169) Wst_at1_001 30 28 27 36 34 24 17
7 (16.62104, 48.21312) Scint2_3_013 15 13 16 26 18 4 6
8 (16.62104, 48.21312) Scint2_2_032 24 22 23 31 28 18 14
9 (16.623, 48.21022) Scint2_1_029 32 30 35 40 37 24 21
.
.
.
.
136 (16.62036, 48.20904) Ru104_005 20 18 34 42 38 20 17
137 (16.62246, 48.21033) Ru103_070 19 19 33 40 36 22 20
138 (16.621, 48.20857) Ru103_057 22 21 35 42 39 24 22
The current code.
#load data#
kriging<-fread(...)
options(max.print=1000000)
min_LON = min(kriging$LON)
min_LAT = min(kriging$LAT)
LON_length = max(kriging$LON - min_LON)
LAT_length = max(kriging$LAT - min_LAT)
cellsize = 50
ncol = round(LON_length/cellsize, 0)
nrow = round(LAT_length/cellsize, 0)
coordinates(kriging) = ~LON+LAT
proj4string(kriging) = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
head(kriging[,c(1,3,4,5,6,7,8,9)], 138)
Probably the problem is cellsize, ncol and nrow. How would I know how big my cellsize should be?
Unfortunately I am stuck with this last part of code. I cant figure out how to finish my goal.
grid = GridTopology(cellcentre.offset = c(min_LON,min_LAT), cellsize = c(cellsize, cellsize), cells.dim = c(ncol,nrow))
The error I get is:
Error in validObject(.Object) :
invalid class “GridTopology” object: cells.dim has incorrect dimension
grid = SpatialPixelsDataFrame(grid,
data=data.frame(id=1:prod(ncol,nrow)),
proj4string=CRS(proj4string(kriging)))
plot(grid)
##autokriging
kriging_result = autoKrige(X04.Apr~1, kriging, kriging.grid)
plot(kriging_result)
I tried to follow the instructions according to this post:
‘https://gis.stackexchange.com/questions/158021/plotting-map-resulted-from-kriging-in-r/164421’
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP