Skip to content
Prev 27858 / 29559 Next

Help

On 09/01/2020 21:02, Bakary Faty wrote:
I made two small changes in your code, and it works fine:

 ?* First I used the suggestion (earlier in this thread) to create your
 ?? grid.
 ?* Then, you had an error in your call to autoKrige.
 ?* After getting that right, I created isohyetal lines with 
rasterToContour

Here's my version


library(automap)
library(raster)
library(rgdal)

## READ INPUT FILE
rain_data <- read.csv(file="DATAFILE_FOR_ISOHYET.csv")

point_coords <- rain_data[c("Lon","Lat")]
coordinates(rain_data) <- point_coords
p4str <- CRS("+init=epsg:4326")
proj4string(rain_data) <- p4str

## CONVERTION TO UTM
p4str_UTM <- CRS("+init=epsg:32633")
rain_data_UTM <- spTransform(rain_data, p4str_UTM)


bb <- bbox(rain_data_UTM)
minx <- bb[1,1]
maxx <- bb[1,2]
miny <- bb[2,1]
maxy <- bb[2,2]

## EACH PIXEL WILL BE 1000 METERS
pixel <- 1000
grd <- expand.grid(x=seq(minx, maxx, by=pixel), y=seq(miny, maxy, 
by=pixel))
coordinates(grd) <- ~x+y
gridded(grd) <- TRUE
proj4string(grd) <- p4str_UTM


## KRIGING, USING AUTOKRIGE WHICH CREATES A BEST GUESS VARIOGRAM
# OK_rain <- autoKrige(Rainfall_value ~ 1, rain_data_UTM, grd)
# There is no variable "Rainfall_value" in your data,
# It is called RAIN_DATA
OK_rain = autoKrige(formula = RAIN_DATA ~1,
 ??????????????????? input_data = rain_data_UTM,
 ??????????????????? new_data = grd)

## TRASFORM TO RASTER
rain_rast <- raster(OK_rain$krige_output)

summary(rain_rast)


# Minimumn is 540, max is 1735
# So create isohyetal lines about every 100 mm and plot

isohyets = rasterToContour(rain_rast, nlevels = 12)
plot(rain_rast)
lines(isohyets, add = TRUE)