Subscribe to R-sig-geo feed
This is an archive for R-sig-geo, not a forum. It is not possible to post through Nabble - you may not start a new thread nor follow up an existing thread. If you wish to post, but are not subscribed to the list through the list homepage, subscribe first (through the list homepage, not via Nabble) and post from your subscribed email address. Until 2015-06-20, subscribers could post through Nabble, but policy has been changed as too many non-subscribers misunderstood the interface.
Updated: 41 min 31 sec ago

Re: how to find projections supported by sf::st_crs()?

Tue, 01/21/2020 - 09:54
Please note that anything I am going to write concerns current sf (on
CRAN), and that things will change pretty strongly in a month or so,
driven by GDAL 3.x and PROJ 6.2.x releases.

st_crs() takes a proj4string or an EPSG number, and tries to resolve
this using GDAL. Although it uses PROJ, GDAL is more strict than PROJ in
accepting proj4strings, and as you noted +wintri is not accepted.

sf_project() was written exactly for the reason that it communicates to
PROJ without going through GDAL, and hence is more flexible; all
proj4strings PROJ accepts are fine as input, but those that are not
accepted by GDAL can only be passed as character strings for the reason
mentioned above.

Hth,

On 1/21/20 3:45 PM, Daniel Kelley wrote:
> I am pondering the use of `sf::sf_project()` for map-projection calculations within the `oce` package. My impression is that `sf::st_crs()` ought to be used to process projection strings, instead of handing strings directly to `sf::sf_project()`.  (I think the idea is that this will lead to the addition of extra information about the ellipse model, etc., and perhaps do some checks for errors in the string.)
>
> However, my tests show that `sf::st_crs()` balks at some projections, such as `wintri`, as shown in the code snippets given near the end of this email.  (Those snippets also show that `sf::sf_project()` handles the wintri projection, as does `rgdal::project()`, so the warning does not mean that `sf` will refuse to do the projection.)
>
> This leads me to three questions, that I'm hoping others can shed light on.
>
> 1. Am I right in thinking that I ought to use `sf::sf_crs()` to "clean up" my projection specifications?
>
> 2. Is there a way to find the list of projections that `sf::sf_crs()` accepts without producing `NA` results and warnings?
>
> 3. Should I avoid using projections that `sf::sf_crs()` warns about?
>
> I apologize if I ought to have gathered the results from my reading.  I could be looking in the wrong places.
>
> Thanks in advance for any advice!  -- Dan.
>
> ```R
>> library(sf)
> Linking to GEOS 3.8.0, GDAL 2.4.2, PROJ 6.2.1
>
>> sf::sf_project("+proj=lonlat", "+proj=wintri", cbind(0,0))
>      [,1] [,2]
> [1,]    0    0
>
>> sf::st_crs("+proj=wintri")
> Coordinate Reference System: NA
> Warning message:
> In CPL_crs_from_proj4string(x) :
>   GDAL cannot import PROJ.4 string `+proj=wintri': returning missing CRS
> ```
>
> Dan E. Kelley [he/him/his 314ppm]
> Dalhousie University
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

pEpkey.asc (3K) Download Attachment

how to find projections supported by sf::st_crs()?

Tue, 01/21/2020 - 08:45
I am pondering the use of `sf::sf_project()` for map-projection calculations within the `oce` package. My impression is that `sf::st_crs()` ought to be used to process projection strings, instead of handing strings directly to `sf::sf_project()`.  (I think the idea is that this will lead to the addition of extra information about the ellipse model, etc., and perhaps do some checks for errors in the string.)

However, my tests show that `sf::st_crs()` balks at some projections, such as `wintri`, as shown in the code snippets given near the end of this email.  (Those snippets also show that `sf::sf_project()` handles the wintri projection, as does `rgdal::project()`, so the warning does not mean that `sf` will refuse to do the projection.)

This leads me to three questions, that I'm hoping others can shed light on.

1. Am I right in thinking that I ought to use `sf::sf_crs()` to "clean up" my projection specifications?

2. Is there a way to find the list of projections that `sf::sf_crs()` accepts without producing `NA` results and warnings?

3. Should I avoid using projections that `sf::sf_crs()` warns about?

I apologize if I ought to have gathered the results from my reading.  I could be looking in the wrong places.

Thanks in advance for any advice!  -- Dan.

```R
> library(sf)
Linking to GEOS 3.8.0, GDAL 2.4.2, PROJ 6.2.1

> sf::sf_project("+proj=lonlat", "+proj=wintri", cbind(0,0))
     [,1] [,2]
[1,]    0    0

> sf::st_crs("+proj=wintri")
Coordinate Reference System: NA
Warning message:
In CPL_crs_from_proj4string(x) :
  GDAL cannot import PROJ.4 string `+proj=wintri': returning missing CRS
```

Dan E. Kelley [he/him/his 314ppm]
Dalhousie University

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Testing --- please ignore.

Mon, 01/20/2020 - 03:40

Sorry for the noise. I'm just trying to test a message filter on my mailer.

cheers,

Rolf Turner

--
Honorary Research Fellow
Department of Statistics
University of Auckland
Phone: +64-9-373-7599 ext. 88276

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: write to netcdf

Fri, 01/17/2020 - 12:32
Hi,

Have you tried, for example, the methods in the RNetCDF package described here?
https://boris.unibe.ch/47220/7/Michna_Woods.pdf

You are likely to get more useful (and faster) responses if you tell
us what you've tried, and what specifically didn't work.

Sarah

On Tue, Jan 14, 2020 at 9:23 AM Abdoulaye Sarr <[hidden email]> wrote:
>
> I have a datasets with below profile I want to write/export to netcdf :
>
> .@ dates      : Date[1:2], format: "2019-03-15" "2019-07-30"
>   ..@ data       :List of 2
>   .. ..$ 2019-03-15:'data.frame': 47005 obs. of  11 variables:
>   .. .. ..$ MeanTemperature     : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ MinTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ MaxTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ Precipitation       : num [1:47005] 0 0 0 0 0 0 0 0 0 0 ...
>
>   .. ..$ 2019-07-30:'data.frame': 47005 obs. of  11 variables:
>   .. .. ..$ MeanTemperature     : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ MinTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ MaxTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>   .. .. ..$ Precipitation       : num [1:47005] NA NA NA NA NA NA NA NA NA
> NA ...
>
>   ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
>   .. .. ..@ cellcentre.offset: Named num [1:2] 289868 1432838
>   .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
>   .. .. ..@ cellsize         : num [1:2] 902 922
>   .. .. ..@ cells.dim        : int [1:2] 395 119
>   ..@ bbox       : num [1:2, 1:2] 289417 1432377 645707 1542095
>   .. ..- attr(*, "dimnames")=List of 2
>   .. .. ..$ : chr [1:2] "s1" "s2"
>   .. .. ..$ : chr [1:2] "min" "max"
>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
>   .. .. ..@ projargs: chr "+init=epsg:32628 +proj=utm +zone=28 +datum=WGS84
> +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
>
> Any hint on this much appreciated.
>
> AS
--
Sarah Goslee (she/her)
http://www.numberwright.com

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: indexing multi-layer rasters

Fri, 01/17/2020 - 10:07
Ahhhh, the lightbulb just went off!  That makes perfect sense.

Thanks!
Ben

On Fri, Jan 17, 2020 at 10:45 AM Bede-Fazekas Ákos <[hidden email]> wrote:
>
> Hello Ben,
>
> I think you are absolutely right about "raster's implementation of `[[`
> is different than base R". But I don't agree on your interpretation on
> the logical indexing of rasters ('first logical index is used'). It is
> not the first one. The logical vector is coerced to integer, and c(TRUE,
> TRUE, TRUE) is treated as c(1, 1, 1), while c(TRUE, FALSE, TRUE) is
> treated as c(1, 0, 1). This is why it cause a 'not a valid subset' error
> (there is no sense of searching for the 0th layer of the RasterStack).
> If the first logical index was used, c(TRUE, TRUE, TRUE) and c(TRUE,
> FALSE, TRUE) would give exactly the same result, i.e a RasterLayer
> created from the first layer of the RasterStack.
>
> Have a nice weekend,
> Ákos
>
>
> 2020.01.17. 15:40 keltezéssel, Ben Tupper írta:
> > Hi,
> >
> > Thanks for this.  I think the root of my muddle is the mish-mash model
> > of a raster that I have in my head.  Depending upon what is most
> > convenient, I sometimes view a raster as a multi-dimensional array and
> > at other times as a list of single layer matrices.   If we step back
> > from logical indexing and use integers it is easier to identify
> > raster's slight variation on base R recursive indexing.  The example
> > below uses integer indices on a list and on a raster.
> >
> > Back to logical indexing, in a way it makes perfect sense as just the
> > first logical index is used; just as if() does.  But what is different
> > is that it uses that first logical index for each element in the index
> > vector.  That's why logo[[c(TRUE, TRUE, TRUE)]] yields red.1, red.2
> > and red.3.
> >
> > Thanks again and cheers,
> > Ben
> >
> >
> > library(raster)
> > x <- list(red = "R", green = "G", blue = "B")
> > logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))
> >
> > # `[` index a list, get a list back (with NULL for not found)
> > x[c(1,3)]
> > # $red
> > # [1] "R"
> > #
> > # $green
> > # [1] "G"
> >
> > # `[` index a raster, get a matrix back (or vector for single layer)
> > logo[c(1,3)]
> > #      red green blue
> > # [1,] 255   255  255
> > # [2,] 255   255  255
> >
> > # `[[` recursive index of a list fails
> > x[[c(1,3)]]
> > # Error in x[[c(1, 3)]] : subscript out of bounds
> >
> > # `[[` index of a raster yields raster
> > # so raster's implementation of `[[` is different than base R
> > logo[[c(1,3)]]
> > # class      : RasterStack
> > # dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
> > # resolution : 1, 1  (x, y)
> > # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> > # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > # names      : red, blue
> > # min values :   0,    0
> > # max values : 255,  255
> >
> > On Fri, Jan 17, 2020 at 7:38 AM Bede-Fazekas Ákos <[hidden email]> wrote:
> >>
> >> Dear Ben,
> >> Although I cannot answer your question on why logical subsetting was not
> >> implemented in package raster, there is a very easy workaround:
> >> logo[[(1:nlayers(logo))[c(TRUE, TRUE, TRUE)]]]
> >> logo[[(1:nlayers(logo))[c(TRUE, FALSE, TRUE)]]]
> >>
> >> Also note that in case of lists '[[' does recursive indexing, and this
> >> type of logical indexing you are asking about works only with '['.
> >> HTH,
> >> Ákos Bede-Fazekas
> >> Hungarian Academy of Sciences
> >>
> >> 2020.01.16. 17:50 keltezéssel, Ben Tupper írta:
> >>> Hi,
> >>>
> >>> I usually avoid using logical indexes with multilayer rasters (stacks
> >>> and bricks), but I have never understood why indexing ala `[[` with
> >>> logicals isn't supported.  Below is an example that shows how it
> >>> doesn't work like the typical indexing with logicals.  I'm not asking
> >>> to have logical indices considered (although it would be handy), but
> >>> instead I really just want to understand why it is the way it is.  I
> >>> scanned over the introductory vignette (https://rspatial.org/raster)
> >>> and issues (https://github.com/rspatial/raster/issues) but found
> >>> nothing there.
> >>>
> >>> Thanks and cheers,
> >>> Ben
> >>>
> >>> ### START
> >>> library(raster)
> >>>
> >>> logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))
> >>>
> >>> # red-red-red
> >>> logo[[c(TRUE, TRUE, TRUE)]]
> >>> # class      : RasterStack
> >>> # dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
> >>> # resolution : 1, 1  (x, y)
> >>> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> >>> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >>> # names      : red.1, red.2, red.3
> >>> # min values :     0,     0,     0
> >>> # max values :   255,   255,   255
> >>>
> >>> # red-red
> >>> logo[[c(TRUE, TRUE)]]
> >>> # class      : RasterStack
> >>> # dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
> >>> # resolution : 1, 1  (x, y)
> >>> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> >>> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >>> # names      : red.1, red.2
> >>> # min values :     0,     0
> >>> # max values :   255,   255
> >>>
> >>> # red
> >>> logo[[c(TRUE)]]
> >>> # class      : RasterLayer
> >>> # band       : 1  (of  3  bands)
> >>> # dimensions : 77, 101, 7777  (nrow, ncol, ncell)
> >>> # resolution : 1, 1  (x, y)
> >>> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> >>> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> >>> # source     : /usr/lib64/R/library/raster/external/rlogo.grd
> >>> # names      : red
> >>> # values     : 0, 255  (min, max)
> >>>
> >>> logo[[c(TRUE, FALSE, TRUE)]]
> >>> #Error in .local(x, ...) : not a valid subset
> >>>
> >>>
> >>> #sessionInfo()
> >>> # R version 3.5.1 (2018-07-02)
> >>> # Platform: x86_64-redhat-linux-gnu (64-bit)
> >>> # Running under: CentOS Linux 7 (Core)
> >>> #
> >>> # Matrix products: default
> >>> # BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
> >>> #
> >>> # locale:
> >>> #   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> >>> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
> >>> # [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
> >>> LC_PAPER=en_US.UTF-8       LC_NAME=C
> >>> # [9] LC_ADDRESS=C               LC_TELEPHONE=C
> >>> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> >>> #
> >>> # attached base packages:
> >>> #   [1] stats     graphics  grDevices utils     datasets  methods   base
> >>> #
> >>> # other attached packages:
> >>> #   [1] raster_3.0-7 sp_1.3-2
> >>> #
> >>> # loaded via a namespace (and not attached):
> >>> #   [1] compiler_3.5.1   rgdal_1.4-8      tools_3.5.1      yaml_2.2.0
> >>>        Rcpp_1.0.3       codetools_0.2-15
> >>> # [7] grid_3.5.1       lattice_0.20-35
> >>>
> >>> ### END
> >>>
> >>>
> >>>
> >>>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> [hidden email]
> >> 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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: indexing multi-layer rasters

Fri, 01/17/2020 - 09:12
Hmmm... maybe the below might be easier to remember (at least it is for me)
then:

logo[[which(c(TRUE, TRUE, TRUE))]]
logo[[which(c(TRUE, FALSE, TRUE))]]

Thanks Akos!

On Fri, Jan 17, 2020 at 7:38 AM Bede-Fazekas Ákos <[hidden email]>
wrote:

>
> Dear Ben,
> Although I cannot answer your question on why logical subsetting was not
> implemented in package raster, there is a very easy workaround:
> logo[[(1:nlayers(logo))[c(TRUE, TRUE, TRUE)]]]
> logo[[(1:nlayers(logo))[c(TRUE, FALSE, TRUE)]]]
>
> Also note that in case of lists '[[' does recursive indexing, and this
> type of logical indexing you are asking about works only with '['.
> HTH,
> Ákos Bede-Fazekas
> Hungarian Academy of Sciences
>
> 2020.01.16. 17:50 keltezéssel, Ben Tupper írta:
> > Hi,
> >
> > I usually avoid using logical indexes with multilayer rasters (stacks
> > and bricks), but I have never understood why indexing ala `[[` with
> > logicals isn't supported.  Below is an example that shows how it
> > doesn't work like the typical indexing with logicals.  I'm not asking
> > to have logical indices considered (although it would be handy), but
> > instead I really just want to understand why it is the way it is.  I
> > scanned over the introductory vignette (https://rspatial.org/raster)
> > and issues (https://github.com/rspatial/raster/issues) but found
> > nothing there.
> >
> > Thanks and cheers,
> > Ben
> >
> > ### START
> > library(raster)
> >
> > logo <- raster::brick(system.file("external/rlogo.grd",
> package="raster"))
> >
> > # red-red-red
> > logo[[c(TRUE, TRUE, TRUE)]]
> > # class      : RasterStack
> > # dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
> > # resolution : 1, 1  (x, y)
> > # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> > # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > # names      : red.1, red.2, red.3
> > # min values :     0,     0,     0
> > # max values :   255,   255,   255
> >
> > # red-red
> > logo[[c(TRUE, TRUE)]]
> > # class      : RasterStack
> > # dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
> > # resolution : 1, 1  (x, y)
> > # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> > # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > # names      : red.1, red.2
> > # min values :     0,     0
> > # max values :   255,   255
> >
> > # red
> > logo[[c(TRUE)]]
> > # class      : RasterLayer
> > # band       : 1  (of  3  bands)
> > # dimensions : 77, 101, 7777  (nrow, ncol, ncell)
> > # resolution : 1, 1  (x, y)
> > # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> > # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> > # source     : /usr/lib64/R/library/raster/external/rlogo.grd
> > # names      : red
> > # values     : 0, 255  (min, max)
> >
> > logo[[c(TRUE, FALSE, TRUE)]]
> > #Error in .local(x, ...) : not a valid subset
> >
> >
> > #sessionInfo()
> > # R version 3.5.1 (2018-07-02)
> > # Platform: x86_64-redhat-linux-gnu (64-bit)
> > # Running under: CentOS Linux 7 (Core)
> > #
> > # Matrix products: default
> > # BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
> > #
> > # locale:
> > #   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> > LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
> > # [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
> > LC_PAPER=en_US.UTF-8       LC_NAME=C
> > # [9] LC_ADDRESS=C               LC_TELEPHONE=C
> > LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> > #
> > # attached base packages:
> > #   [1] stats     graphics  grDevices utils     datasets  methods   base
> > #
> > # other attached packages:
> > #   [1] raster_3.0-7 sp_1.3-2
> > #
> > # loaded via a namespace (and not attached):
> > #   [1] compiler_3.5.1   rgdal_1.4-8      tools_3.5.1      yaml_2.2.0
> >       Rcpp_1.0.3       codetools_0.2-15
> > # [7] grid_3.5.1       lattice_0.20-35
> >
> > ### END
> >
> >
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

--
Vijay Lulla

Assistant Professor,
Dept. of Geography, IUPUI
425 University Blvd, CA-207C.
Indianapolis, IN-46202
[hidden email]
ORCID: https://orcid.org/0000-0002-0823-2522
Online: http://vijaylulla.com | https://github.com/vlulla

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: indexing multi-layer rasters

Thu, 01/16/2020 - 11:53

Dear Ben,
Although I cannot answer your question on why logical subsetting was not
implemented in package raster, there is a very easy workaround:
logo[[(1:nlayers(logo))[c(TRUE, TRUE, TRUE)]]]
logo[[(1:nlayers(logo))[c(TRUE, FALSE, TRUE)]]]

Also note that in case of lists '[[' does recursive indexing, and this
type of logical indexing you are asking about works only with '['.
HTH,
Ákos Bede-Fazekas
Hungarian Academy of Sciences

2020.01.16. 17:50 keltezéssel, Ben Tupper írta:
> Hi,
>
> I usually avoid using logical indexes with multilayer rasters (stacks
> and bricks), but I have never understood why indexing ala `[[` with
> logicals isn't supported.  Below is an example that shows how it
> doesn't work like the typical indexing with logicals.  I'm not asking
> to have logical indices considered (although it would be handy), but
> instead I really just want to understand why it is the way it is.  I
> scanned over the introductory vignette (https://rspatial.org/raster)
> and issues (https://github.com/rspatial/raster/issues) but found
> nothing there.
>
> Thanks and cheers,
> Ben
>
> ### START
> library(raster)
>
> logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))
>
> # red-red-red
> logo[[c(TRUE, TRUE, TRUE)]]
> # class      : RasterStack
> # dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
> # resolution : 1, 1  (x, y)
> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> # names      : red.1, red.2, red.3
> # min values :     0,     0,     0
> # max values :   255,   255,   255
>
> # red-red
> logo[[c(TRUE, TRUE)]]
> # class      : RasterStack
> # dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
> # resolution : 1, 1  (x, y)
> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> # names      : red.1, red.2
> # min values :     0,     0
> # max values :   255,   255
>
> # red
> logo[[c(TRUE)]]
> # class      : RasterLayer
> # band       : 1  (of  3  bands)
> # dimensions : 77, 101, 7777  (nrow, ncol, ncell)
> # resolution : 1, 1  (x, y)
> # extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
> # crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> # source     : /usr/lib64/R/library/raster/external/rlogo.grd
> # names      : red
> # values     : 0, 255  (min, max)
>
> logo[[c(TRUE, FALSE, TRUE)]]
> #Error in .local(x, ...) : not a valid subset
>
>
> #sessionInfo()
> # R version 3.5.1 (2018-07-02)
> # Platform: x86_64-redhat-linux-gnu (64-bit)
> # Running under: CentOS Linux 7 (Core)
> #
> # Matrix products: default
> # BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
> #
> # locale:
> #   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
> # [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
> LC_PAPER=en_US.UTF-8       LC_NAME=C
> # [9] LC_ADDRESS=C               LC_TELEPHONE=C
> LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
> #
> # attached base packages:
> #   [1] stats     graphics  grDevices utils     datasets  methods   base
> #
> # other attached packages:
> #   [1] raster_3.0-7 sp_1.3-2
> #
> # loaded via a namespace (and not attached):
> #   [1] compiler_3.5.1   rgdal_1.4-8      tools_3.5.1      yaml_2.2.0
>       Rcpp_1.0.3       codetools_0.2-15
> # [7] grid_3.5.1       lattice_0.20-35
>
> ### END
>
>
>
>
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

indexing multi-layer rasters

Thu, 01/16/2020 - 10:50
Hi,

I usually avoid using logical indexes with multilayer rasters (stacks
and bricks), but I have never understood why indexing ala `[[` with
logicals isn't supported.  Below is an example that shows how it
doesn't work like the typical indexing with logicals.  I'm not asking
to have logical indices considered (although it would be handy), but
instead I really just want to understand why it is the way it is.  I
scanned over the introductory vignette (https://rspatial.org/raster)
and issues (https://github.com/rspatial/raster/issues) but found
nothing there.

Thanks and cheers,
Ben

### START
library(raster)

logo <- raster::brick(system.file("external/rlogo.grd", package="raster"))

# red-red-red
logo[[c(TRUE, TRUE, TRUE)]]
# class      : RasterStack
# dimensions : 77, 101, 7777, 3  (nrow, ncol, ncell, nlayers)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names      : red.1, red.2, red.3
# min values :     0,     0,     0
# max values :   255,   255,   255

# red-red
logo[[c(TRUE, TRUE)]]
# class      : RasterStack
# dimensions : 77, 101, 7777, 2  (nrow, ncol, ncell, nlayers)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# names      : red.1, red.2
# min values :     0,     0
# max values :   255,   255

# red
logo[[c(TRUE)]]
# class      : RasterLayer
# band       : 1  (of  3  bands)
# dimensions : 77, 101, 7777  (nrow, ncol, ncell)
# resolution : 1, 1  (x, y)
# extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
# crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
# source     : /usr/lib64/R/library/raster/external/rlogo.grd
# names      : red
# values     : 0, 255  (min, max)

logo[[c(TRUE, FALSE, TRUE)]]
#Error in .local(x, ...) : not a valid subset


#sessionInfo()
# R version 3.5.1 (2018-07-02)
# Platform: x86_64-redhat-linux-gnu (64-bit)
# Running under: CentOS Linux 7 (Core)
#
# Matrix products: default
# BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so
#
# locale:
#   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
LC_PAPER=en_US.UTF-8       LC_NAME=C
# [9] LC_ADDRESS=C               LC_TELEPHONE=C
LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base
#
# other attached packages:
#   [1] raster_3.0-7 sp_1.3-2
#
# loaded via a namespace (and not attached):
#   [1] compiler_3.5.1   rgdal_1.4-8      tools_3.5.1      yaml_2.2.0
     Rcpp_1.0.3       codetools_0.2-15
# [7] grid_3.5.1       lattice_0.20-35

### END




--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

write to netcdf

Tue, 01/14/2020 - 08:23
I have a datasets with below profile I want to write/export to netcdf :

.@ dates      : Date[1:2], format: "2019-03-15" "2019-07-30"
  ..@ data       :List of 2
  .. ..$ 2019-03-15:'data.frame': 47005 obs. of  11 variables:
  .. .. ..$ MeanTemperature     : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ MinTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ MaxTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ Precipitation       : num [1:47005] 0 0 0 0 0 0 0 0 0 0 ...

  .. ..$ 2019-07-30:'data.frame': 47005 obs. of  11 variables:
  .. .. ..$ MeanTemperature     : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ MinTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ MaxTemperature      : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...
  .. .. ..$ Precipitation       : num [1:47005] NA NA NA NA NA NA NA NA NA
NA ...

  ..@ grid       :Formal class 'GridTopology' [package "sp"] with 3 slots
  .. .. ..@ cellcentre.offset: Named num [1:2] 289868 1432838
  .. .. .. ..- attr(*, "names")= chr [1:2] "s1" "s2"
  .. .. ..@ cellsize         : num [1:2] 902 922
  .. .. ..@ cells.dim        : int [1:2] 395 119
  ..@ bbox       : num [1:2, 1:2] 289417 1432377 645707 1542095
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "s1" "s2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+init=epsg:32628 +proj=utm +zone=28 +datum=WGS84
+units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"

Any hint on this much appreciated.

AS

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: setting spatial units to miles in variogramST

Mon, 01/13/2020 - 09:33
That’s perfect.

Thanks,
Erin


On Mon, Jan 13, 2020 at 3:14 AM Dr. Benedikt Gräler <[hidden email]>
wrote:

> Dear Erin,
>
> the spatial unit is determined by the CRS of the sp slot of the STIDF
> object that you provide to variogramST; details are e.g. in ?spDists.
>
> HTH,
>
>   Ben
>
>
> On 12.01.20 02:39, Erin Hodgess wrote:
> > Hello!
> >
> > In the variogramST function, I set "tunits=mins".  Is there a way to
> change
> > from km to miles, please?
> >
> > Thanks,
> > Erin
> >
> >
> > Erin Hodgess, PhD
> > mailto: [hidden email]
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Dr. Benedikt Gräler
> 52°North Initiative for Geospatial Open Source Software GmbH
> Martin-Luther-King-Weg 24
> 48155 Muenster, Germany
>
> E-Mail: [hidden email]
> Fon: +49-(0)-251/396371-39
> Fax: +49-(0)-251/396371-11
>
> http://52north.org/
> Twitter: @FiveTwoN
>
> General Managers:
> Prof. Dr. Albert Remke, Prof. Dr. Andreas Wytzisk-Arens
> Local Court Muenster HRB 10849
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Erin Hodgess, PhD
mailto: [hidden email]

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: setting spatial units to miles in variogramST

Mon, 01/13/2020 - 04:14
Dear Erin,

the spatial unit is determined by the CRS of the sp slot of the STIDF
object that you provide to variogramST; details are e.g. in ?spDists.

HTH,

  Ben


On 12.01.20 02:39, Erin Hodgess wrote:
> Hello!
>
> In the variogramST function, I set "tunits=mins".  Is there a way to change
> from km to miles, please?
>
> Thanks,
> Erin
>
>
> Erin Hodgess, PhD
> mailto: [hidden email]
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: [hidden email]
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

http://52north.org/
Twitter: @FiveTwoN

General Managers:
Prof. Dr. Albert Remke, Prof. Dr. Andreas Wytzisk-Arens
Local Court Muenster HRB 10849

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Mon, 01/13/2020 - 03:39
Dear Silver,

Thank for your help.

I think there had been a mix-up in the Coordinate Rference System. When I
plot this shapefile I get
it in the UTM Coordinate.
Le Senegal is in Zone 28 North. The watershed on what I work is shared by
trois Countries (Senegal, Gambie and Guinea).
But in the choose of UTM 33 in my script is an error.
Thank very much.

Best regards


Le lun. 13 janv. 2020 à 07:22, Micha Silver <[hidden email]> a écrit :

> 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 <[hidden email]> 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 <[hidden email]
>> > <mailto:[hidden email]>> 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 <[hidden email]
>> >     <mailto:[hidden email]>> 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
>> >         <[hidden email] <mailto:[hidden email]>> 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
>> >             [hidden email] <mailto:[hidden email]>
>> >             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
>> > [hidden email]
>> > 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
>
> --
> Micha Silver
> Ben Gurion Univ
> Sde Boker, Remote Sensing Lab
> +972-523-665918
>
>
--



Bakary

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Mon, 01/13/2020 - 01:23

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 <[hidden email]> 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 <[hidden email]
> <mailto:[hidden email]>> 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 <[hidden email]
>     <mailto:[hidden email]>> 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
>         <[hidden email] <mailto:[hidden email]>> 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
>             [hidden email] <mailto:[hidden email]>
>             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
> [hidden email]
> 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
-- Micha Silver Ben Gurion Univ Sde Boker, Remote Sensing Lab +972-523-665918
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

do_R_code_for_isohyet_map[R-sig-Geo].R (1K) Download Attachment

setting spatial units to miles in variogramST

Sat, 01/11/2020 - 19:39
Hello!

In the variogramST function, I set "tunits=mins".  Is there a way to change
from km to miles, please?

Thanks,
Erin


Erin Hodgess, PhD
mailto: [hidden email]

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Fri, 01/10/2020 - 03:49
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 <[hidden email]> 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 <[hidden email]
> > <mailto:[hidden email]>> 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 <[hidden email]
> >     <mailto:[hidden email]>> 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
> >         <[hidden email] <mailto:[hidden email]>> 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
> >             [hidden email] <mailto:[hidden email]>
> >             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
> > [hidden email]
> > 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Fri, 01/10/2020 - 03:41

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 <[hidden email]
> <mailto:[hidden email]>> 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 <[hidden email]
>     <mailto:[hidden email]>> 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
>         <[hidden email] <mailto:[hidden email]>> 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
>             [hidden email] <mailto:[hidden email]>
>             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
> [hidden email]
> 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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Thu, 01/09/2020 - 13:20
On Thu, 9 Jan 2020, Ben Tupper wrote:

> Your assignments that look like...
>
> minx <-  rain_data_UTM at bbox[1,1]
>
> are not valid R statements - and that will cause an error.

I think that this might be the mail or some other system replacing
commercial at by at. Better syntax might be:

minx <-  slot(rain_data_UTM, "bbox")[1,1]

if this was the real intention, as slot() is the preferred script-level S4
slot accessor, although one finds lots of untidy examples using
commercial at in blogs and elsewhere on the net. However, your version is
obviously much easier to read:

> Instead,
> obtain a matrix of the bounding box using the bbox() function.  Then
> extract your coordinates from that. I think you want...
>
> bb <- bbox(rain_data_UTM)
> minx <- bb[1,1]
> maxx <- bb[1,2]
> miny <- bb[2,1]
> maxy <- bb[2,2]

And yes, this is very good advice!

Roger

>
> If you haven't seen this, https://www.rspatial.org/ , it is well worth your
> time.  There are a lot of other great resources about using spatial data in
> R.  Try searching the Rseek.org. Like this
> https://rseek.org/?q=spatial+tutorials  It is a gold mine.
>
>
> On Thu, Jan 9, 2020 at 1:42 PM Bakary Faty <[hidden email]> 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, my data file and my sahefile of
>> watershed.
>>
>> Best regards
>>
>>
>> Le jeu. 9 janv. 2020 à 17:47, Ben Tupper <[hidden email]> 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 <[hidden email]> 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
>>>> [hidden email]
>>>> 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
>>
>
>
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

Re: Help

Thu, 01/09/2020 - 13:13
Thanks very much angain,
 I'll check them out.

Best regards

Le jeu. 9 janv. 2020 à 19:05, Ben Tupper <[hidden email]> a écrit :

> Your assignments that look like...
>
> minx <-  rain_data_UTM at bbox[1,1]
>
> are not valid R statements - and that will cause an error.  Instead,
> obtain a matrix of the bounding box using the bbox() function.  Then
> extract your coordinates from that. I think you want...
>
> bb <- bbox(rain_data_UTM)
> minx <- bb[1,1]
> maxx <- bb[1,2]
> miny <- bb[2,1]
> maxy <- bb[2,2]
>
> If you haven't seen this, https://www.rspatial.org/ , it is well worth
> your time.  There are a lot of other great resources about using spatial
> data in R.  Try searching the Rseek.org. Like this
> https://rseek.org/?q=spatial+tutorials  It is a gold mine.
>
>
> On Thu, Jan 9, 2020 at 1:42 PM Bakary Faty <[hidden email]> 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, my data file and my sahefile of
>> watershed.
>>
>> Best regards
>>
>>
>> Le jeu. 9 janv. 2020 à 17:47, Ben Tupper <[hidden email]> 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 <[hidden email]>
>>> 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
>>>> [hidden email]
>>>> 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
>>
>
>
> --
> Ben Tupper
> Bigelow Laboratory for Ocean Science
> West Boothbay Harbor, Maine
> http://www.bigelow.org/
> https://eco.bigelow.org
>
>
--



Bakary

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Thu, 01/09/2020 - 13:05
Your assignments that look like...

minx <-  rain_data_UTM at bbox[1,1]

are not valid R statements - and that will cause an error.  Instead,
obtain a matrix of the bounding box using the bbox() function.  Then
extract your coordinates from that. I think you want...

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

If you haven't seen this, https://www.rspatial.org/ , it is well worth your
time.  There are a lot of other great resources about using spatial data in
R.  Try searching the Rseek.org. Like this
https://rseek.org/?q=spatial+tutorials  It is a gold mine.


On Thu, Jan 9, 2020 at 1:42 PM Bakary Faty <[hidden email]> 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, my data file and my sahefile of
> watershed.
>
> Best regards
>
>
> Le jeu. 9 janv. 2020 à 17:47, Ben Tupper <[hidden email]> 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 <[hidden email]> 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
>>> [hidden email]
>>> 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
>

--
Ben Tupper
Bigelow Laboratory for Ocean Science
West Boothbay Harbor, Maine
http://www.bigelow.org/
https://eco.bigelow.org

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: Help

Thu, 01/09/2020 - 13:02
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 <[hidden email]> 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 <[hidden email]> 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 <- 20maxx <- 25miny <- 31maxy <- 36pixel <- 1grd <- 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 <[hidden email]> 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
[hidden email]
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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

DATAFILE_FOR_ISOHYET.csv (1K) Download Attachment
R_code_for_rainfall_isohyet_map.R (2K) Download Attachment

Pages