Hello Bakary:
Please keep questions on the list, so that others can benefit from the
conversation.
I see that the Basins shapefile you sent is declared to be in Long/Lat
WGS84 Coordinate Reference System, however the coordinates (the X-Y
locations) look like it is already projected to UTM 33. Is it possible
there is some mixup in the Coordinate Reference System of the Basins
shapefile?
(That would explain why the basins layer does not appear on your plot)
Furthermore, I think that you are working in Senegal/Gambia (from the
coordinates in the DATAFILE). If so, why did you choose to project to UTM
33 in your script? Isn't Senegal in UTM zone 28?
Attached is the corrected script (changed to UTM28). The corrected basin
shapefile is sent privately. (too big for the list)
Regards, Micha
On 12/01/2020 21:29, Bakary Faty wrote:
Dear Micha Silver
I'm writing you this email following your contribution on my isohyet map
code, in 'r-sig-geo'.
I really appreciated your help which helped me a lot. I can draw the
isohyets but I can't draw the contour of my study basin.
You can see where the problem is in the code to make a change again?
You will find the files of the watershed and my code R attached.
Thank you very much by advance
Best regards
Le ven. 10 janv. 2020 ? 09:41, Micha Silver <tsvibar at gmail.com> 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
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 <bakaryfaty at gmail.com
<mailto:bakaryfaty at gmail.com>> 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 <btupper at bigelow.org
<mailto:btupper at bigelow.org>> 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
<bakaryfaty at gmail.com <mailto:bakaryfaty at gmail.com>> 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