|
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
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 <bakaryf...@gmail.com
> <mailto:bakaryf...@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 <btup...@bigelow.org
> <mailto:btup...@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
> <bakaryf...@gmail.com
<mailto:bakaryf...@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
>
_______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
<mailto:R-sig-Geo@r-project.org>
> 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
> R-sig-Geo@r-project.org
> 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
--
--
Micha Silver
Ben Gurion Univ
Sde Boker, Remote Sensing Lab
+972-523-665918
|
