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 <[email protected]> a écrit : > > On 09/01/2020 21:02, Bakary Faty wrote: > > 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 > > > > > 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) > > > > Le jeu. 9 janv. 2020 à 18:41, Bakary Faty <[email protected] > > <mailto:[email protected]>> a écrit : > > > > 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, my data file and my sahefile > > of watershed. > > > > Best regards > > > > > > Le jeu. 9 janv. 2020 à 17:47, Ben Tupper <[email protected] > > <mailto:[email protected]>> a écrit : > > > > Welcome to r-sig-geo! > > > > I don't think that you haven't provided us enough information > > so that we can help. On the other hand, does the example > > below using expand.grid help? > > > > minx <- 20 > > maxx <- 25 > > miny <- 31 > > maxy <- 36 > > pixel <- 1 > > grd <- expand.grid(x = seq(minx, maxx, by=pixel), y = > > seq(miny, maxy, by=pixel)) > > > > Ben > > > > On Thu, Jan 9, 2020 at 11:41 AM Bakary Faty > > <[email protected] <mailto:[email protected]>> wrote: > > > > > > Dear, > > > > I'm writing to express my wish to join R-sig-geo list users. > > I was doing a search on the net to know how to build an > > isohyet map and I came across this R code. > > However, I stumbled upon a problem from the line : > > grd <- expand.grid(x=seq(minx, maxx, by=pixel), > > y=seq(miny, maxy, by=pixel)), > > I get the following error message: > > default method not implemented for type 'S4'. I want to > > know how resolve this error. > > > > Also, I would like to ask you only at the line level: > > minx <- rain_data_UTM at bbox[1,1] > > maxx <- rain_data_UTM at bbox[1,2] > > miny <- rain_data_UTM at bbox[2,1] > > maxy <- rain_data_UTM at bbox[2,2], > > if I should limit myself to "rain_data_UTM" or write > > completely: > > rain_data_UTM at bbox[,]. > > By the way, this is the pointfile I reconstructed. > > You can find it attached to the mail. > > > > Thanks in advance > > > > Best regards > > > > > > > > Bakary > > _______________________________________________ > > R-sig-Geo mailing list > > [email protected] <mailto:[email protected]> > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > > > > > -- > > Ben Tupper > > Bigelow Laboratory for Ocean Science > > West Boothbay Harbor, Maine > > http://www.bigelow.org/ > > https://eco.bigelow.org > > > > > > > > -- > > > > > > > > Bakary > > > > > > > > -- > > > > > > > > Bakary > > > > _______________________________________________ > > R-sig-Geo mailing list > > [email protected] > > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > -- > Micha Silver > Ben Gurion Univ. > Sde Boker, Remote Sensing Lab > cell: +972-523-665918 > https://orcid.org/0000-0002-1128-1325 > > -- Bakary [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
