Skip to content

Help

3 messages · Micha Silver, Bakary Faty

#
Thank you for appreciated reply,

I explane you exactly what I want to do with this R code attached.
I want to adapt this code to my data to build an isohyet map.
But i have some difficulties to adapt it to my case.
I will be very happy when you will help my to adapt this R code (attached)
to my case.
You can find attached the you R code and my data file.
Best regards

Le jeu. 9 janv. 2020 ? 18:41, Bakary Faty <bakaryfaty at gmail.com> a ?crit :

  
    
#
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)

  
    
#
Thank you very much for help,
The code works fine but when I add the shapefile of my watershed
I can not plot it.

Le ven. 10 janv. 2020 ? 09:41, Micha Silver <tsvibar at gmail.com> a ?crit :