Re: [R-sig-Geo] as.im in spatstat::ppm not working as before
On 11/05/2024 10:21, Roozbeh Valavi wrote: This might also help: the solution below will not raise the error if x has geographic (unprojected) coordinates, a case for which you should not use spatstat. library(spatstat) library(terra) .raster_to_im <- function(x){ r <- as.data.frame(x, xy = TRUE) im <- spatstat.geom::as.im <http://as.im>(r) return(im) } # to generalise the spatstat.geom::as.im <http://as.im> to SpatRaster and RasterLayer as.im.SpatRaster <- function(x){ .raster_to_im(x) } as.im.RasterLayer <- function(x){ .raster_to_im(x) } I haven’t used it in a while so there might be a direct support in the spatstat.geom now! Roozbeh Valavi, Biodiversity Modeller & Spatial Ecologist Senior Research Scientist | CSIRO Environment, VIC, Australia *Twitter*:@ValaviRoozbeh | ResearchGate <https://www.researchgate.net/profile/Roozbeh-Valavi> | Google Scholar <https://scholar.google.co.uk/citations?user=3m0jCHwJ=en=ao> On Wed, 8 May 2024 at 7:41 PM, Roger Bivand <mailto:roger.biv...@nhh.no>> wrote: as.im <http://as.im>() was in maptools, which was retired last year. If you need to recreate old work, install maptools from source from the CRAN package archive https://cran.r-project.org/src/contrib/Archive/maptools/ <https://cran.r-project.org/src/contrib/Archive/maptools/>. For non-archival work, update the workflow. See also https://github.com/r-spatial/evolution/issues/8 <https://github.com/r-spatial/evolution/issues/8>. Hope this clarifies, Roger --- Roger Bivand Emeritus Professor Department of Economics Norwegian School of Economics, Bergen, Norway Fra: R-sig-Geo mailto:r-sig-geo-boun...@r-project.org>> på vegne av Edzer Pebesma mailto:edzer.pebe...@uni-muenster.de>> Sendt: onsdag, mai 8, 2024 9:28:54 a.m. Til: r-sig-geo@r-project.org <mailto:r-sig-geo@r-project.org> mailto:r-sig-geo@r-project.org>> Emne: Re: [R-sig-Geo] as.im <http://as.im> in spatstat::ppm not working as before r.grid.world.SGDF.F[[1]] is a numerical vector, not even a matrix, so clearly as.im.default doesn't know what it should do with it. Maybe there was once a method as.im <http://as.im> for r.grid.world.SGDF.F, which is of class SpatialGridDataFrame, but no longer - software evolves. A way you could approach this is converting that object to a stars object: > library(stars) Loading required package: abind Loading required package: sf Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE > as.im <http://as.im>(st_as_stars(r.grid.world.SGDF.F)) Error: Only projected coordinates may be converted to spatstat class objects and that points you to another issue you have: spatstat works with data in Cartesian coordiantes (R2), not with geodetic coordinates (S2). So that's a helpful error message. On 08/05/2024 01:49, Alexandre Santos via R-sig-Geo wrote: > Dear Members, > > I'll calculate some old codes for a paper review, and I'm astonished that as.im <http://as.im>() for fitting the Poisson model is not working as before. In my example: > > library(raster) > library(spatstat) > > > # Elevation > r.grid.world <- raster('https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fraw.githubusercontent.com%2FLeprechault%2Ftrash%2Fmain%2Fwc2.1_10m_elev.tif=05%7C02%7CRoger.Bivand%40nhh.no%7Cd2a2b64634a7462f860708dc6f3082d5%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638507501339842669%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C0%7C%7C%7C=cyducyKF16yIdtfuykXUSifCR8clFtxDhiG1eY4waFU%3D=0 <https://eur02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fraw.githubusercontent.com%2FLeprechault%2Ftrash%2Fmain%2Fwc2.1_10m_elev.tif=05%7C02%7CRoger.Bivand%40nhh.no%7Cd2a2b64634a7462f860708dc6f3082d5%7C33a15b2f849941998d56f20b5aa91af2%7C0%7C0%7C638507501339842669%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C0%7C%7C%7C=cyducyKF16yIdtfuykXUSifCR8clFtxDhiG1eY4waFU%3D=0><https://raw.githubusercontent.com/Leprechault/trash/main/wc2.1_10m_elev.tif <https://raw.githubusercontent.com/Leprechault/trash/main/wc2.1_10m_elev.tif>>') > proj4string(r.grid.world) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") > r.grid.world.SGDF.F <- as(r.grid.world, "SpatialGridDataFrame") > > > #Fitting Poisson model > fit.c <- ppm(nztrees, ~ elev, > covariates=list(elev=as.im <http://as.im>(r.grid.world.SGDF.F[[1]]))) > # > # Error
Re: [R-sig-Geo] as.im in spatstat::ppm not working as before
r.grid.world.SGDF.F[[1]] is a numerical vector, not even a matrix, so clearly as.im.default doesn't know what it should do with it. Maybe there was once a method as.im for r.grid.world.SGDF.F, which is of class SpatialGridDataFrame, but no longer - software evolves. A way you could approach this is converting that object to a stars object: > library(stars) Loading required package: abind Loading required package: sf Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE > as.im(st_as_stars(r.grid.world.SGDF.F)) Error: Only projected coordinates may be converted to spatstat class objects and that points you to another issue you have: spatstat works with data in Cartesian coordiantes (R2), not with geodetic coordinates (S2). So that's a helpful error message. On 08/05/2024 01:49, Alexandre Santos via R-sig-Geo wrote: Dear Members, I'll calculate some old codes for a paper review, and I'm astonished that as.im() for fitting the Poisson model is not working as before. In my example: library(raster) library(spatstat) # Elevation r.grid.world <- raster('https://raw.githubusercontent.com/Leprechault/trash/main/wc2.1_10m_elev.tif') proj4string(r.grid.world) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0") r.grid.world.SGDF.F <- as(r.grid.world, "SpatialGridDataFrame") #Fitting Poisson model fit.c <- ppm(nztrees, ~ elev, covariates=list(elev=as.im(r.grid.world.SGDF.F[[1]]))) # # Error in as.im.default(r.grid.world.SGDF.F[[1]]) : # Can't convert X to a pixel image I tried as matrix and rotate, and the final real-valued pixel image didn't work as a covariate in the model. Please, any help with it? Thanks in advance Alexandre ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma (he/him) Universität Münster, Institute for Geoinformatics Heisenbergstrasse 2, 48149 Münster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Unit of parameters in the expandBB argument
On 28/02/2024 06:04, Xiang Ye via R-sig-Geo wrote: Dear community, I recently revisited the expandBB argument in the plot() function for sf objects. It conveys a numeric vector of length 4 to expand the default canvas (in the order of bottom, left, top, right) when drawing an sf object. A quick refresh is here: https://r.geocompx.org/spatial-class#:~:text=to%20geographic%20data.-,expandBB,-%2C%20for%20example%2C%20can [https://r.geocompx.org/images/cover.png]<https://r.geocompx.org/spatial-class#:~:text=to%20geographic%20data.-,expandBB,-%2C%20for%20example%2C%20can> Chapter 2 Geographic data in R | Geocomputation with R<https://r.geocompx.org/spatial-class#:~:text=to%20geographic%20data.-,expandBB,-%2C%20for%20example%2C%20can> Prerequisites This is the first practical chapter of the book, and therefore it comes with some software requirements. You need access to a computer with a recent version of R installed (R 4.3.2... r.geocompx.org I am curious about the unit for the four values of the expandBB argument. In the help document, it only says the values are "fraction values to expand the bounding box with": https://r-spatial.github.io/sf/reference/plot.html#:~:text=fractional%20values%20to%20expand%20the%20bounding%20box%20with [https://r-spatial.github.io/sf/logo.png]<https://r-spatial.github.io/sf/reference/plot.html#:~:text=fractional%20values%20to%20expand%20the%20bounding%20box%20with> plot sf object �� plot<https://r-spatial.github.io/sf/reference/plot.html#:~:text=fractional%20values%20to%20expand%20the%20bounding%20box%20with> plot one or more attributes of an sf object on a map Plot sf object r-spatial.github.io However, it did not explicitly mention the unit adopted for the values. My guess: For the 1st and 3rd values in the argument, the unit is the height of the bounding box of the sf object; therefore a value of 0.2 expands (downward and upward, respectively) the canvas by 20% of the height of the bounding box. For the 2nd and 4th values in the argument, the unit is the width of the bounding box of the sf object; therefore a value of 0.2 expands (leftward and rightward, respectively) the canvas by 20% of the width of the bounding box. Is my understanding correct? Yes. It takes 4 values for each of the sides (bottom, left, top, right). You can use xlim and ylim to set the plot region in spatial coordinates. Thank you in advance! Ҷ�� YE, Xiang THINKING SPATIALLY<http://www.linkedin.com/in/spatialyexiang>. Ph.D. in Spatial Statistics [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma (he/him) Universität Münster, Institute for Geoinformatics Heisenbergstrasse 2, 48149 Münster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Spatial Data Science across Languages, Muenster Sept 18+19
Dear all, We are organising a small workshop on “Spatial Data Science across Languages” on Sept 18+19 in Muenster, Germany. The goal is to bring developers and users from R-Spatial, GeoPython and GeoJulia together to discuss shared challenges and interests. The website is https://r-spatial.org/sdsl/ . If you are interested in participating please register through the form found on the website before Aug 7 UTC. -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Lo Cape projections in sf and terra
rm: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19044) Matrix products: default locale: [1] LC_COLLATE=English_South Africa.utf8 LC_CTYPE=English_South Africa.utf8 [3] LC_MONETARY=English_South Africa.utf8 LC_NUMERIC=C [5] LC_TIME=English_South Africa.utf8 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] terra_1.6-47 sf_1.0-9 loaded via a namespace (and not attached): [1] Rcpp_1.0.9 magrittr_2.0.3 units_0.8-1tidyselect_1.2.0 [5] R6_2.5.1 rlang_1.0.6fansi_1.0.3dplyr_1.0.10 [9] tools_4.2.2grid_4.2.2 KernSmooth_2.23-20 utf8_1.2.2 [13] cli_3.5.0 e1071_1.7-12 DBI_1.1.3 class_7.3-20 [17] assertthat_0.2.1 tibble_3.1.8 lifecycle_1.0.3codetools_0.2-18 [21] vctrs_0.5.1glue_1.6.2 proxy_0.4-27 compiler_4.2.2 [25] pillar_1.8.1 generics_0.1.3 classInt_0.4-8 pkgconfig_2.0.3 Thank you Dr David Furniss, Wits University. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
Yes, the pipeline approach bypasses GDAL, and doesn't result in an object with an appropriate CRS as a consequence. On 31/01/2023 14:47, Andrea Gilardi wrote: Thank you very much for your example! I briefly checked the examples reported in ?st_transform and the docs at the PROJ website around "Unit conversion" (https://proj.org/operations/conversions/unitconvert.html) and I found that the following code also seems to work (although with a warning message that I'm not sure I understand): library(sp) demo(meuse, echo = FALSE, ask = FALSE) library(sf) #> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE sf_proj_network(TRUE) #> [1] "https://cdn.proj.org; st_as_sf(meuse) |> st_bbox() #> xmin ymin xmax ymax #> 178605 329714 181390 333611 meuse2 <- st_as_sf(meuse) |> st_transform(pipeline = "+proj=pipeline +step +proj=unitconvert +xy_in=m +xy_out=km") #> Warning in st_transform.sfc(st_geometry(x), crs, ...): pipeline not found in #> PROJ-suggested candidate transformations st_bbox(meuse2) #>xminyminxmaxymax #> 178.605 329.714 181.390 333.611 I think this might be easier than manually adjusting the WKT representation of the CRS, but I'm not sure that this is really a recommended way to transform units since, for some reason, it does not modify the CRS of the objects which makes any analysis very difficult: all.equal(st_crs(meuse), st_crs(meuse2)) #> [1] TRUE library(mapview) mapview(meuse) # seems right mapview(meuse2) # completely off Andrea Il giorno mar 31 gen 2023 alle ore 12:33 Edzer Pebesma ha scritto: On 31/01/2023 12:12, Andrea Gilardi wrote: st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that. Dear all, I'm sorry if this is a trivial or stupid question but I was wondering if you could provide an example of using st_transform to apply a unit conversion transformation. For example, if I define something like library(sf) pt <- st_sfc(st_point(c(150, 500)), crs = 3003) # some place in North Italy st_crs(3003) # units are expressed in metres can I use st_transform to convert the unit measurement of pt from metres to km and preserve that information in the CRS? It seems so, using proj4strings (not recommended, but simple): library(sp) demo(meuse, echo = FALSE, ask = FALSE) library(sf) # Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE st_as_sf(meuse) |> st_bbox() # xmin ymin xmax ymax # 178605 329714 181390 333611 st_crs(meuse)$proj4string # [1] "+proj=sterea +lat_0=52.156160556 +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs" st_as_sf(meuse) |> st_transform("+proj=sterea +lat_0=52.156160556 +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=km +no_defs") |> st_bbox() #xminyminxmaxymax # 178.605 329.714 181.390 333.611 The better way would be to modify WKT representations of CRSs. Thanks Andrea Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma ha scritto: On 25/01/2023 18:42, Josiah Parry wrote: I wonder if `units::set_units()` is a better fit for the job https://r-quantities.github.io/units/articles/measurement_units_in_R.html It could definitely tell you which number to choose when going from, say, US feet to km. Coordinates in sf geometries are not stored as units objects, unit info is encoded in the CRS. st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that. On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter wrote: Is it safe to re-scale sf geometry coordinates from meters to kilometers using, for example: sfobj$geometry <- sfobj$geometry / 1000 It seems to work, but I understand too little about spatial data to know whether that practice is actually safe. I am working with spatial data in a continental-scale equidistant conic projection, and the units are meters. It seems that kilometers would be better suited stochastic partial differential equation modeling on a finite-element mesh, but maybe that is a moot point. Any and all advice will be much appreciated. Thanks! -- Steve Gutreuter [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
On 31/01/2023 12:12, Andrea Gilardi wrote: st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that. Dear all, I'm sorry if this is a trivial or stupid question but I was wondering if you could provide an example of using st_transform to apply a unit conversion transformation. For example, if I define something like library(sf) pt <- st_sfc(st_point(c(150, 500)), crs = 3003) # some place in North Italy st_crs(3003) # units are expressed in metres can I use st_transform to convert the unit measurement of pt from metres to km and preserve that information in the CRS? It seems so, using proj4strings (not recommended, but simple): library(sp) demo(meuse, echo = FALSE, ask = FALSE) library(sf) # Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE st_as_sf(meuse) |> st_bbox() # xmin ymin xmax ymax # 178605 329714 181390 333611 st_crs(meuse)$proj4string # [1] "+proj=sterea +lat_0=52.156160556 +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs" st_as_sf(meuse) |> st_transform("+proj=sterea +lat_0=52.156160556 +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=km +no_defs") |> st_bbox() #xminyminxmaxymax # 178.605 329.714 181.390 333.611 The better way would be to modify WKT representations of CRSs. Thanks Andrea Il giorno mer 25 gen 2023 alle ore 18:58 Edzer Pebesma ha scritto: On 25/01/2023 18:42, Josiah Parry wrote: I wonder if `units::set_units()` is a better fit for the job https://r-quantities.github.io/units/articles/measurement_units_in_R.html It could definitely tell you which number to choose when going from, say, US feet to km. Coordinates in sf geometries are not stored as units objects, unit info is encoded in the CRS. st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that. On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter wrote: Is it safe to re-scale sf geometry coordinates from meters to kilometers using, for example: sfobj$geometry <- sfobj$geometry / 1000 It seems to work, but I understand too little about spatial data to know whether that practice is actually safe. I am working with spatial data in a continental-scale equidistant conic projection, and the units are meters. It seems that kilometers would be better suited stochastic partial differential equation modeling on a finite-element mesh, but maybe that is a moot point. Any and all advice will be much appreciated. Thanks! -- Steve Gutreuter [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
On 25/01/2023 18:42, Josiah Parry wrote: I wonder if `units::set_units()` is a better fit for the job https://r-quantities.github.io/units/articles/measurement_units_in_R.html It could definitely tell you which number to choose when going from, say, US feet to km. Coordinates in sf geometries are not stored as units objects, unit info is encoded in the CRS. st_transform can do transformations and conversions between different CRSs, unit conversions is a special case of that. On Wed, Jan 25, 2023 at 12:37 PM Steve Gutreuter wrote: Is it safe to re-scale sf geometry coordinates from meters to kilometers using, for example: sfobj$geometry <- sfobj$geometry / 1000 It seems to work, but I understand too little about spatial data to know whether that practice is actually safe. I am working with spatial data in a continental-scale equidistant conic projection, and the units are meters. It seems that kilometers would be better suited stochastic partial differential equation modeling on a finite-element mesh, but maybe that is a moot point. Any and all advice will be much appreciated. Thanks! -- Steve Gutreuter [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Re-scale units of coordinates in sf geometry?
I guess the question is: "safe for what?" - if numerical errors are a problem without rescaling, then you should do it. You might be able to choose the rescaling factor such that these errors are minimized. On 25/01/2023 18:37, Steve Gutreuter wrote: Is it safe to re-scale sf geometry coordinates from meters to kilometers using, for example: sfobj$geometry <- sfobj$geometry / 1000 It seems to work, but I understand too little about spatial data to know whether that practice is actually safe. I am working with spatial data in a continental-scale equidistant conic projection, and the units are meters. It seems that kilometers would be better suited stochastic partial differential equation modeling on a finite-element mesh, but maybe that is a moot point. Any and all advice will be much appreciated. Thanks! -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] geodata package: CMIP6 climate model data downloading error
On 25/01/2023 15:48, Manuel Spínola wrote: Thank you very much Ákos. It didn't work either. Is there any other package that allows to download future climate change scenarios? CMIP6 data is being uploaded as Zarr files on the google cloud; a blog post that gives examples of downloading some parts of it, using packages stars and sf, is here: https://r-spatial.org/r/2022/09/13/zarr.html Manuel El mié, 25 ene 2023 a las 1:27, Bede-Fazekas Ákos () escribió: Dear Manuel, There are some broken links in the WorldClim website. If you find one, you can write to i...@worldclim.org. Although I have negative experiences, you should try this way. Anyway, the global, non-tiled raster can be downloaded from this link: https://geodata.ucdavis.edu/cmip6/30s/ACCESS-CM2/ssp126/wc2.1_30s_bioc_ACCESS-CM2_ssp126_2021-2040.tif Have a nice week, Ákos _ Ákos Bede-Fazekas Centre for Ecological Research, Hungary 2023.01.24. 23:29 keltezéssel, Manuel Spínola írta: Dear list members, I am trying to download CMIP6 climate model data using the R package geodata but I got this error message: bio <- cmip6_tile(lon = -84, lat = 10,"ACCESS-CM2", "126", "2021-2040", var="bioc", path = "climate_ckange") trying URL ' https://geodata.ucdavis.edu/cmip6/tiles/ACCESS-CM2/ssp126/wc2.1_30s_bioc_ACCESS-CM2_ssp126_2021-2040_tile-28.tif ' Error in utils::download.file(url = url, destfile = filename, quiet = quiet, : cannot open URL ' https://geodata.ucdavis.edu/cmip6/tiles/ACCESS-CM2/ssp126/wc2.1_30s_bioc_ACCESS-CM2_ssp126_2021-2040_tile-28.tif ' download failed I cannot even download the data using the worldclim website. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] retirement of rgdal, rgeos and maptools: second blog post
News and developments with regard to the retirement of R packages rgdal, rgeos and maptools (scheduled 2023) is found in the following blog post: https://r-spatial.org/r/2022/12/14/evolution2.html Follow-up questions can be directed to Roger or me, or to this list, or raised as issues on https://github.com/r-spatial/evolution Many regards, -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Convert geojson file to R
Interestingly, what seems to works is readLines('countrymasks.geojson') |> st_read() -> r with a warning: Warning message: In readLines("countrymasks.geojson") : incomplete final line found on 'countrymasks.geojson' On 29/11/2022 00:58, Miluji Sb wrote: Thank you. I will report this bug (I did not have the confidence to call this a bug before). Even using your code, I get the same output. structure(list(X = c(-67.3804401, -67.36091, -67.38058, -67.339709998, -67.3780199, -67.32211), Y = c(-55.565569996, -55.584009899, -55.600410001, -55.614969997, -55.63521, -55.640089997), L1 = c(1, 1, 1, 1, 1, 1), L2 = c(1, 1, 1, 1, 1, 1), L3 = c(1, 1, 1, 1, 1, 1)), row.names = c(NA, 6L), class = "data.frame") Thank you again. On Mon, Nov 28, 2022 at 11:13 PM Barry Rowlingson wrote: This seems to be a weird bug in `st_read`. If you read it with an SQL query that matches every row it works: js = st_read("./countrymasks.geojson", query="select * from countrymasks where 1 = 1") Reading query `select * from countrymasks where 1 = 1' from data source `/home/rowlings/Downloads/countrymasks.geojson' using driver `GeoJSON' Simple feature collection with 214 features and 15 fields Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: -180 ymin: -55.79439 xmax: 180 ymax: 83.62742 Geodetic CRS: WGS 84 But leave out the query and you get that C code level error. Another equivalent query would be "select * from countrymasks" (without the "where" clause) but this triggers the error too. Very odd. Worth reporting as a bug? Barry On Mon, Nov 28, 2022 at 8:08 PM Miluji Sb wrote: Thank you for reply. When I try using sf, I get the following error; Error in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : attempt to set index 210/210 in SET_STRING_ELT. Thanks again! On Mon, Nov 28, 2022 at 1:50 PM Josiah Parry wrote: You're going to want to read the file with sf. Try object <- sf::st_read("~countrymasks.geojson") On Mon, Nov 28, 2022 at 7:09 AM Miluji Sb wrote: Greetings everyone, I would like to convert the geojson file ( https://drive.google.com/file/d/18h3sOjZg5jp5euLTWRi5mC40Sja8TZDN/view?usp=sharing ) to a dataframe - essentially obtain which has coordinates matched to a country. I have tried the following; ### states <- geojsonsf::geojson_sf("~/countrymasks.geojson") geo <- geojsonsf::sf_geojson(states) sf <- sf::st_read(geo, quiet = T ) df <- as.data.frame(sf::st_coordinates(sf) ) ## But I get the following output, I am a bit lost. Any help will be highly appreciated. Best, structure(list(X = c(-67.3804401, -67.36091, -67.38058, -67.339709998, -67.3780199, -67.32211), Y = c(-55.565569996, -55.584009899, -55.600410001, -55.614969997, -55.63521, -55.640089997), L1 = c(1, 1, 1, 1, 1, 1), L2 = c(1, 1, 1, 1, 1, 1), L3 = c(1, 1, 1, 1, 1, 1)), row.names = c(NA, 6L), class = "data.frame") [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] LINESTRING to vector in R
st_as_text(st_geometry(rivers[1,])) On 22/10/2022 13:06, Nick Wray wrote: Hello I have downloaded a large shapefile dataset of UK rivers and I want to isolate (as an ordinary R string) the LINESTRING values for particular lines, corresponding to rivers Looking at the first line I can isolate the geometry by Hello I have downloaded a large shapefile dataset of UK rivers and I want to isolate (as an ordinary R string) the LINESTRING values for particular lines, corresponding to rivers Looking at the first line I can isolate the geometry by st_geometry(rivers[1,8]) Geometry set for 1 feature Geometry type: LINESTRING Dimension: XYZ Bounding box: xmin: 462010.6 ymin: 1213039 xmax: 462306.5 ymax: 1213199 z_range: zmin: 0 zmax: 0 Projected CRS: OSGB 1936 / British National Grid LINESTRING Z (462306.5 1213048 0, 462275.4 1213... What I need is all the values in the LINESTRING as a common or garden R vector, but I cannot find a way to do this. Does anyone know how? Thanks, Nick Wray [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Customizing levelplot coloring scheme in r
Look up the col.regions argument. On 12/10/2022 16:44, rain1290--- via R-sig-Geo wrote: I have climate data that I have plotted on a global map using a levelplot map. However, I am trying to adjust the default coloring scheme to allow for only blue and red shades to appear. Ideally, I am trying to showcase blue for the positive values, and red for the negative values, along with a whitish coloring near and at 0. The default coloring scheme isn't bad, but I think a blue-red-white coloring would allow the plot to appear more visually pleasing. Here is the code that I have now to create my current levelplot: #packages installed library(raster) library(ncdf4) library(maps) library(maptools) library(rasterVis) library(ggplot2) library(rgdal) library(sp) library(gridExtra) library(grid) library(RColorBrewer) #Using levelplot FPlot10 <- levelplot(FDifference5,margin=F,at=c(seq(-50,150,10)),pretty=TRUE, par.settings=mapTheme, main="Higher end") FM10 <- FPlot10 + latticeExtra::layer(sp.lines(world.outlines.sp)) This yields a range of nice default colors, but I'd like to adopt something with only shades of red and blue (with white at and around 0). I can do it easily in ggplot, but levelplot appears completely different in using color commands (if they exist). Is that even possible in levelplot? Thank you! [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Raster Data Management Advice
STAC is clearly the future of catalogues for spatial data, but not everyone has gotten there yet. Searching or browsing available STACs is helped by stac index, https://stacindex.org/ On 07/10/2022 18:43, Zivan Karaman wrote: Hi, Perhaps STAC <https://stacspec.org/en/> could help you? Best, Zivan On Fri, Oct 7, 2022 at 6:35 PM Alexander Ilich wrote: Hi, I was wondering if anyone has some advice on how to organize raster data so that it is easily queryable by various attributes (e.g. find me all the rasters of data type bathymetry, collected by this organization with 10m resolution or finer ). Currently we have data on a server organized often by when/where it was collected but that can make it difficult to find specific rasters that meet a certain criteria. I've created a table as a csv file on github <https://github.com/ailich/WFS_Multibeam_Metadata> where each row is a raster and it has various column attributes describing it (e.g. who collected it, what sonar was used, resolution, coordinate system, etc) and a path to the filename as a temporary solution, but I think some type of spatial database that would allow for querying and then reading into R as terra objects, as well as into QGIS and ArcGIS as layers for visualization would be optimal as multiple project members use these data. Tools I've come across that seem potentially useful include PostGIS and Geopackage, but I'm not entirely sure how to properly set them up or if they'd suit my needs. Any advice would be greatly appreciated. Thanks, Alex [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] cokriging simulation
Thanks Gerard, that has now been fixed in https://github.com/r-spatial/gstat/issues/106 The problem was probably introduced in 2016 when I swapped out the shipped matrix library with R's native one. The two differed in how a Choleski decomposition was returned . On 27/04/2022 21:29, Heuvelink, Gerard wrote: Dear List, We wish to create simulations from a cokriging model. I used to be able to do that in the past, but now I am struggling to get the correlations preserved. I dimplified the problem to a simple case with cross-correlation 0.8 (see script below), but simulations show no correlation at all. What do I do wrong? Thank you for your help, Gerard library(gstat) set.seed(12345) d <- data.frame(0,0,0,0) names(d) <- c("x","y","lmCres","lmNres") coordinates(d) = ~x+y nugC <- 1; nugN <- 1; nugCN <- 0.8 psillC <- 5; psillN <- 5; psillCN <- 4 range <- 10 g_l <- gstat(id=c("lmCres"),formula=lmCres~1,data=d,beta=0,nmax=100, model=vgm(psillC,"Sph",range=10,nugC)) g_l <- gstat(g_l,id="lmNres",formula=lmNres~1,data=d,beta=0,nmax=100, model=vgm(psillN,"Sph",range=10,nugN)) g_l <- gstat(g_l,id=c("lmCres","lmNres"),model=vgm(psillCN,"Sph",range=10,nugCN)) df <- data.frame(x=1,y=1) coordinates(df) <- ~x+y nsim <- 250 sim <- predict(g_l,df,nsim=nsim) simdf <- as.data.frame(sim) plot(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)])) cor(as.numeric(simdf[1,3:(nsim+2)]),as.numeric(simdf[1,(nsim+3):((2*nsim)+2)])) [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
On 14/11/2021 13:59, Edzer Pebesma wrote: Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5 datasets, one containing the raster of the longitude and latitude, the other rasters with the attributres (LST) but no coordinates. I could get a (rough: factor 50 downsampling) plot with stars using: library(stars) lat = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude") lon = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude") var = "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST" lst = read_stars(var, curvilinear = list(x = lon, y = lat)) plot(lst, downsample = 50, axes = TRUE, reset = FALSE) library(rnaturalearth) ne = ne_countries(returnclass = "sf") plot(st_geometry(ne), add = TRUE, border = 'yellow') # california border seems to make sense but I'm not sure how to make a regular raster (in some CRS) out of this for exporting to GeoTIFF. This might be useful: https://github.com/r-spatial/stars/issues/386 but converting to polygons will be incredibly slow and memory hungry. On 14/11/2021 11:00, Edzer Pebesma wrote: This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs : NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. hth On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) r and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it as a raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer and if possible save it as GeoTIFF file format to be further processed with the raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier -- Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-
Re: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
Gabriel shared the data with me (off-list, 1Gb) and it contains two hdf5 datasets, one containing the raster of the longitude and latitude, the other rasters with the attributres (LST) but no coordinates. I could get a (rough: factor 50 downsampling) plot with stars using: library(stars) lat = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/latitude") lon = read_stars("HDF5:ECOSTRESS_L1B_GEO_11712_012_20200730T061029_0601_01.h5://Geolocation/longitude") var = "HDF5:ECOSTRESS_L2_LSTE_11712_012_20200730T061029_0601_01.h5://SDS/LST" lst = read_stars(var, curvilinear = list(x = lon, y = lat)) plot(lst, downsample = 50, axes = TRUE, reset = FALSE) library(rnaturalearth) ne = ne_countries(returnclass = "sf") plot(st_geometry(ne), add = TRUE, border = 'yellow') # california border seems to make sense but I'm not sure how to make a regular raster (in some CRS) out of this for exporting to GeoTIFF. On 14/11/2021 11:00, Edzer Pebesma wrote: This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs : NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. hth On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) r and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it as a raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer and if possible save it as GeoTIFF file format to be further processed with the raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier -- Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Question on hierarchical data format 5 .h5 file format
This might be because the grid is not regular but possibly curvilinear, having two raster layers with the lon & lat values for each pixel. Can you share a sample data set? On 14/11/2021 10:47, Gabriel Cotlier wrote: Dear Michael, Following your advice I have tired : library(terra) library(raster) library(dplyr) r = terra::rast(file.choose()) plot(r$LST) r_lst = r$LST r_lst %>% raster() plot(r_lst) Then I effectively get a raster layer format from the raster package as I wanted, but looking at the plot and the raster layer information--see bellow-- the scene appears to need both latitude and longitude for correct geolocation, how can I do this procedure? r_lst %>% raster() class : RasterLayer dimensions : 5632, 5400, 30412800 (nrow, ncol, ncell) resolution : 1, 1 (x, y) extent : 0, 5400, 0, 5632 (xmin, xmax, ymin, ymax) crs: NA source : LST names : LST values : 0, 65535 (min, max) The extent and the resolution are not the correct ones. Maybe in the same .h5 is all the data for georeferencing it although I could not see how to get to it and apply it to the layer. There is also another file with complementary data related to geolocation but how can I query it, and get the data for geolocation and apply it to the raster layer?. The objective is to save it as a GeoTIFF file to be opened in QGIS with correct geolocation. What do you suggest for converting the .h5 file to a georeferenced Raster Layer in GeoTIFF format ? Thanks a lot. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:22 PM Gabriel Cotlier wrote: Dear Michael Summer, Thanks a lot for your swift reply. I will give a try to the procedure you mentioned. Kind regards, Gabriel On Sat, Nov 13, 2021 at 12:14 PM Michael Sumner wrote: try terra::rast() or stars::read_stars() on the file, both use gdal to interrogate and read. What you get depends on the structure of the files, variables interpreted as subdatasets. stars is more general but sees variables as 2D arrays with bands, terra calls these bands layers. Alternatively try raster::raster() or RNetCDF or tidync which might be more straightforward than rhdf5 but all are quite different to each other and to the gdal view. hth On Sat, 13 Nov 2021, 18:56 Gabriel Cotlier, wrote: Hello, I would like to be able to visualize / plot a specific data layer within the structure a hierarchical data file format .h5 in R. I could observe the structure as follows : library(rhdf5) r= h5ls(file.choose()) r and also in this way identify the Group and Name of the data to be plotted, however I've tried I could not access the data itself, nor plot it as a raster layer with its corresponding geographic coordinates. I would like to be able to select such data, plot it as a raster layer and if possible save it as GeoTIFF file format to be further processed with the raster package as a common raster layer. If any small guidance example, reprex, of this procedure is possible it would be appreciated. Thanks a lot in advance Kind regards, Gabriel Cotlier -- Gabriel Cotlier, PhD Haifa Research Center for Theoretical Physics and Astrophysics (HCTPA) University of Haifa Israel [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Deprecated packages
Hi Erin, I am also thinking about it. With Roger's retirement, rgdal, rgeos and maptools will effectively go away, by the end of 2023. We are working on contingency plans and instructions for package maintainers, and will share on this list as soon as there is something to share. For automap, I have no clue; it hasn't been updated for a long time but doesn't depend on one of these three. On 02/11/2021 20:38, Erin Hodgess wrote: Hello there! Are both maptools and automap going away, please? Also, is anyone putting together a sort of substitution list for old packages going to new packages, please? Just wondering. I was thinking about it if no one else is. Thanks, Erin -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Aggregation and disaggregation in sf
I'm not so familiar with the terminology you use, but the function you'd most likely want to look into is st_interpolate_aw(). On 14/10/2021 12:11, Sean Trende wrote: Thank you all for developing this package and taking a huge amount of time to answer questions. I had a question that I found related answers to, but nothing directly on point. Perhaps I'm just not being sufficiently creative with my R coding, and I apologize if I missed the "explainer" elsewhere. A common task is aggregating and disaggregating within a larger polygon set. I believe that I know how to do the following in the geomander package, but am wondering if there is a "pure" sf solution. There is a simple version of the task, and a more complex one. The simple version of task: You have two shapefiles, one of census blocks and one of precinct results. You have a polygon for a precinct that split 100 votes to 100 votes for candidates A and B in the preceding election. The precinct wholly contains three census blocks with respective populations of 100, 200 and 300. The task is to identify those blocks that exist within the precinct from the broader block shapefile and then to apportion the votes among these census blocks proportional to population, such that you would have 16.6, 33.3, and 50 votes recorded at the block level for each candidate, respectively. The more difficult version: Same facts, but now there is a 4th block that is split between precincts. So there has to be some sort of way to acknowledge these blocks and split them, likely in an arbitrary way between the portion within the precinct and without (e.g., 50-50 or by land area within each precinct) I wouldn't impose upon you to actually write code, but if I could be directed toward the most relevant functions (or told to go play around in ArcGIS) it would be incredibly helpful. Best regards, Sean ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Problem to retrieve/plot selected bands in a stars proxy object
Thanks for the clear report, Loïc; this should now be fixed in the version on github. On 14/09/2021 14:55, Loïc Valéry wrote: Dear list members, I hope this message finds you well. I submit to you a problem with which I have been struggling for several days without finding the beginning of a solution. I would like to select some bands of a 'stars proxy' object that I have previously created by merging two images with the function read_stars(c(image1, image2), proxy = true, along = "band") in order to plot them. My problem is that R returns the following error message when trying to plot the selected bands: Error in x[[i]][, , 10:12, , drop = FALSE] : subscript out of bounds Here is a short REPREX: library(stars) tif_1 <- system.file("tif/L7_ETMs.tif", package = "stars") tif_2 <- system.file("tif/L7_ETMs.tif", package = "stars") tif_merge <- read_stars(c(tif_1,tif_2), proxy = TRUE, along = "band") plot(tif_merge[,,,10:12], rgb = 1:3) #> Error in x[[i]][, , 10:12, , drop = FALSE]: indice hors limites I guess the problem is that the "tif_merge" object has different dimensions before and after the promise evaluation: dim(tif_merge) #>xy band #> 349 352 12 dim(st_as_stars(tif_merge)) #> x yband new_dim #> 349 352 6 2 However, this did not help me to find any solution! Thank you in advance for your help. Cheers, Loïc ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Fwd: GDAL Error 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326
I would say this indicates a broken installation of PROJ on your system, or you have more than one instance (and different versions) of PROJ installed. On 10/09/2021 12:20, Anna Mujal Colilles wrote: Hi, I am a new R user trying to transform lat-lon data to utm. I've been doing so in two different ways and both give me the same error. My machine is a Linux Ubuntu 18.04 system. The first way I wanted to transform them was: sp <- SpatialPoints(df[,c("longitude","latitude")],proj4string=CRS("+init=epsg:4326 +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 ")) sp2<-spTransform(sp, CRS("+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")) the second option was sf_sample <- sf::st_as_sf(coord_sample, coords = c("longitude", "latitude"), crs = 4326) However, I keep having the same problem over and over again related to the crs init 4326 Warning message: In CPL_crs_from_input(x) : GDAL Error 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326: SQLite error on SELECT name, ellipsoid_auth_name, ellipsoid_code, prime_meridian_auth_name, prime_meridian_code, area_of_use_auth_name, area_of_use_code, publication_date, deprecated FROM geodetic_datum WHERE auth_name = ? AND code = ?: no such column: publication_date can anyone help me please? Thank you, Anna [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Issue with st_warp, gdal and NAs
Dear Julian, you can do that by passing a non-NA value as no_data_value, e.g. by ag3 <- st_warp(src = s1, dest = st_as_stars(st_bbox(s1), dx = 1, dy = 1), use_gdal = TRUE, method = "average", no_data_value = -) On 16/07/2021 12:17, Julian M. Burgos wrote: Dear list, My apologies for posting this again. I am not sure if it went through the first time. I am having a bit of trouble using st_warp() to aggregate star objects with the option "use_gdal = TRUE". I am using the "average" method, and I want cells with NA values to be ignored when computing the averages. To demonstrate, I will make a small raster for testing, and will include cells with NA values: #- library(raster) # Create raster for testing myfile <- "~/myraster.tif" myraster <- raster(xmn = 0, ymn = 0, xmx = 10, ymx = 10, resolution = 0.1, crs = 4326, vals = sample(1:20, size = 1, replace = TRUE)) myraster[sample(1:1, 200)] <- NA writeRaster(myraster, filename = myfile, overwrite = TRUE) #- With the raster package, I can aggregate the cells ignoring the NA values like this: r1 <- raster(myfile) ag1 <- raster::aggregate(x = raster(myfile), fact = 10, fun = mean, na.rm = TRUE) ...which is what I want. I can also aggregate the cells including the NAs, like this: ag2 <- raster::aggregate(x = raster(myfile), fact = 10, fun = mean, na.rm = FALSE) which produces a raster with a lots of gaps. Now, to do a similar aggregation with the stars package I can do this: s1 <- read_stars(myfile) ag3 <- st_warp(src = s1, dest = st_as_stars(st_bbox(s1), dx = 1, dy = 1), use_gdal = TRUE, method = "average") The resulting raster (ag3) has gaps like ag2, which means that the NAs are not being ignored when computing the averages. With use_gdal_TRUE and method="average" we are using gdalwarp with the average resampling method, which according to the documentation "computes the weighted average of all non-NODATA contributing pixels.". So I guess the problem is that the NAs are not being recognized as NODATA pixels. Does anybody know how to solve this? Many thanks, Julian -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] problems installing stars and tmap in fedora 34
Hi Vero, I would start with reinstalling package lwgeom. On 21/06/2021 21:30, Veronica Andreo wrote: Dear all, I have just updated my laptop to fedora 34 and the only problems I've encountered so far are related to stars and tmap updates. This is the error I get in both cases: Error in dyn.load(file, DLLpath = DLLpath, ...) : unable to load shared object '/home/veroandreo/R/x86_64-redhat-linux-gnu-library/4.0/lwgeom/libs/lwgeom.so': libproj.so.15: cannot open shared object file: No such file or directory Calls: ... asNamespace -> loadNamespace -> library.dynam -> dyn.load Execution halted ERROR: lazy loading failed for package ‘stars’ * removing ‘/home/veroandreo/R/x86_64-redhat-linux-gnu-library/4.0/stars’ Warning in install.packages : installation of package ‘stars’ had non-zero exit status Searching around I found a (maybe related) issue reported for sf ( https://github.com/r-spatial/sf/issues/1670), but here sf installed just fine and also rgdal. GRASS also compiles well, btw, and gdal-devel and proj-devel are in place. I tried the solution proposed in the issue (`remotes::install_github("r-spatial/stars", configure.args= '--with-proj-lib=/lib64')`), but it does not work either. I guess it's because of this one: https://bugzilla.redhat.com/show_bug.cgi?id=1958191. Right? Is there anything else I could try/test/report? I really appreciate your help :) Best, Vero [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] stars st_transform gives NA to the offset and delta parameters
On 07/06/2021 01:16, Manuel Spínola wrote: Dear list members, When using st_transform in a stars object the offset and delta parameters become NA, is this an expected behavior? Yes, it also says that you got a curvilinear grid as a result, which is a grid with non-constant cell size (delta) and non-constant coordinate of the side (offset); see https://keen-swartz-3146c4.netlify.app/intro.html#raster-types If you want a regular grid in the new CRS, use st_warp(), preferably with the target grid template. geomatrix = system.file("tif/geomatrix.tif", package = "stars") x = read_stars(geomatrix) new = st_crs(4326) y = st_transform(x, new) y *Manuel Spínola, Ph.D.* Instituto Internacional en Conservación y Manejo de Vida Silvestre Universidad Nacional Apartado 1350-3000 Heredia COSTA RICA mspin...@una.cr mspinol...@gmail.com Teléfono: (506) 8706 - 4662 Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/> Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/> [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Using facet_wrap with geom_stars for a stacked raster
On 01/06/2021 19:15, Manuel Spínola wrote: Thank you very much Edzer. What about the labels of the stripes? If I try to change the labels, they come out as NA r <- read_stars(system.file("external/test.grd", package="raster")) s1 <- c(r, r*2) s2 <- merge(s1) n <- c("raster 01", "raster 02") ggplot() + geom_stars(data = s2) + facet_wrap(~attributes, labeller = labeller(attributes = n)) + coord_equal() s3 = st_set_dimensions(s2, "attributes", values = c("label_A", "label_B")) ggplot() + geom_stars(data = s3) + facet_wrap(~attributes) + coord_equal() El mar, 1 jun 2021 a las 10:54, Edzer Pebesma (mailto:edzer.pebe...@uni-muenster.de>>) escribió: On 01/06/2021 18:43, Manuel Spínola wrote: > Dear list members, > > I am trying to use facet_wrap with geom_stars but I don't know how to > specify the variable argument in the facet_wrap function (I know is not > "band"), and also, how to change the labels of the strip names? > > r <- read_stars(system.file("external/test.grd", package="raster")) > > s1 <- c(r, r*2) > > ggplot() + > geom_stars(data = s1) + > coord_equal() + > facet_wrap(~ band) > Try: ggplot() + geom_stars(data = merge(s1)) + f facet_wrap(~attributes) + coord_equal() c() binds attributes, merge() merges them into a dimension. -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ 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 <https://stat.ethz.ch/mailman/listinfo/r-sig-geo> -- *Manuel Spínola, Ph.D.* Instituto Internacional en Conservación y Manejo de Vida Silvestre Universidad Nacional Apartado 1350-3000 Heredia COSTA RICA mspin...@una.cr <mailto:mspin...@una.ac.cr> mspinol...@gmail.com <mailto:mspinol...@gmail.com> Teléfono: (506) 8706 - 4662 Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/> Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/> -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Using facet_wrap with geom_stars for a stacked raster
On 01/06/2021 18:43, Manuel Spínola wrote: Dear list members, I am trying to use facet_wrap with geom_stars but I don't know how to specify the variable argument in the facet_wrap function (I know is not "band"), and also, how to change the labels of the strip names? r <- read_stars(system.file("external/test.grd", package="raster")) s1 <- c(r, r*2) ggplot() + geom_stars(data = s1) + coord_equal() + facet_wrap(~ band) Try: ggplot() + geom_stars(data = merge(s1)) + f facet_wrap(~attributes) + coord_equal() c() binds attributes, merge() merges them into a dimension. -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Identifying adjacent polygons within larger polygon groupings
On 01/04/2021 17:01, Sean Trende wrote: Hello, I am using sf for the project below, but am open to non-sf solutions. Given an sf data.frame of polygons sorted into n contiguous groups (e.g., a column in the df identifies each polygon as a member of group (1, . . ., n), where within each group all polygons are neighbors), is there a simple way to identify all polygons in, say, Group 5 that are neighbors to polygons in, say, Group 8 (assuming groups 5 and 8 are neighbors)? E.g., if I have an sf data.frame of U.S. counties with a column that identifies the state in which each county is located, is there a way to identify all counties on the border between, say, Kansas and Nebraska? Yes, look into using st_intersects (or st_touches) with x the counties of Kansas and y the counties of Nebraska. What you get back is a list of length length(x), with integer vectors, the i-th vector containing the indexes of Nebraska counties touching/intersecting county i of Kansas. Thank you so much, Sean ___ R-sig-Geo mailing list R-sig-Geo@r-project.org 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial layers for Europe at 30-m available as Cloud Optimized GeoTiffs (corrected)
Thanks, Tom. How can one use the .qml file you point to in an R workflow to get category labels and colors for the COG data you shared originally? On 16/03/2021 12:24, Tomislav Hengl wrote: Hi Edzer thanks for spotting this. Cloud Optimized GeoTiffs seems to be somewhat complicated when it comes to adding metadata etc in the tif file directly. For example, the COG validator reports: # rio cogeo validate /data/raster/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif he following errors were found: - The offset of the main IFD should be < 300. It is 1681443308 instead - The offset of the IFD for overview of index 0 is 750, whereas it should be greater than the one of the main image, which is at byte 1681443308 /data/raster/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif is NOT a valid cloud optimized GeoTIFF But it seems that the tif fails the check, but works fine on the end. We are still only testing if this has any effect on speed of access etc. The worst case scenario is that users can not access the data because of some small glitch in the tif header. To produce a legend for land cover maps, we recommend using: https://gitlab.com/geoharmonizer_inea/spatial-layers/-/blob/master/map-style/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2000_eumap_epsg3035_v0.1.qml You can of course always visually check what you get via: https://maps.opendatascience.eu HTH, T. Hengl https://opengeohub.org/about On 3/15/21 4:40 PM, Edzer Pebesma wrote: Hi Tom, great work! When reading this either with terra or stars, as in in.tif = "/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif; library(terra) # terra version 1.1.5 plot(rast(in.tif)) #library(raster) #plot(raster(in.tif)) library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1 plot(read_stars(in.tif)) I see that both plots show a continuous raster, not a categorical. We spent quite a bit of time recently in stars & tmap development to get handling and plotting of categorical rasters right (including those having a color table), but this file doesn't give proper access to the categories. The metadata tags do have them: stars::gdal_metadata(in.tif) # [1] "1=111 - Urban fabric" # [2] "10=212 - Permanently irrigated arable land" # [3] "11=213 - Rice fields" # [4] "12=221 - Vineyards" # [5] "13=222 - Fruit trees and berry plantations" # [6] "14=223 - Olive groves" # [7] "15=231 - Pastures" # [8] "16=311 - Broad-leaved forest" # [9] "17=312 - Coniferous forest" # [10] "18=321 - Natural grasslands" # [11] "19=322 - Moors and heathland" # [12] "2=122 - Road and rail networks and associated land" # [13] "20=323 - Sclerophyllous vegetation" # [14] "21=324 - Transitional woodland-shrub" # [15] "22=331 - Beaches, dunes, sands" # [16] "23=332 - Bare rocks" # [17] "24=333 - Sparsely vegetated areas" # [18] "25=334 - Burnt areas" # [19] "26=335 - Glaciers and perpetual snow" # [20] "27=411 - Inland wetlands" # [21] "28=421 - Maritime wetlands" # [22] "29=511 - Water courses" # [23] "3=123 - Port areas" # [24] "30=512 - Water bodies" # [25] "31=521 - Coastal lagoons" # [26] "32=522 - Estuaries" # [27] "33=523 - Sea and ocean" # [28] "4=124 - Airports" # [29] "5=131 - Mineral extraction sites" # [30] "6=132 - Dump sites" # [31] "7=133 - Construction sites" # [32] "8=141 - Green urban areas" # [33] "9=211 - Non-irrigated arable land" # [34] "AREA_OR_POINT=Area" but GDAL doesn't give access to those in a programmatic way. I've tried to add a .aux.xml file with the table, this worked locally (for both stars - after using droplevels() - and terra) and might as well work over the /vsicurl connection. File attached. Many regards, On 02/03/2021 18:30, Tomislav Hengl wrote: We have mapped land cover classes for the 2000-2019 period for continental Europe at 30-m resolution using spatiotemporal Machine Learning (we used R and python for modeling). Explore the dynamic EU landscapes on your palm using the ODS-Europe viewer: https://maps.opendatascience.eu To access almost 10TB of data using R you use the terra or similar packages e.g.: R> library(terra) R> in.tif = "/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif; R> tif = rast(in.tif) From here you can use any native operation e.g. to crop some polygon or resample / aggregate values (there is no need
Re: [R-sig-Geo] Spatial layers for Europe at 30-m available as Cloud Optimized GeoTiffs (corrected)
Hi Tom, great work! When reading this either with terra or stars, as in in.tif = "/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif; library(terra) # terra version 1.1.5 plot(rast(in.tif)) #library(raster) #plot(raster(in.tif)) library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1 plot(read_stars(in.tif)) I see that both plots show a continuous raster, not a categorical. We spent quite a bit of time recently in stars & tmap development to get handling and plotting of categorical rasters right (including those having a color table), but this file doesn't give proper access to the categories. The metadata tags do have them: stars::gdal_metadata(in.tif) # [1] "1=111 - Urban fabric" # [2] "10=212 - Permanently irrigated arable land" # [3] "11=213 - Rice fields" # [4] "12=221 - Vineyards" # [5] "13=222 - Fruit trees and berry plantations" # [6] "14=223 - Olive groves" # [7] "15=231 - Pastures" # [8] "16=311 - Broad-leaved forest" # [9] "17=312 - Coniferous forest" # [10] "18=321 - Natural grasslands" # [11] "19=322 - Moors and heathland" # [12] "2=122 - Road and rail networks and associated land" # [13] "20=323 - Sclerophyllous vegetation" # [14] "21=324 - Transitional woodland-shrub" # [15] "22=331 - Beaches, dunes, sands" # [16] "23=332 - Bare rocks" # [17] "24=333 - Sparsely vegetated areas" # [18] "25=334 - Burnt areas" # [19] "26=335 - Glaciers and perpetual snow" # [20] "27=411 - Inland wetlands" # [21] "28=421 - Maritime wetlands" # [22] "29=511 - Water courses" # [23] "3=123 - Port areas" # [24] "30=512 - Water bodies" # [25] "31=521 - Coastal lagoons" # [26] "32=522 - Estuaries" # [27] "33=523 - Sea and ocean" # [28] "4=124 - Airports" # [29] "5=131 - Mineral extraction sites" # [30] "6=132 - Dump sites" # [31] "7=133 - Construction sites" # [32] "8=141 - Green urban areas" # [33] "9=211 - Non-irrigated arable land" # [34] "AREA_OR_POINT=Area" but GDAL doesn't give access to those in a programmatic way. I've tried to add a .aux.xml file with the table, this worked locally (for both stars - after using droplevels() - and terra) and might as well work over the /vsicurl connection. File attached. Many regards, On 02/03/2021 18:30, Tomislav Hengl wrote: We have mapped land cover classes for the 2000-2019 period for continental Europe at 30-m resolution using spatiotemporal Machine Learning (we used R and python for modeling). Explore the dynamic EU landscapes on your palm using the ODS-Europe viewer: https://maps.opendatascience.eu To access almost 10TB of data using R you use the terra or similar packages e.g.: R> library(terra) R> in.tif = "/vsicurl/http://s3.eu-central-1.wasabisys.com/eumap/lcv/lcv_landcover.hcl_lucas.corine.rf_p_30m_0..0cm_2019_eumap_epsg3035_v0.1.tif; R> tif = rast(in.tif) From here you can use any native operation e.g. to crop some polygon or resample / aggregate values (there is no need to download whole data sets). A detailed tutorial on how to work with Cloud Optimized GeoTiffs is available here: https://gitlab.com/openlandmap/global-layers/-/blob/master/tutorial/OpenLandMap_COG_tutorial.md. Complete list of Cloud Optimized GeoTiffs we produced so far for Europe is available here: https://gitlab.com/geoharmonizer_inea/eumap/-/blob/master/gh_raster_layers.csv If not otherwise specified, the data available on this portal is licensed under the Open Data Commons Open Database License <https://opendatacommons.org/licenses/odbl/> (ODbL) and/or Creative Commons Attribution-ShareAlike 4.0 <https://creativecommons.org/licenses/by-sa/4.0/legalcode> and/or Creative Commons Attribution 4.0 <https://creativecommons.org/licenses/by/4.0/legalcode> International license (CC BY). Read more in: https://opengeohub.medium.com/europe-from-above-space-time-machine-learning-reveals-our-changing-environment-1b05cb7be520 If you experience any technical problems or if you discover a bug, please report via: https://gitlab.com/geoharmonizer_inea/spatial-layers/-/issues T. Hengl https://opengeohub.org/about ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 111 - Urban fabric 212 - Permanently irrigated arable land 213 - Rice fields 221 - Vineyards 222 - Fruit trees and berry plan
Re: [R-sig-Geo] Species distribution modeling with presence-only data using log-Gaussian Cox process
A way out of this may be to make a local copy of the lgcp sources, remove the rpanel dependency from DESCRIPTION, and install it from there. You won't get the functions that need rpanel, but those might not be essential. On 09/02/2021 20:09, Manuel Spínola wrote: Thank you very much Edzer. I tried to install lgcp, but I got an error message regarding BWidget. I think that Bwidget is an external software and I don't know how to install it in MacOS (Big Sur 11.2). Manuel This is the error that I get when trying to install lgcp in R 4.0.3: xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun Warning message: In system2("/usr/bin/otool", c("-L", shQuote(DSO)), stdout = TRUE) : running command ''/usr/bin/otool' -L '/Library/Frameworks/R.framework/Resources/library/tcltk/libs//tcltk.so'' had status 1 Error in structure(.External(.C_dotTcl, ...), class = "tclObj") : [tcl] can't find package BWidget. Error: unable to load R code in package ‘rpanel’ Execution halted El mar, 9 feb 2021 a las 12:51, Edzer Pebesma (mailto:edzer.pebe...@uni-muenster.de>>) escribió: Have you tried package lgcp? On 09/02/2021 17:26, Manuel Spínola wrote: > Dear list members, > > I would like to know of any R package that can be used for modeling species > distribution with presence-only data using log-Gaussian Cox process. > > Manuel > ___ 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 -- *Manuel Spínola, Ph.D.* Instituto Internacional en Conservación y Manejo de Vida Silvestre Universidad Nacional Apartado 1350-3000 Heredia COSTA RICA mspin...@una.cr <mailto:mspin...@una.ac.cr> mspinol...@gmail.com <mailto:mspinol...@gmail.com> Teléfono: (506) 8706 - 4662 Personal website: Lobito de río <https://sites.google.com/site/lobitoderio/> Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/> ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Species distribution modeling with presence-only data using log-Gaussian Cox process
Have you tried package lgcp? On 09/02/2021 17:26, Manuel Spínola wrote: Dear list members, I would like to know of any R package that can be used for modeling species distribution with presence-only data using log-Gaussian Cox process. Manuel ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial joins between sf and stars objects
On 03/02/2021 20:14, Julian M. Burgos wrote: Although if I use bind_cols the result is not an sf object right, this is because bind_cols is not a generic (its first argument is ... , which you can't dispatch on) so doesn't have an associated sf method. dplyr offers a mechanism to work around this by providing an sf method for dplyr_reconstruct - something we could try, but then the dplyr authors mention that this is experimental and considered a stop-gap measure, so that might work now but stop working any time. However, sf3 <- cbind(st_extract(st1, sf1), st_drop_geometry(sf1)) does return an sf object (note the changed order of arguments) Julian M. Burgos writes: Thanks Edzer, I should have checked if st_extract kept the order of the rows. In this case, a way to do it more generally is: sf1 <- bind_cols(st_drop_geometry(sf1), st_extract(st1, sf1)) Edzer Pebesma writes: st_join does a spatial join, which comes at a cost; since st_extract doesn't change the row order of the sf object, you could as well use st_extract(st1, sf1) %>% mutate(nums = sf1$nums) On 03/02/2021 17:10, Julian M. Burgos wrote: Dear list, What would be the best way to do a spatial join between a sf object (of POINT geometry) and a star object, so the sf object gets the values of the corresponding pixels in the star object? I have done using a combination of st_join and st_extract, but it feels a bit clunky. Is there a better way? Here is what I have done: ## -- ## Load a stars object st1 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>% slice(band, 1) ## Create an sf object set.seed(100) bb <- st_bbox(st1) sf1 <- tibble(lon = runif(n = 10, min = bb[1], max = bb[3]), lat = runif(n = 10, min = bb[2], max = bb[4]), nums = 1:10) %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(st1)) ## Do the spatial join sf1 <- st_extract(st1, sf1) %>% st_join(sf1) ## -- Thanks, Julian -- Julian Mariano Burgos, PhD Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/ Marine and Freshwater Research Institute Botnsjávarsviðs / Demersal Division Fornubúðir 5, IS-220 Hafnarfjörður, Iceland http://www.hafogvatn.is/ Sími/Telephone : +354-5752037 Netfang/Email: julian.bur...@hafogvatn.is ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geodata=04%7C01%7C%7C7b9230c283dd4883458a08d8c876e7d1%7C8e105b94435e4303a61063620dbe162b%7C0%7C0%7C637479760450605041%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000sdata=d%2FY%2FuSVIou5juUoKgSwCZc4EZ5OdEJRCg1YZHH%2FEYGc%3Dreserved=0 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geodata=04%7C01%7C%7C7b9230c283dd4883458a08d8c876e7d1%7C8e105b94435e4303a61063620dbe162b%7C0%7C0%7C637479760450605041%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000sdata=d%2FY%2FuSVIou5juUoKgSwCZc4EZ5OdEJRCg1YZHH%2FEYGc%3Dreserved=0 -- Julian Mariano Burgos, PhD Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/ Marine and Freshwater Research Institute Botnsjávarsviðs / Demersal Division Fornubúðir 5, IS-220 Hafnarfjörður, Iceland www.hafogvatn.is Sími/Telephone : +354-5752037 Netfang/Email: julian.bur...@hafogvatn.is ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial joins between sf and stars objects
st_join does a spatial join, which comes at a cost; since st_extract doesn't change the row order of the sf object, you could as well use st_extract(st1, sf1) %>% mutate(nums = sf1$nums) On 03/02/2021 17:10, Julian M. Burgos wrote: Dear list, What would be the best way to do a spatial join between a sf object (of POINT geometry) and a star object, so the sf object gets the values of the corresponding pixels in the star object? I have done using a combination of st_join and st_extract, but it feels a bit clunky. Is there a better way? Here is what I have done: ## -- ## Load a stars object st1 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>% slice(band, 1) ## Create an sf object set.seed(100) bb <- st_bbox(st1) sf1 <- tibble(lon = runif(n = 10, min = bb[1], max = bb[3]), lat = runif(n = 10, min = bb[2], max = bb[4]), nums = 1:10) %>% st_as_sf(coords = c("lon", "lat"), crs = st_crs(st1)) ## Do the spatial join sf1 <- st_extract(st1, sf1) %>% st_join(sf1) ## -- Thanks, Julian -- Julian Mariano Burgos, PhD Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/ Marine and Freshwater Research Institute Botnsjávarsviðs / Demersal Division Fornubúðir 5, IS-220 Hafnarfjörður, Iceland www.hafogvatn.is Sími/Telephone : +354-5752037 Netfang/Email: julian.bur...@hafogvatn.is ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] st_segmentize across east and west hemispheres
You may want to look into st_wrap_dateline(), which cuts LINESTRING and POLYGON geometries in multi-part equivalents where parts do not cross the antimeridian. As in: seg2 <- st_segmentize(st_wrap_dateline(sf), units::set_units(1000, km)) On 29/12/2020 00:24, Amanda Rehbein wrote: Dear r-sig-geo list, I have a package called raytracing for calculating atmospheric Rossby wave paths. I need to get segments of the great circle or routes from some geographical coordinates. st_segmentize is calculating them correctly. However, when I need to connect two points in different hemispheres, east and west, it creates an unwanted horizontal line, as shown in the following example. Is it possible (and correct) to avoid or remove this horizontal line? library(sf) m <- rbind(c(100,-50), c(-100,50)) sf <- st_sf(a=1, geom=st_sfc(st_linestring(m)), crs = 4326) seg <- st_segmentize(sf, units::set_units(1000, km)) plot(seg, axes = TRUE, reset = FALSE, type = "p", pch = 16) plot(seg$geom, add = TRUE, col = "red") text(x = m[, 1], y = m[, 2] - 7, label = 1:2, col = "blue") Many thanks. [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Accurate spatial database for country administrative boundaries
Have you looked at package rnaturalearth? There's a web site for the dataset it interfaces: https://www.naturalearthdata.com/ - it looks like there's more information about its origins and processing that went on than behind GADM. On 12/1/20 4:26 PM, Manuel Spínola wrote: Dear list members, I am looking for a global accurate spatial database for countries' administrative boundaries. I tried the raster package, using the getData function and the GADM database and rgeoboundaries for the geoboundaries database, but I get differences between them. Is there any other spatial database for administrative boundaries accessible from R, so I can compare and decide which one is the most accurate of the open access spatial database for the countries I am interested? Thank you very much in advance. Manuel -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] raster::levels() not working in packaged function.
On 10/31/20 1:24 PM, nevil amos wrote: Apologies, I cannot see how to make a rero for this issue. I have a function that uses levels(r) tor return the RAT of a raster "r" when the function is sourced from a script source(".\R\function.r") it works fine. when the function is built into a package and sourced from there library(mypackage) using the same script file to make the package levels(r)[[1]] the same line throws an error, as levels(myraster returns NULL If I modify the script to include the raster namespace: raster::levels(r)[[1]] Then I get the error Error: 'levels<-' is not an exported object from 'namespace:raster' The error suggests that you are trying to assign levels, rather than retrieve them. I have also tried just using levels(r) and putting raster as a depends rather than an import in the DESCRIPTION file for the package, this does not solve the error. Any suggestions on how to overcome the problem? many thanks Nevil Amos [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] stars objects and case_when
Thank you for pointing to this possibility, I'll add it to the stars docs somewhere. This works but only slightly adapted, e.g. as out = read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>% slice(band, 1) %>% setNames("x") %>% mutate(x = case_when (x<100 ~ NA_real_, x>100 & x<150 ~ 2)) where: - I uses setNames("x") so that the attribute is renamed to "x", and you can use x in the case_when expressions (rather than the lengthy "L7_ETMs.tif") - I chnaged NA into NA_real_ : the first RHS of the formulas need to be all of the same type; as typeof(NA) is "logical", it breaks on the second RHS which returns a numeric; if you would TRUE or FALSE rather than 2 in the second case_when formula, using NA would be the right thing to do in the first. The examples of case_when document the need for typed NA in RHS: this is intended behavior. On 9/14/20 4:39 PM, Julian M. Burgos wrote: Dear list, I am wondering if there is a way to use logical statements to replace values of a stars object, for example the case_when function or some other "tidyverse friendly" approach. For example, I can do something like this: st1 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>% slice(band, 1) st1[st1 < 100] <- NA st1[st1 > 100 & st1 < 150] <- 2 ... and so for. But I am wondering if there is a way to do this as part of a pipe, looking something like this: st1 <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars")) %>% slice(band, 1) %>% mutate(x <- case_when (x<100 ~ NA, x>100 & x<150 ~ 2)) Any ideas? Takk, Julian -- Julian Mariano Burgos, PhD Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/ Marine and Freshwater Research Institute Botnsjávarsviðs / Demersal Division Fornubúðir 5, IS-220 Hafnarfjörður, Iceland www.hafogvatn.is Sími/Telephone : +354-5752037 Netfang/Email: julian.bur...@hafogvatn.is ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] stars analogous of raster::aggregate
It now also works for multi-band rasters (fixed in sf when installed from github): https://github.com/r-spatial/stars/issues/320 On 9/9/20 1:10 PM, Julian M. Burgos wrote: Thanks Edzer and Micha, this was very helpful. st_warp with use_gdal=TRUE worked like a charm. My best, Julian Edzer Pebesma writes: On 9/8/20 8:47 PM, Micha Silver wrote: On 08/09/2020 19:33, Julian M. Burgos wrote: Dear all, The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10: library(raster) rst1 <- raster(system.file("tif/L7_ETMs.tif", package = "stars")) rst2 <- aggregate(rst1, fact = 10, fun = mean) res(rst1) [1] 28.5 28.5 res(rst2) [1] 285 285 I am trying to do the same thing with a stars object. The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons. I could do something like this: Probably st_warp() is what you want: l7_file = system.file("tif/L7_ETMs.tif", package = "stars") l7 = read_stars(l7_file) l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs = st_crs(l7)) stars::st_dimensions(l7) from to offset delta refsys point values x 1 349 288776 28.5 +proj=utm +zone=25 +south... FALSE NULL [x] y 1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL [y] band1 6 NANA NANA NULL stars::st_dimensions(l7_lowres) from to offset delta refsys point values x 1 111 28877690 +proj=utm +zone=25 +south...NA NULL [x] y 1 112 9120761 -90 +proj=utm +zone=25 +south...NA NULL [y] band1 6 NANA NANA NULL Yes; if you want to get similar behaviour as raster, choose a cell size that is an exact multiple of the origin's cell size, use_gdal = TRUE, and method = "average"; this currently seems to only work for single band rasters; along these lines: library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0 l7_file = system.file("tif/L7_ETMs.tif", package = "stars") l7 = read_stars(l7_file) (dest = st_as_stars(st_bbox(l7), dx = 90, dy = 90)) # stars object with 2 dimensions and 1 attribute # attribute(s): # values # Min. :0 # 1st Qu.:0 # Median :0 # Mean :0 # 3rd Qu.:0 # Max. :0 # dimension(s): # from to offset delta refsys point values x/y # x1 111 28877690 UTM Zone 25, Southern Hem...NA NULL [x] # y1 112 9120761 -90 UTM Zone 25, Southern Hem...NA NULL [y] (l7_lowres = st_warp(src = l7[,,,1], dest = dest, use_gdal = TRUE, method = "average")) # stars object with 2 dimensions and 1 attribute # attribute(s): # file14bb5c7321dc.tif # Min. : 56.81 # 1st Qu.: 68.75 # Median : 79.00 # Mean : 79.25 # 3rd Qu.: 88.62 # Max. :207.44 # dimension(s): # from to offset delta refsys point values x/y # x1 111 28877690 UTM Zone 25, Southern Hem... FALSE NULL [x] # y1 112 9120761 -90 UTM Zone 25, Southern Hem... FALSE NULL [y] -- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918 https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Forcid.org%2F-0002-1128-1325data=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294sdata=K5kMOHf6dka1eNJrVRZsgDe0UmxTEhBRU40C%2FPpRXVE%3Dreserved=0 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://eur01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-geodata=02%7C01%7C%7Cb8b83838082f41ad457c08d854354350%7C8e105b94435e4303a61063620dbe162b%7C0%7C1%7C637351935026744294sdata=KwymhTSAHdcBDVffZHhp2xZ%2BQUzDKh%2B%2BN2DP1yLtGY4%3Dreserved=0 -- Julian Mariano Burgos, PhD Hafrannsóknastofnun, rannsókna- og ráðgjafarstofnun hafs og vatna/ Marine and Freshwater Research Institute Botnsjávarsviðs / Demersal Division Fornubúðir 5, IS-220 Hafnarfjörður, Iceland www.hafogvatn.is Sími/Telephone : +354-5752037 Netfang/Email: julian.bur...@hafogvatn.is -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] stars analogous of raster::aggregate
On 9/8/20 8:47 PM, Micha Silver wrote: On 08/09/2020 19:33, Julian M. Burgos wrote: Dear all, The raster package has the raster::aggregate function that can be used to reduce the resolution of a raster by aggregating cells by a specific factor. For example, this reduces the resolution of the L7_ETMs.tif raster by a factor of 10: library(raster) rst1 <- raster(system.file("tif/L7_ETMs.tif", package = "stars")) rst2 <- aggregate(rst1, fact = 10, fun = mean) res(rst1) [1] 28.5 28.5 res(rst2) [1] 285 285 I am trying to do the same thing with a stars object. The stars package has the stars::aggregate function, but for spatial aggregation it takes an object of class sf or sfc, so it is meant to be used for aggregation over polygons. I could do something like this: Probably st_warp() is what you want: l7_file = system.file("tif/L7_ETMs.tif", package = "stars") l7 = read_stars(l7_file) l7_lowres = st_warp(src = l7, cellsize = c(90, 90), crs = st_crs(l7)) stars::st_dimensions(l7) from to offset delta refsys point values x 1 349 288776 28.5 +proj=utm +zone=25 +south... FALSE NULL [x] y 1 352 9120761 -28.5 +proj=utm +zone=25 +south... FALSE NULL [y] band 1 6 NA NA NA NA NULL stars::st_dimensions(l7_lowres) from to offset delta refsys point values x 1 111 288776 90 +proj=utm +zone=25 +south... NA NULL [x] y 1 112 9120761 -90 +proj=utm +zone=25 +south... NA NULL [y] band 1 6 NA NA NA NA NULL Yes; if you want to get similar behaviour as raster, choose a cell size that is an exact multiple of the origin's cell size, use_gdal = TRUE, and method = "average"; this currently seems to only work for single band rasters; along these lines: library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0 l7_file = system.file("tif/L7_ETMs.tif", package = "stars") l7 = read_stars(l7_file) (dest = st_as_stars(st_bbox(l7), dx = 90, dy = 90)) # stars object with 2 dimensions and 1 attribute # attribute(s): # values # Min. :0 # 1st Qu.:0 # Median :0 # Mean :0 # 3rd Qu.:0 # Max. :0 # dimension(s): # from to offset delta refsys point values x/y # x1 111 28877690 UTM Zone 25, Southern Hem...NA NULL [x] # y1 112 9120761 -90 UTM Zone 25, Southern Hem...NA NULL [y] (l7_lowres = st_warp(src = l7[,,,1], dest = dest, use_gdal = TRUE, method = "average")) # stars object with 2 dimensions and 1 attribute # attribute(s): # file14bb5c7321dc.tif # Min. : 56.81 # 1st Qu.: 68.75 # Median : 79.00 # Mean : 79.25 # 3rd Qu.: 88.62 # Max. :207.44 # dimension(s): # from to offset delta refsys point values x/y # x1 111 28877690 UTM Zone 25, Southern Hem... FALSE NULL [x] # y1 112 9120761 -90 UTM Zone 25, Southern Hem... FALSE NULL [y] -- Micha Silver Ben Gurion Univ. Sde Boker, Remote Sensing Lab cell: +972-523-665918 https://orcid.org/-0002-1128-1325 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Error when saving an sf (data) object to file as a shapefile
On 6/20/20 12:28 PM, Lom Navanyo wrote: > Unknown field name `Act_Depthtwtr': updating a layer with improper field > name(s)? Shapefiles have funny restrictions as to the length and form of variable names, you probably ran into that, as the message suggests. A general advice is to not use shapefiles for writing; there are better alternatives, geopackage comes to mind. Shapefiles make life unnecessary complicated, you've now run into two of the many problems they have. The days of creating shapefiles are over. If someone asks you to create shapfiles, explain them it is a bad idea. -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal release candidate 1.5-9 rev. 1000 ready for testing
On 6/6/20 9:15 PM, Patrick Schratz wrote: > Since I use and contribute to R (6 years for now) I was always on the > side of "help everyone", "things can get better/easier", "don't create > your own island solutions". However, I was hitting some (in my view, not > understandable) "walls" recently and I am not sure if I want to continue > with this attitude. Dear Patrick, We should keep in mind that the system requirements of the R spatial packages is quite likely the most complicated there is on CRAN, that it has a long legacy, and that GDAL/PROJ came very dynamic recently, and that many people (though not enough) are involved. I think that you're learning that working with people is much harder than working with technology. I highly appreciate the contributions you made to CI for R packages [*], and the way you helped making it work for sf; this helps robust development and finding bugs early. However, the moment something breaks in the CI setup, I have no clue what to do and have to ask your help. If I wouldn't get that help a couple of times, I would drop the whole thing and throw it out of the window. For people using sf this is the same: they have no clue how it works; if it doesn't work a couple of times and help doesn't come quickly, they throw it out of the window and look for something else that does work. What you have to realize is that what looks simple for you (e.g. CI) looks incredibly complicated for other people, for the simple reason that they have been doing very different things for the last 6 years, and have neither time nor ambition to make up for that. Convincing others to change something by arguing the change is "simple and makes things better" often does not work because the receiver has a different perception of "simple" and is not convinced that the status quo is not good enough, or is convinced the change brings a lot of work and/or added complexity. I don't think it helps to get emotional publicly in such a case and move to twitter [&], nor to become disappointed in the R project. If you'd move to another project, you'll find yourself again confronted with people, and discover that communication is always difficult. [*] see e.g. https://github.com/ropensci/tic [&] https://twitter.com/pjs_228/status/1269301044481339392 -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal release candidate 1.5-9 rev. 1000 ready for testing
ic because their production system is blocked >> until the new versions are ready (PROJ and GDAL are C++11 or more, and >> take an order of magnitude longer to compile than just a few years >> ago). >> >> So for Windows and macOS, waiting for the CRAN binaries is a >> reasonable choice. >> >> Beyond this, we need to find ways of providing share/proj and >> share/gdal metadata files for all of the packages now using the PROJ >> and GDAL libraries, and of navigating the content download network for >> geodetic transformation grids available from PROJ 7. But that is >> another story ... >> >> Roger >> >>> >>> Best, >>> Ista >>> >>>> >>>> Manuel >>>> >>>> El vie., 5 jun. 2020 a las 14:31, Patrick Schratz (< >>>> patrick.schr...@gmail.com>) escribió: >>>> >>>>> I am not sure if the part with >>>>> >>>>> use --with-proj_api="proj_api.h" for deprecated API >>>>> >>>>> Is of much help since c/p won’t work but the text let’s people >>>>> assume that c/p could/should work. >>>>> In fact, a full path to “proj_api.h” is required? >>>>> >>>>> I still do not like this blocker and I still do not know if this >>>>> combination causes serious issues in production or just limits new >>>>> features. >>>>> >>>>> For the time being, using and linking osgeo-gdal (3.0.1) and >>>>> osgeo-proj >>>>> (7.0.1) works and can be used as a workaround until homebrew-core >>>>> formulas catch up. >>>>> >>>>>> checks OK on PROJ 7.0.1 and GDAL 2.2.4 >>>>> >>>>> Again, since it was maybe caused by my typo a few mails ago: The >>>>> homebred-core gdal version is at 2.4.4 and not 2.2.4. >>>>> >>>>> On 5 Jun 2020, at 13:29, Roger Bivand wrote: >>>>> >>>>>> The release candidate of rgdal_1.5-9 is ready for testing on >>>>>> R-forge: >>>>>> >>>>>> https://r-forge.r-project.org/R/?group_id=884 >>>>>> >>>>>> Those insisting on installing on PROJ >= 6 and GDAL < 3 must use >>>>>> configure argument --with-proj_api="proj_api.h"; with this used, >>>>>> this >>>>>> version builds with --no-build-vignettes and installs and checks >>>>>> OK on >>>>>> PROJ 7.0.1 and GDAL 2.2.4 with --with-proj_api="proj_api.h". >>>>>> >>>>>> Otherwise checked OK with PROJ 4.8.0, 4.9.2, 4.9.3 and 5.2.0 with >>>>>> GDAL >>>>>> 1.11.4; with PROJ 5.2.0 and GDAL 2.2.4, 2.3.2 and 2.4.2; with PROJ >>>>>> 6.3.2 and GDAL 3.0.4; with PROJ 7.0.1 and GDAL 3.0.4 and 3.1.0. >>>>>> >>>>>> All who have indicated issues with source installs are asked to >>>>>> try >>>>>> the release candidate and to report back here by midnight CEST >>>>>> Monday >>>>>> 8 June. If no indications are forthcoming, I'll assume that >>>>>> problems >>>>>> with 1.5-8 are resolved. >>>>>> >>>>>> Roger >>>>>> >>>>>> -- >>>>>> Roger Bivand >>>>>> Department of Economics, Norwegian School of Economics, >>>>>> Helleveien 30, N-5045 Bergen, Norway. >>>>>> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no >>>>>> https://orcid.org/-0003-2392-6140 >>>>>> https://scholar.google.no/citations?user=AWeghB0J=en >>>>>> >>>>>> ___ >>>>>> R-sig-Geo mailing list >>>>>> R-sig-Geo@r-project.org >>>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>>> ___ >>>>> R-sig-Geo mailing list >>>>> R-sig-Geo@r-project.org >>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>>>> >>>> >>>> >>>> -- >>>> *Manuel Spínola, Ph.D.* >>>> Instituto Internacional en Conservación y Manejo de Vida Silvestre >>>> Universidad Nacional >>>> Apartado 1350-3000 >>>> Heredia >>>> COSTA RICA >>>> mspin...@una.cr >>>> mspinol...@gmail.com >>>> Teléfono: (506) 8706 - 4662 >>>> Personal website: Lobito de río >>>> <https://sites.google.com/site/lobitoderio/> >>>> Institutional website: ICOMVIS <http://www.icomvis.una.ac.cr/> >>>> >>>> [[alternative HTML version deleted]] >>>> >>>> ___ >>>> R-sig-Geo mailing list >>>> R-sig-Geo@r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> -- >> Roger Bivand >> Department of Economics, Norwegian School of Economics, >> Helleveien 30, N-5045 Bergen, Norway. >> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no >> https://orcid.org/-0003-2392-6140 >> https://scholar.google.no/citations?user=AWeghB0J=en___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Package lwgeom (0.2-4) fails to install on Shinyapps.io with an "upgrade GEOS to 3.6.0 or later" error
R package lwgeom had to upgrade to a newer version of liblwgeom (now a subdirectory in the PostGIS source tree, really) because of the new PROJ versions and deprecation of the old PROJ interface; PostGIS requires the newer GEOS version. I have no control over versions of GEOS that are run on the shinyapps.io machines, and don't know which OS at which version they run. I also have no clue how one could debug this issue. One could try to patch the install file, https://github.com/rstudio/shinyapps-package-dependencies/blob/master/packages/lwgeom/install such that it installs a newer GEOS (maybe binary from another ppa, otherwise from source). On 6/3/20 6:33 PM, Rich Shepard wrote: > On Wed, 3 Jun 2020, Erin Hodgess wrote: > >> Could this possibly be due to the upgrade for GDAL and PROJ libraries, >> please? > > Erin, > > I learned a while ago that when any of proj, geos, and gdal are upgraded I > need to rebuid tools higher in that chain. If you have a new geos version > then a new gdal build is in order. > > If this does not address your issue please excuse my responding. > > Regards, > > Rich > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] R 4.0.0 Windows binaries for devel sf w/GDAL3
Similar to sp and rgdal, upcoming sf CRAN windows binaries will use GDAL 3.0.4, PROJ 6.3.1, and GEOS 3.8.0. You can try a beta-version of sf 0.9-3 now by installing it with: install.packages("sf", repos = "https://edzer.github.io/drat;, type = "win.binary") Comments are welcome here, or as issues on github. -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Convert rectilinear to regular grid in R (stars and raster)
The stars part of this question has also been sent as an issue on github, and that is where it has been answered (partially solved): https://github.com/r-spatial/stars/issues/264 please feel free to send follow-up's here, or there. On 3/23/20 6:45 PM, Frederico Faleiro wrote: > Dear all, > > I am trying convert netCDF data from rectilinear to regular grid in R, but > I am not sure what is the best approach for it. > I would like to convert the grid modifying the original resolution only > when necessary. In some cases I must change resolution to keep regular cell > size and cover all the world extent. > I have tried without success st_warp from stars and I find a way with > rasterize from raster package (see script bellow). > > I have the following questions: > 1. What am I doing wrong in stars package (see script bellow)? > 2. Is there a way to rasterize points to raster using metrics that consider > distance (e.g. bilinear, nearest neighbor, etc)? Note that I need > use points in raster package because it do not handle with rectilinear grid. > 3. Do you recommend any other approach? > > The data are available in: > GDrive: > https://drive.google.com/file/d/1HKxFpAwvY9k_cFbasagKHAwi1luzLv90/view?usp=sharing > CMIP5 portal (you must make a login): > http://aims3.llnl.gov/thredds/fileServer/cmip5_css01_data/cmip5/output1/BCC/bcc-csm1-1/rcp45/mon/atmos/Amon/r1i1p1/v20120705/pr/pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc > > # script > library(stars) > library(raster) > > # prepare data > bc <- read_ncdf("pr_Amon_bcc-csm1-1_rcp45_r1i1p1_200601-209912.nc") > # get coordinates and data from the first time period > bc1 <- as.data.frame(bc[, , , 1]) > # convert longitude to variate between -180 and 180 > bc1$lon <- ifelse(bc1$lon > 180, bc1$lon - 360, bc1$lon) > > # original resolution > bc1$lon %>% diff %>% table > bc1$lat %>% diff %>% table > # new resolution > new_res <- c(2.8125, 2.8125) > # Is it multiple of 180 or 360? > 180 %% new_res > > # create a new standard grid > wgs84 <- "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs" > grid_r <- raster(res = new_res, crs = wgs84) > grid_st <- st_as_stars(grid_r) > > # change from rectilinear to regular grid > # stars package: use st_warp to change from rectilinear to regular raster > with different methods > bc_reg <- bc %>% st_warp(grid_st, method = "average") > bc_reg <- bc %>% st_warp(grid_st, method = "near") > bc_reg <- bc %>% st_warp(grid_st, method = "bilinear") > # In all cases I have this error: Error in 1:prod(dims[dxy]) : NA/NaN > argument > > # raster package: rasterize the points values to the new regular grid > bc1_sp <- bc1 > coordinates(bc1_sp) <- ~ lon + lat > bc1_reg <- rasterize(bc1_sp, grid_r, field = "pr", fun = mean) > > # plot the first time period > dev.new(title = "rectilinear") > plot(bc[, , , 1]) > dev.new(title = "regular with raster package") > plot(bc1_mn_reg) > > Best regards, > > Frederico Faleiro > Postdoctoral Researcher in the INCT-EECBio (https://www.eecbio.ufg.br/) > Federal University of Goiás | Brazil > RG: https://www.researchgate.net/profile/Frederico_Faleiro > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] issue when opening raster image in R, the pixel values a rescaled somehow
I see at the end of the output obtained with running gdalinfo on the file: Band 1 Block=1200x3 Type=UInt16, ColorInterp=Gray NoData Value=0 Unit Type: K Offset: 0, Scale:0.02 Metadata: DESCRIPTION=8-day daytime 1km grid Land-surface Temperature This means that Kelving values are obtained by computing 0 (offset) + 0.02 (scale) * pixel_value packages stars and readGDAL do this transformation for you, so you get actual Kelvin values, and represent them with doubles so there is no rounding problem. Other, lower level tools may give you INT2 values where you need to take care of the offset & scale yourself. On 3/24/20 10:28 AM, Roger Bivand wrote: > Opening in stars shows that the read data are interpreted as degrees > Kelvin: > >> library(stars) >> ss <- read_stars("MOD11A2.A249.h17v05.006.2015058135048.tif") >> print(ss, n=1e8) > stars object with 2 dimensions and 1 attribute > attribute(s): > MOD11A2.A249.h17v05.006.2015058135048.tif [K] > Min. :272.1 > 1st Qu.:293.3 > Median :296.5 > Mean :296.1 > 3rd Qu.:298.8 > Max. :309.2 > NA's :926190 > dimension(s): > from to offset delta refsys point values > x 1 1200 -951 926.625 unnamed FALSE NULL [x] > y 1 1200 4447802 -926.625 unnamed FALSE NULL [y] > > The summary agrees with > >> library(rgdal) >> o <- readGDAL("MOD11A2.A249.h17v05.006.2015058135048.tif") >> summary(o) > Object of class SpatialGridDataFrame > Coordinates: > min max > x -951 0 > y 3335852 4447802 > Is projected: TRUE > proj4string : > [+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs] > Grid attributes: > cellcentre.offset cellsize cells.dim > x -487 926.6254 1200 > y 3336315 926.6254 1200 > Data attributes: > band1 > Min. :272.1 > 1st Qu.:293.3 > Median :296.5 > Mean :296.1 > 3rd Qu.:298.8 > Max. :309.2 > NA's :926190 > > and > >> library(raster) >> r <- raster("MOD11A2.A249.h17v05.006.2015058135048.tif") >> summary(values(r)) > Min. 1st Qu. Median Mean 3rd Qu. Max. NA's > 272.1 293.3 296.5 296.1 298.8 309.2 926190 > > So my feeling is that the representation is correct in degrees Kelvin. > > Others please join in, I'm not familiar with MODIS HDF data. I think > something gives the metadata in the GeoTIFF the wrong data bounds - they > do seem to be UInt16. The output of gdalinfo (command line) includes: > > ... > units=K > valid_range=7500, 65535 > ... > Band 1 Block=1200x3 Type=UInt16, ColorInterp=Gray > NoData Value=0 > Unit Type: K > ... > > Roger > > > On Tue, 24 Mar 2020, Victor F. Rodriguez Galiano wrote: > >> Hi Roger, the link to he image is given below. This image was generated >> from a .hdf file using the function gdal_translate. >> >> gdal_translate(sds[1], dst_dataset=paste0(substr(filename, 1, >> nchar(filename)-4) ,".tif", )) >> >> https://wetransfer.com/downloads/4569e9a85c02faad5016b83d0b50217820200324084015/f2eafdfafaae74de9ac8e25017fc5cc520200324084015/e46e80 >> >> >> >> Thanks >> >> Victor >> >> El lun., 23 mar. 2020 a las 21:07, Roger Bivand () >> escribió: >> >>> Please make the file in question available for download (url), or the >>> url >>> you downloaded it from (no login), on this list. It looks as though the >>> unsigned 16 bit integer is not right. >>> >>> Roger Bivand >>> Norwegian School of Economics >>> Bergen, Norway >>> >>> -- >>> *Fra:* R-sig-Geo på vegne av Victor F. >>> Rodriguez Galiano >>> *Sendt:* mandag 23. mars 2020, 21.02 >>> *Til:* r-sig-geo@r-project.org >>> *Emne:* [R-sig-Geo] issue when opening raster image in R, the pixel >>> values a rescaled somehow >>> >>> I opened a Geotiff image in R using the raster function. However, the >>> image seems to be rescaled. The minimum and maximum values should be >>> 13607 >>> and 15461 but are 275 and 305. The Geotiff image when opened in a GIS is >>> correct but not in R. This is my code: Script: library(raster) >>> trial<-raster("MOD11A2.A249.h17v05.006.2015058135048.tif", >>> datatype = >>> "INT2U") trial plot(trial) [[alternative HTML version deleted]] >>> ___ R-sig-Geo mailing list >>> R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> > > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal travis error: configure: error: gdal-config not found or not executable.
Maybe this: https://github.com/edzer/sf_dep/blob/master/.travis.yml On 3/19/20 8:26 PM, Joe Lewis wrote: > Hi, > > Thanks for the help. > > I've incorporated some from sf's travis file, but it's still failing > unfortunately. > > # R for travis: see documentation at > https://docs.travis-ci.com/user/languages/r > > language: R > sudo: false > cache: packages > > r: > - oldrel > - release > - devel > > repos: > CRAN: https://cran.rstudio.com > rforge: http://R-Forge.R-project.org > > before_install: > - R -e 'install.packages("rgdal", repos="http://R-Forge.R-project.org;)' > > apt_packages: > - gdal-bin > - libgdal-dev > - libxml2-dev > - libproj-dev > - libprotobuf-dev > - protobuf-compiler > - libv8-3.14-dev > - libjq-dev > - libudunits2-dev > - libproj-dev > - libgeos-dev > - libspatialite-dev > - libgdal-dev > > - libjson-c-dev > > > https://travis-ci.org/github/josephlewis/leastcostpath/builds/664541703 > > Any ideas? > > Thanks. > > Kind regards, > Joseph Lewis > > On Wed, Mar 18, 2020 at 4:04 PM Edzer Pebesma > wrote: > >> At least libgdal-dev and libproj-dev are missing. You can have a look at >> https://github.com/r-spatial/sf/blob/master/.travis.yml (but that tries >> to do quite a bit more). >> >> On 3/18/20 3:21 PM, Joe Lewis wrote: >>> Hi, >>> >>> I'm testing a package using travis and I keep getting an error message >>> about rgdal not being found or executable. I've been trying to piece >>> together things I've found online but to no avail. >>> >>> https://travis-ci.org/github/josephlewis/leastcostpath/builds/663944787 >>> >>> travis.yml >>> >>> # R for travis: see documentation at >>> https://docs.travis-ci.com/user/languages/r >>> >>> language: R >>> sudo: false >>> cache: packages >>> >>> repos: >>> CRAN: https://cran.rstudio.com >>> rforge: http://R-Forge.R-project.org >>> >>> apt_packages: >>> - gdal-bin >>> >>> before_install: >>> install.packages("rgdal", repos="http://R-Forge.R-project.org;) >>> >>> >>> Any help would be appreciated. >>> >>> Thanks. >>> >>> Kind regards, >>> Joseph Lewis >>> >>> [[alternative HTML version deleted]] >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> -- >> Edzer Pebesma >> Institute for Geoinformatics >> Heisenbergstrasse 2, 48149 Muenster, Germany >> Phone: +49 251 8333081 >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Extract VIIRS data in R
If the product info is in this document: https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/viirs/NASAVIIRSL1BUGAug2019.pdf then it sounds like these data are distributed as curvilinear grids, meaning that there are two additional raster with for each pixel the longitude and the latitude. The file you pointed to does not contain these two grids, but you should be able to get them at the place where you got the data. With that, you can use package stars to load them, and e.g. resample to some (more) regular grid. See e.g. https://r-spatial.github.io/stars/articles/stars4.html#curvilinear-grids On 3/18/20 3:29 PM, Michael Sumner wrote: > You won't find much on this kind of file in the R community, this is > relatively new NetCDF-4 with "groups", and contains a satellite "granule" - > a kind of pre-mapping raw data array with only technical information about > where the scanning occurred. > > In short, raster and stars package (and RNetCDF and ncdf4) will read the > data from this file, but give very little or no help for how to treat it > like a map. The spatial resolution and extent are purely nominal as shown > by raster, a spacing of 1 in two dimensions and extending from 0 to nrows, > 0 to ncols (with a half cell shift). There aren't any coordinate arrays in > the file, or any simple way to transform from this matrix-space to > geography as far as I know. (Would love to be corrected on how to do > that). > > You'd need to pursue domain-specific expertise to discover the mapping of > this, unless some kind soul turns up to help here. I expect you'll want to > find a product that has been converted into simpler form for these data. > > Cheers, Mike. > > > On Wed, Mar 18, 2020 at 1:39 AM Miluji Sb wrote: > >> Dear all, >> >> Hope everyone is keeping safe. >> >> I am trying to extract/read VIIRS nighttime lights data but the output >> seems rather strange. >> >> I have uploaded the file here >> <https://drive.google.com/open?id=1zpqXJ8AlEcnk6ApRb75tfQc4-Y8AfK5h> so as >> not exceed the size limit of the email. >> >> ## Code ## >> library(ncdf4) >> library(rgdal) >> >> netcdf_file <- c("VNP02DNB_NRT.A2020069.1048.001.nc") >> nl <- brick(netcdf_file, lvar=0, values=TRUE, >> varname="observation_data/DNB_observations") >> >> ## Detail ## >> class : RasterLayer >> dimensions : 3232, 4064, 13134848 (nrow, ncol, ncell) >> resolution : 1, 1 (x, y) >> extent : 0.5, 4064.5, 0.5, 3232.5 (xmin, xmax, ymin, ymax) >> crs: NA >> source : C:/Users/shour/Desktop/VNP02DNB_NRT.A2020069.1048.001.nc >> names : DNB.observations.at.pixel.locations >> zvar : observation_data/DNB_observations >> >> The spatial resolution is supposed to be 750m but this shows 1°. What am I >> doing wrong? Any help will be greatly appreciated. Thank you. >> >> Sincerely, >> >> Millu >> >> [[alternative HTML version deleted]] >> >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal travis error: configure: error: gdal-config not found or not executable.
At least libgdal-dev and libproj-dev are missing. You can have a look at https://github.com/r-spatial/sf/blob/master/.travis.yml (but that tries to do quite a bit more). On 3/18/20 3:21 PM, Joe Lewis wrote: > Hi, > > I'm testing a package using travis and I keep getting an error message > about rgdal not being found or executable. I've been trying to piece > together things I've found online but to no avail. > > https://travis-ci.org/github/josephlewis/leastcostpath/builds/663944787 > > travis.yml > > # R for travis: see documentation at > https://docs.travis-ci.com/user/languages/r > > language: R > sudo: false > cache: packages > > repos: > CRAN: https://cran.rstudio.com > rforge: http://R-Forge.R-project.org > > apt_packages: > - gdal-bin > > before_install: > install.packages("rgdal", repos="http://R-Forge.R-project.org;) > > > Any help would be appreciated. > > Thanks. > > Kind regards, > Joseph Lewis > > [[alternative HTML version deleted]] > > _______ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] using the same legend for different scenarios in sf plots
With the breaks argument you can specify the numbers at which colors break. On 2/19/20 3:05 PM, Paulo Flores Ribeiro wrote: > Hello, > I want to compare the spatial distribution of a variable in 7 different > scenarios by running the same plot code successively 7 times (within sf > package), each time changing the values of the variable for each > scenario. In each run, I copy-paste the map into a word document. In the > end, I want to visually compare the 7 maps, but for that I need the > scale range of the legend values to be the same on the seven maps > (including the colour ramp). > How to force the scale interval (and the corresponding colour ramp) to > be the same on the seven maps, within the sf package? > Thanks, > PauloFR > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48149 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] how to find projections supported by sf::st_crs()?
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,]00 > >> 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 > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] PPA ubuntugis-unstable update to Proj6/GDAL3 breaks SF / rgdal
Yes, I had that too. It turned out that my distribution wouldn't update gdal all the way due to all kind of dependencies, and kept it back at some version, resulting in this error. After installing it manually with sudo apt-get install gdal-dev gdal3 got installed. gdal3 does not use pcs.csv. On 12/12/19 5:58 AM, Edouard LEGOUPIL wrote: > Hi, > > Has anyone found how to get sf/rgdal to work with the recent push on > ubuntugis-unstable PPA of the last version of Proj6 & Gdal3 ? > > It seems that sf is looking for pcs.csv that is not here anymore in gdal3 - > configure: error: pcs.csv not found in GDAL data directory. > and rgdal looks for proj_def.dat wich does not exist anymore in proj6 > > Thanks, > Edouard > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] PROJ6 rgdal::project/sp::spTransform - issue of points that cannot be inverted
I don't know why Roger wants to deprecate rgdal::project, but for various reasons I implemented a similar function in sf: library(sf) # Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1 ll <- cbind(rep(180, 3), c(-89, -90, -89)) xy = sf_project("+proj=longlat", "+proj=moll", ll) sf_project("+proj=moll", "+proj=longlat", xy) # [,1] [,2] # [1,] 180 -89 # [2,] Inf Inf # [3,] 180 -89 # Warning message: # In CPL_proj_direct(as.character(c(from[1], to[1])), as.matrix(pts)) : # one or more projected point(s) not finite As it stands now, I want to keep this function in sf, and would be happy to add an argument to suppress the warning msg. On 12/6/19 3:25 PM, Daniel Kelley wrote: > From prior discussions on this thread and elsewhere, I am given to understand > that rgdal::project() will not be carried over into PROJ6, and that we are > recommendd to use sp::spTransform() instead. > > Accordingly, I started work on switching the 'oce' package from > rgdal::project() to sp::spTransform(). However, I have a problem with points > that cannot be inverted. With rgdal::project(), we get warning messages if > the input contains such points (and the output is Inf for the associated > projected points). However, it seems that with sp::spTransform(), we get an > error message and no results, if any of the data contain points that cannot > be inverted. I want to be able to transform a lot of points at the same time > (owing to data-set size and speed considerations), so handling the data > point-by-point is not an option. > > My question is simple: is there a way I can make sp::spTransform() return an > equivalent to that of rgdal::project(), with finite data where possible and > Inf (or similar) where an inverse cannot be done? > > Possibly I am just something. For anyone who has the patience to look, I > have some thoughts at > https://github.com/dankelley/oce/issues/1599#issuecomment-562578641 but the > gist is as follows: note that 'PROJ' consists of three points, one of which > is (of course) Inf, but that 'TRAN' consists only of an error message. Maybe > there is an argument to sp::spTransform() that I'm unaware of, that will > cause it to return data when it can, and Inf when it cannot? > > ``` > R Under development (unstable) (2019-11-07 r77386) -- "Unsuffered > Consequences" > Copyright (C) 2019 The R Foundation for Statistical Computing > Platform: x86_64-apple-darwin15.6.0 (64-bit) > > R is free software and comes with ABSOLUTELY NO WARRANTY. > You are welcome to redistribute it under certain conditions. > Type 'license()' or 'licence()' for distribution details. > > Natural language support but running in an English locale > > R is a collaborative project with many contributors. > Type 'contributors()' for more information and > 'citation()' on how to cite R or R packages in publications. > > Type 'demo()' for some demos, 'help()' for on-line help, or > 'help.start()' for an HTML browser interface to help. > Type 'q()' to quit R. > >> library("rgdal") >> library("sp") >> ll <- cbind(rep(180, 3), c(-89, -90, -89)) >> lonlat <- sp::SpatialPoints(ll, sp::CRS("+proj=longlat")) >> xy <- sp::spTransform(lonlat, sp::CRS("+proj=moll")) >> # rgdal::project returns a mixture of results and Inf values >> dput(rgdal::project(sp::coordinates(xy), proj="+proj=moll", inv=TRUE)) > structure(c(179.9998, Inf, 179.9998, -89.1, > Inf, -89.1), .Dim = 3:2, .Dimnames = list(NULL, c("coords.x1", > "coords.x2"))) >> ## but sp::spTransform returns no data at all >> dput(try(sp::spTransform(xy, sp::CRS("+proj=longlat")), silent=TRUE)) > non finite transformation detected: > coords.x1 coords.x2 > 1.104637e-09 -9.020048e+06 Inf Inf > structure("Error in sp::spTransform(xy, sp::CRS(\"+proj=longlat\")) : \n > failure in points 2\n", class = "try-error", condition = structure(list( > message = "failure in points 2", call = sp::spTransform(xy, > sp::CRS("+proj=longlat"))), class = c("simpleError", > "error", "condition"))) >> > ``` > > > > Dan Kelley > Department of Oceanography > Dalhousie University > Halifax, NS, Canada > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Installation of proj4 v6 on Ubuntu 18.04 Error: proj/epsg not found
Would you mind giving the ouptut of: pkg-config proj --libs pkg-config proj --cflags and also try: install.packages("rgdal", repos="http://R-Forge.R-project.org;) On 10/29/19 7:11 AM, Tomislav Hengl wrote: > > Thanks Edzer, > > I removed my gdal (which removes also all QGIS and everything connect > unfortunately) then I followed your steps and I get (somehow the proj > version is still 5 and not 6): > >> install.packages("rgdal") > % Total % Received % Xferd Average Speed Time Time Time > Current > Dload Upload Total Spent Left > Speed > 100 1631k 100 1631k 0 0 713k 0 0:00:02 0:00:02 --:--:-- > 713k > * installing *source* package ‘rgdal’ ... > ** package ‘rgdal’ successfully unpacked and MD5 sums checked > checking for g++... g++ > checking whether the C++ compiler works... yes > checking for C++ compiler default output file name... a.out > checking for suffix of executables... > checking whether we are cross compiling... no > checking for suffix of object files... o > checking whether we are using the GNU C++ compiler... yes > checking whether g++ accepts -g... yes > configure: CC: gcc -std=gnu99 > configure: CXX: g++ > configure: rgdal: 1.3-3 > checking for /usr/bin/svnversion... yes > configure: svn revision: 759 > checking whether g++ supports C++11 features by default... yes > configure: C++11 support available > checking for gdal-config... /usr/bin/gdal-config > checking gdal-config usability... yes > configure: GDAL: 2.4.2 > checking C++11 support for GDAL >= 2.3.0... yes > checking GDAL version >= 1.11.4... yes > checking gdal: linking with --libs only... yes > checking GDAL: /usr/share/gdal/pcs.csv readable... yes > configure: pkg-config proj exists, will use it > configure: PROJ version: 5.2.0 > checking proj_api.h presence and usability... no > configure: error: proj_api.h not found in standard or given locations. > ERROR: configuration failed for package ‘rgdal’ > * removing ‘/opt/microsoft/ropen/3.5.1/lib64/R/library/rgdal’ > * restoring previous ‘/opt/microsoft/ropen/3.5.1/lib64/R/library/rgdal’ > > The downloaded source packages are in > ‘/data/RTMP/Rtmpl5wj21/downloaded_packages’ > Updating HTML index of packages in '.Library' > Making 'packages.html' ... done > Warning message: > In install.packages("rgdal") : > installation of package ‘rgdal’ had non-zero exit status > > what am I doing wrong? > > > On 10/28/19 10:36 PM, Edzer Pebesma wrote: >> Hi Tom, >> >> formerly PROJ.4 is now called PROJ. >> >> Have you also tried the usual ubuntu install suggestions we give to >> R/ubuntu users, e.g. on https://github.com/r-spatial/sf ? >> >> They are: >> >> sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable >> sudo apt-get update >> sudo apt-get install libudunits2-dev libgdal-dev libgeos-dev libproj-dev >> >> On 10/28/19 4:34 PM, Tomislav Hengl wrote: >>> >>> Dear R-sig-geo, >>> >>> I am trying to install proj4 v6 on Ubuntu 18.04 bionic. I though I can >>> install it simply by following the ubuntugis steps >>> (https://mothergeo-py.readthedocs.io/en/latest/development/how-to/gdal-ubuntu-pkg.html) >>> >>> but this does not lead to newest version of proj4 v6. >>> >>> Then I have tried to follow the proj6 steps >>> (https://proj.org/install.html#install) which indicate that I should do >>> it via the Anaconda installer. After I installed proj4 v6 from conda >>> I get: >>> >>> $ which proj >>> /home/tomislav/anaconda2/bin/proj >>> >>> $ proj -V >>> Rel. 6.2.0, September 1st, 2019 >>> : >>> projection initialization failure >>> cause: no arguments in initialization list >>> program abnormally terminated >>> >>> After that I've added the 'PROJ_LIB = /home/tomislav/anaconda2/bin/proj >>> ' to my .Renviron file and then I get: >>> >>> install.packages("rgdal") >>> Installing package into >>> ‘/home/tomislav/R/x86_64-pc-linux-gnu-library/3.5’ >>> (as ‘lib’ is unspecified) >>> ... >>> checking PROJ.4: epsg found and readable... no >>> Error: proj/epsg not found >>> Either install missing proj support files, for example >>> the proj-nad and proj-epsg RPMs on systems using RPMs, >>> or if installed but not autodetected, set PROJ_LIB to the >>> correct path, and if need be use the --with-proj-share= >>> configure argument. >>> ERROR: configuration fa
Re: [R-sig-Geo] Installation of proj4 v6 on Ubuntu 18.04 Error: proj/epsg not found
Hi Tom, formerly PROJ.4 is now called PROJ. Have you also tried the usual ubuntu install suggestions we give to R/ubuntu users, e.g. on https://github.com/r-spatial/sf ? They are: sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable sudo apt-get update sudo apt-get install libudunits2-dev libgdal-dev libgeos-dev libproj-dev On 10/28/19 4:34 PM, Tomislav Hengl wrote: > > Dear R-sig-geo, > > I am trying to install proj4 v6 on Ubuntu 18.04 bionic. I though I can > install it simply by following the ubuntugis steps > (https://mothergeo-py.readthedocs.io/en/latest/development/how-to/gdal-ubuntu-pkg.html) > but this does not lead to newest version of proj4 v6. > > Then I have tried to follow the proj6 steps > (https://proj.org/install.html#install) which indicate that I should do > it via the Anaconda installer. After I installed proj4 v6 from conda I get: > > $ which proj > /home/tomislav/anaconda2/bin/proj > > $ proj -V > Rel. 6.2.0, September 1st, 2019 > : > projection initialization failure > cause: no arguments in initialization list > program abnormally terminated > > After that I've added the 'PROJ_LIB = /home/tomislav/anaconda2/bin/proj > ' to my .Renviron file and then I get: > > install.packages("rgdal") > Installing package into ‘/home/tomislav/R/x86_64-pc-linux-gnu-library/3.5’ > (as ‘lib’ is unspecified) > ... > checking PROJ.4: epsg found and readable... no > Error: proj/epsg not found > Either install missing proj support files, for example > the proj-nad and proj-epsg RPMs on systems using RPMs, > or if installed but not autodetected, set PROJ_LIB to the > correct path, and if need be use the --with-proj-share= > configure argument. > ERROR: configuration failed for package ‘rgdal’ > * removing ‘/home/tomislav/R/x86_64-pc-linux-gnu-library/3.5/rgdal’ > Warning in install.packages : > installation of package ‘rgdal’ had non-zero exit status > > $gdalinfo --version > GDAL 2.4.0, released 2018/12/14 > > What am I doing wrong? Which tutorial / documentation do you advise to > install proj4 v6 and rgdal package on Ubuntu? > > thank you (especially you Roger!), > > T. Hengl > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Implementing a rolling window for stars object
On 10/26/19 4:43 AM, Micha Silver wrote: > > On 25/10/2019 15:55, Edzer Pebesma wrote: >> Thanks for reporting Micha, and creating a reprex on github Andy; the >> problem seems to be fixed, now. > > > Great, thanks for the quick response. > > I guess reinstalling stars from CRAN should get the fix? Sorry, I meant no: use remotes::install_github("r-spatial/stars") > > >> Best regards, >> >> On 10/23/19 2:22 PM, Andy Teucher wrote: >>> I’m fairly certain this is a bug in stars; I opened an issue here: >>> https://github.com/r-spatial/stars/issues/223 >>> <https://github.com/r-spatial/stars/issues/223> >>> >>> Cheers, >>> Andy Teucher >>> >>>> On Oct 23, 2019, at 10:42 AM, Andy Teucher >>>> wrote: >>>> >>>> Interesting - there’s some rlang/tidy evaluation trickery going on >>>> there that I couldn’t quite figure out (I think it might be >>>> searching for yr in the wrong environment), but defining your range >>>> as a single variable, and putting that in the square brackets seems >>>> to work for me: >>>> >>>> rng <- yr:last_yr >>>> stars_window = ci_stars[,,,rng] >>>> >>>> >>>>> On Oct 23, 2019, at 10:20 AM, Micha Silver wrote: >>>>> >>>>> >>>>> On 23/10/2019 19:30, Andy Teucher wrote: >>>>>> Hi Micha, >>>>>> >>>>>> I can see two problems immediately with your code: >>>>>> 1. you are using a double-colon (yr::last_yr) - the double colon >>>>>> is used for looking for an object in a package, so it is looking >>>>>> for object ‘yrs’ in package ‘yr’, which obviously doesn’t make >>>>>> sense. Use a single colon to create a range (like you did with 2:6) >>>>>> 2. the object ‘last_yr’ is never defined, so even if you used a >>>>>> single colon to define the range yr:last_yr, it would fail as it >>>>>> would not be able to find object ‘last_yr’ >>>>>> >>>>> Thanks, >>>>> >>>>> I fixed those typos (corrected script attached) and I still get >>>>> this error: >>>>> >>>>> >>>>> micha@tp480:R$ Rscript stars_window.R >>>>> Loading required package: abind >>>>> Loading required package: sf >>>>> Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0 >>>>> Error in eval(rlang::expr(x[[i]][!!!args])) : object 'yr' not found >>>>> Calls: RunMK -> [ -> [.stars -> structure -> eval -> eval >>>>> Execution halted >>>>> >>>>> >>>>> >>>>>> Cheers, >>>>>> Andy Teucher >>>>>> >>>>>>> On Oct 23, 2019, at 8:28 AM, Micha Silver >>>>>> <mailto:tsvi...@gmail.com>> wrote: >>>>>>> >>>>>>> I am trying to run a function (mk.test to find MannKendall >>>>>>> trends) using st_apply over a "rolling" window for a time series >>>>>>> of rasters in a stars object. >>>>>>> When I use subscript notation to slice out the window dimension >>>>>>> with a looping variable I get an error: >>>>>>> >>>>>>> Error in loadNamespace(name) : there is no package called ‘yr’ >>>>>>> Calls: [ ... tryCatch -> tryCatchList -> tryCatchOne -> >>>>>>> Execution halted >>>>>>> >>>>>>> However If I replace the subscript with integers it works fine. >>>>>>> (see attached) >>>>>>> What is the correct way to work this out? >>>>>>> >>>>>>> Attached is a reprex with a small subset of my data. (The script >>>>>>> starts with a long structure, code is at the end) >>>>>>> >>>>>>> Thanks >>>>>>> -- >>>>>>> Micha Silver >>>>>>> Ben Gurion Univ. >>>>>>> Sde Boker, Remote Sensing Lab >>>>>>> cell: +972-523-665918 >>>>>>> ___ >>>>>>> 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 >>>>> -- >>>>> Micha Silver >>>>> Ben Gurion Univ. >>>>> Sde Boker, Remote Sensing Lab >>>>> cell: +972-523-665918 >>>>> >>>>> >>> >>> [[alternative HTML version deleted]] >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Implementing a rolling window for stars object
Thanks for reporting Micha, and creating a reprex on github Andy; the problem seems to be fixed, now. Best regards, On 10/23/19 2:22 PM, Andy Teucher wrote: > I’m fairly certain this is a bug in stars; I opened an issue here: > https://github.com/r-spatial/stars/issues/223 > <https://github.com/r-spatial/stars/issues/223> > > Cheers, > Andy Teucher > >> On Oct 23, 2019, at 10:42 AM, Andy Teucher wrote: >> >> Interesting - there’s some rlang/tidy evaluation trickery going on there >> that I couldn’t quite figure out (I think it might be searching for yr in >> the wrong environment), but defining your range as a single variable, and >> putting that in the square brackets seems to work for me: >> >> rng <- yr:last_yr >> stars_window = ci_stars[,,,rng] >> >> >>> On Oct 23, 2019, at 10:20 AM, Micha Silver wrote: >>> >>> >>> On 23/10/2019 19:30, Andy Teucher wrote: >>>> Hi Micha, >>>> >>>> I can see two problems immediately with your code: >>>> 1. you are using a double-colon (yr::last_yr) - the double colon is used >>>> for looking for an object in a package, so it is looking for object ‘yrs’ >>>> in package ‘yr’, which obviously doesn’t make sense. Use a single colon to >>>> create a range (like you did with 2:6) >>>> 2. the object ‘last_yr’ is never defined, so even if you used a single >>>> colon to define the range yr:last_yr, it would fail as it would not be >>>> able to find object ‘last_yr’ >>>> >>> >>> Thanks, >>> >>> I fixed those typos (corrected script attached) and I still get this error: >>> >>> >>> micha@tp480:R$ Rscript stars_window.R >>> Loading required package: abind >>> Loading required package: sf >>> Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0 >>> Error in eval(rlang::expr(x[[i]][!!!args])) : object 'yr' not found >>> Calls: RunMK -> [ -> [.stars -> structure -> eval -> eval >>> Execution halted >>> >>> >>> >>>> Cheers, >>>> Andy Teucher >>>> >>>>> On Oct 23, 2019, at 8:28 AM, Micha Silver >>>> <mailto:tsvi...@gmail.com>> wrote: >>>>> >>>>> I am trying to run a function (mk.test to find MannKendall trends) using >>>>> st_apply over a "rolling" window for a time series of rasters in a stars >>>>> object. >>>>> When I use subscript notation to slice out the window dimension with a >>>>> looping variable I get an error: >>>>> >>>>> Error in loadNamespace(name) : there is no package called ‘yr’ >>>>> Calls: [ ... tryCatch -> tryCatchList -> tryCatchOne -> >>>>> Execution halted >>>>> >>>>> However If I replace the subscript with integers it works fine. (see >>>>> attached) >>>>> What is the correct way to work this out? >>>>> >>>>> Attached is a reprex with a small subset of my data. (The script starts >>>>> with a long structure, code is at the end) >>>>> >>>>> Thanks >>>>> -- >>>>> Micha Silver >>>>> Ben Gurion Univ. >>>>> Sde Boker, Remote Sensing Lab >>>>> cell: +972-523-665918 >>>>> ___ >>>>> 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 >>>> >>> -- >>> Micha Silver >>> Ben Gurion Univ. >>> Sde Boker, Remote Sensing Lab >>> cell: +972-523-665918 >>> >>> >> > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] CSV with Geometry Column to SF object
Without sharing a (minimal) reproducible example, it is unlikely that someone else can help you find out whether this points to a problem in your data, or in the software. On 10/15/19 2:31 PM, argunaw . wrote: > When I run this, I get the following error: > > *Error in CPL_sfc_from_wkt(x) : OGR error* > > When I run it on my data, I get the same error. > > On Tue, Oct 15, 2019 at 2:16 PM Edzer Pebesma > mailto:edzer.pebe...@uni-muenster.de>> > wrote: > > You may try something along these lines: > > # read data.frame with read.csv; here, we create an example by hand: > df = data.frame( > a = 1:3, b = 3:1, geom = c("LINESTRING(0 0, 1 1)", > "LINESTRING(1 1,2 > 2)", "LINESTRING(5 5,6 6)") > ) > > library(sf) > sf = st_sf(df, geom = st_as_sfc(df$geom)) > sf > > > On 10/15/19 1:58 PM, argunaw . wrote: > > I'm not sure how it was exported from postgis- the person who gave > me the > > file wasn't the one who downloaded it unfortunately. > > > > The file is a line file of roads. The files main columns are road ID > > numbers (type integer) and the geometry column (type geometry, > long strong > > of letters and numbers). Only the IT admins where I am have the > postgis > > load/import tools in the pgadmin/sql interface. The rest of us can > download > > from sql and create new tables from other sql databases, but not > create a > > new table from a csv file. > > > > On Tue, Oct 15, 2019 at 1:11 PM Alex Mandel > mailto:tech_...@wildintellect.com>> > > wrote: > > > >> On 10/15/19 8:43 AM, argunaw . wrote: > >>> Hello Everyone, > >>> > >>> I have a csv file with a postgis "geometry" column. I've been > trying to > >>> import it in to R as a SF file, with the goal of exporting it to a > >> postgis > >>> database, but to no avail. I've used the following methods: > >>> > >>> 1. file <- st_read("name.csv", stringsAsFactors=F, > geometry_column=geom) > >>> > >>> 2. file <- fread("name.csv", headers=True) > >>> file <- st_as_sf(file) > >>> > >>> How can I import a csv with a postgis "geometry" column in to R as a > >>> spatial/SF object? > >>> > >>> [[alternative HTML version deleted]] > >>> > >>> ___ > >>> 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 > >>> > >> > >> Can you paste an example somewhere, is it binary data or some kind of > >> plain text column? Do you know how it was exported from Postgis? > >> > >> If it's a dump from a postgis database did you try loading the table > >> directly to postgis with it's own load/import, or sql tools? > >> > >> Thanks, > >> Alex > >> > > > > [[alternative HTML version deleted]] > > > > ___ > > 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 > > > > -- > Edzer Pebesma > Institute for Geoinformatics > Heisenbergstrasse 2, 48151 Muenster, Germany > Phone: +49 251 8333081 > ___ > 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 > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] CSV with Geometry Column to SF object
You may try something along these lines: # read data.frame with read.csv; here, we create an example by hand: df = data.frame( a = 1:3, b = 3:1, geom = c("LINESTRING(0 0, 1 1)", "LINESTRING(1 1,2 2)", "LINESTRING(5 5,6 6)") ) library(sf) sf = st_sf(df, geom = st_as_sfc(df$geom)) sf On 10/15/19 1:58 PM, argunaw . wrote: > I'm not sure how it was exported from postgis- the person who gave me the > file wasn't the one who downloaded it unfortunately. > > The file is a line file of roads. The files main columns are road ID > numbers (type integer) and the geometry column (type geometry, long strong > of letters and numbers). Only the IT admins where I am have the postgis > load/import tools in the pgadmin/sql interface. The rest of us can download > from sql and create new tables from other sql databases, but not create a > new table from a csv file. > > On Tue, Oct 15, 2019 at 1:11 PM Alex Mandel > wrote: > >> On 10/15/19 8:43 AM, argunaw . wrote: >>> Hello Everyone, >>> >>> I have a csv file with a postgis "geometry" column. I've been trying to >>> import it in to R as a SF file, with the goal of exporting it to a >> postgis >>> database, but to no avail. I've used the following methods: >>> >>> 1. file <- st_read("name.csv", stringsAsFactors=F, geometry_column=geom) >>> >>> 2. file <- fread("name.csv", headers=True) >>> file <- st_as_sf(file) >>> >>> How can I import a csv with a postgis "geometry" column in to R as a >>> spatial/SF object? >>> >>> [[alternative HTML version deleted]] >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> >> Can you paste an example somewhere, is it binary data or some kind of >> plain text column? Do you know how it was exported from Postgis? >> >> If it's a dump from a postgis database did you try loading the table >> directly to postgis with it's own load/import, or sql tools? >> >> Thanks, >> Alex >> > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] calculation of kriging range, nugget, sill
If you can read R and C, yes: https://github.com/r-spatial/gstat On 10/7/19 9:08 PM, Daniel Turner via R-sig-Geo wrote: > Hello! > Is there a way to see how the range, nugget, and sill for kriging are > calculated in the gstat package? > Thanks, > Daniel Turner > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] variogramCloud to gstatVariogram in gstat
Not elegant, but this might work: x # is your variogramCloud x$np = 1 class(x) = c("gstatVariogram", "data.frame") On 10/1/19 8:20 PM, Benjamin Hemingway wrote: > Hello, > > I am computing the variogram from points measured along profiles. I do not > want to compute the variogram from pairs of points belonging to the same > profile. To accomplish this I am iterating through each point and computing > the variogram cloud, keeping only the point pairs that I want. However, I > would like to plot the sample variogram and fit a model to it using only > the point pairs in my final variogramCloud data frame. > > Is there an elegant way to convert a variogramCloud data frame to a > gstatVariogram data frame? > > Ben > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal, PROJ6 and "+init=epsg" syntax on Jupyter/conda
CRAN, btw, reports similar errors on all debian testing platforms, of the kind > st_crs("+init=epsg:3857")$epsg proj_create: init=epsg:/init=IGNF: syntax not supported in non-PROJ4 emulation mode see https://cran.r-project.org/web/checks/check_results_edzer.pebesma_at_uni-muenster.de.html#sf I failed to reproduce these because I can't build the CRAN debian testing docker image from https://github.com/jeroen/rcheckserver , reported here: https://github.com/jeroen/rcheckserver/issues/2 I would appreciate any help here. On 9/12/19 11:08 AM, Edzer Pebesma wrote: > > > On 9/12/19 10:55 AM, James Sample wrote: >> Thanks Edzer - that's very helpful! >> >> Would it be possible for you to link/share your Ubuntu Dockerfile, please? >> (I completely understand if you'd rather not, of course). > > https://github.com/r-spatial/sf/tree/master/inst/docker/gdal > >> >> Perhaps if I can see how you're setting thing up in Ubuntu I can modify my >> Jupyter Dockerfile accordingly, and if I can convince myself that this is a >> actually conda-forge/build issue (rather than my own mistake) I can open an >> issue on the conda r-rgdal feedstock ( >> https://github.com/conda-forge/r-rgdal-feedstock). > > Great, please also report back here! > >> >> Thanks again for your help! >> >> Best wishes, >> >> >> James. >> >> >> >> On Wed, 11 Sep 2019 at 23:10, Edzer Pebesma >> wrote: >> >>> On an ubuntu docker image with gdal 3.0.1 and PROJ 6.2.0, I see: >>> >>>> library(rgdal) >>> Loading required package: sp >>> rgdal: version: 1.4-6, (SVN revision (unknown)) >>> Geospatial Data Abstraction Library extensions to R successfully loaded >>> Loaded GDAL runtime: GDAL 3.0.1, released 2019/06/28 >>> Path to GDAL shared files: >>> GDAL binary built with GEOS: FALSE >>> Loaded PROJ.4 runtime: Rel. 6.2.0, September 1st, 2019, [PJ_VERSION: 620] >>> Path to PROJ.4 shared files: (autodetected) >>> Linking to sp version: 1.3-1 >>>> CRS("+init=epsg:3035") >>> CRS arguments: >>> +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 >>> +y_0=321 +ellps=GRS80 +units=m +no_defs >>> >>> but also a warning with sf: >>> >>>> library(sf) >>> Linking to GEOS 3.7.2, GDAL 3.0.1, PROJ 6.2.0 >>>> st_crs("+init=epsg:3035") >>> Coordinate Reference System: >>> No EPSG code >>> proj4string: "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=321 >>> +ellps=GRS80 +units=m +no_defs" >>> Warning message: >>> In CPL_crs_from_proj4string(x) : >>> GDAL Message 1: +init=epsg: syntax is deprecated. It might return >>> a CRS with a non-EPSG compliant axis order. >>> >>> so it feels like a combination of how PROJ has been installed, and how >>> it has been compiled into the R packages. >>> >>> >>> On 9/11/19 9:38 PM, James Sample wrote: >>>> Dear all, >>>> >>>> I am trying to setup R and rgdal within a JupyterLab environment >>> alongside >>>> my usual Python tools. I realise this is probably an unfamiliar setup for >>>> many, but I hope someone might be able to help nevertheless. >>>> >>>> I'm using a Docker container based on Ubuntu 18.04 and derived from the >>>> Jupyter Data Science Notebook ( >>>> >>> https://github.com/jupyter/docker-stacks/tree/master/datascience-notebook >>> ). >>>> I have Python 3.7 and R 3.6 installed, and I'm using "conda" as my >>> package >>>> manager. >>>> >>>> I have successfully installed GDAL and PROJ, together with various Python >>>> and R packages, including 'sp' and 'rgdal'. When I run >>>> >>>> require(rgdal) >>>> >>>> I see the following >>>> >>>> Loading required package: rgdal >>>> Loading required package: sp >>>> rgdal: version: 1.4-4, (SVN revision 833) >>>> Geospatial Data Abstraction Library extensions to R successfully loaded >>>> Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28 >>>> Path to GDAL shared files: /opt/conda/share/gdal >>>> GDAL binary built with GEOS: TRUE >>>> Loaded PROJ.4 runtime: Rel. 6.1.0, May 15th, 2019, [PJ_VERSION: 610] >>>> Path to PROJ.4 shared files: (autodetected) >>>> Linking to sp version: 1.3-1 >>>
Re: [R-sig-Geo] rgdal, PROJ6 and "+init=epsg" syntax on Jupyter/conda
On 9/12/19 10:55 AM, James Sample wrote: > Thanks Edzer - that's very helpful! > > Would it be possible for you to link/share your Ubuntu Dockerfile, please? > (I completely understand if you'd rather not, of course). https://github.com/r-spatial/sf/tree/master/inst/docker/gdal > > Perhaps if I can see how you're setting thing up in Ubuntu I can modify my > Jupyter Dockerfile accordingly, and if I can convince myself that this is a > actually conda-forge/build issue (rather than my own mistake) I can open an > issue on the conda r-rgdal feedstock ( > https://github.com/conda-forge/r-rgdal-feedstock). Great, please also report back here! > > Thanks again for your help! > > Best wishes, > > > James. > > > > On Wed, 11 Sep 2019 at 23:10, Edzer Pebesma > wrote: > >> On an ubuntu docker image with gdal 3.0.1 and PROJ 6.2.0, I see: >> >>> library(rgdal) >> Loading required package: sp >> rgdal: version: 1.4-6, (SVN revision (unknown)) >> Geospatial Data Abstraction Library extensions to R successfully loaded >> Loaded GDAL runtime: GDAL 3.0.1, released 2019/06/28 >> Path to GDAL shared files: >> GDAL binary built with GEOS: FALSE >> Loaded PROJ.4 runtime: Rel. 6.2.0, September 1st, 2019, [PJ_VERSION: 620] >> Path to PROJ.4 shared files: (autodetected) >> Linking to sp version: 1.3-1 >>> CRS("+init=epsg:3035") >> CRS arguments: >> +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 >> +y_0=321 +ellps=GRS80 +units=m +no_defs >> >> but also a warning with sf: >> >>> library(sf) >> Linking to GEOS 3.7.2, GDAL 3.0.1, PROJ 6.2.0 >>> st_crs("+init=epsg:3035") >> Coordinate Reference System: >> No EPSG code >> proj4string: "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=321 >> +ellps=GRS80 +units=m +no_defs" >> Warning message: >> In CPL_crs_from_proj4string(x) : >> GDAL Message 1: +init=epsg: syntax is deprecated. It might return >> a CRS with a non-EPSG compliant axis order. >> >> so it feels like a combination of how PROJ has been installed, and how >> it has been compiled into the R packages. >> >> >> On 9/11/19 9:38 PM, James Sample wrote: >>> Dear all, >>> >>> I am trying to setup R and rgdal within a JupyterLab environment >> alongside >>> my usual Python tools. I realise this is probably an unfamiliar setup for >>> many, but I hope someone might be able to help nevertheless. >>> >>> I'm using a Docker container based on Ubuntu 18.04 and derived from the >>> Jupyter Data Science Notebook ( >>> >> https://github.com/jupyter/docker-stacks/tree/master/datascience-notebook >> ). >>> I have Python 3.7 and R 3.6 installed, and I'm using "conda" as my >> package >>> manager. >>> >>> I have successfully installed GDAL and PROJ, together with various Python >>> and R packages, including 'sp' and 'rgdal'. When I run >>> >>> require(rgdal) >>> >>> I see the following >>> >>> Loading required package: rgdal >>> Loading required package: sp >>> rgdal: version: 1.4-4, (SVN revision 833) >>> Geospatial Data Abstraction Library extensions to R successfully loaded >>> Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28 >>> Path to GDAL shared files: /opt/conda/share/gdal >>> GDAL binary built with GEOS: TRUE >>> Loaded PROJ.4 runtime: Rel. 6.1.0, May 15th, 2019, [PJ_VERSION: 610] >>> Path to PROJ.4 shared files: (autodetected) >>> Linking to sp version: 1.3-1 >>> >>> >>> which seems OK. Most things work as expected, but this >>> >>> CRS("+init=epsg:3035") >>> >>> gives an exception >>> >>> Error in CRS("+init=epsg:3035"): no arguments in initialization list >>> >>> The same code works fine in R-Studio, although I note that my R-Studio >>> installation has PROJ 4.9.2 (with the same versions of rgdal and sp as >>> listed above). Unfortunately I can't downgrade PROJ, since some of my >>> Python packages require version 6.1. >>> >>> I have read that rgdal is compatible with PROJ6 and I haven't been able >> to >>> find (m)any similar issues online, so I assume I'm doing something wrong. >>> As a workaround, this runs successfully >>> >>> showP4(showWKT("+init=epsg:3035")) >>> >>>
Re: [R-sig-Geo] rgdal, PROJ6 and "+init=epsg" syntax on Jupyter/conda
On an ubuntu docker image with gdal 3.0.1 and PROJ 6.2.0, I see: > library(rgdal) Loading required package: sp rgdal: version: 1.4-6, (SVN revision (unknown)) Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 3.0.1, released 2019/06/28 Path to GDAL shared files: GDAL binary built with GEOS: FALSE Loaded PROJ.4 runtime: Rel. 6.2.0, September 1st, 2019, [PJ_VERSION: 620] Path to PROJ.4 shared files: (autodetected) Linking to sp version: 1.3-1 > CRS("+init=epsg:3035") CRS arguments: +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=321 +ellps=GRS80 +units=m +no_defs but also a warning with sf: > library(sf) Linking to GEOS 3.7.2, GDAL 3.0.1, PROJ 6.2.0 > st_crs("+init=epsg:3035") Coordinate Reference System: No EPSG code proj4string: "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=321 +ellps=GRS80 +units=m +no_defs" Warning message: In CPL_crs_from_proj4string(x) : GDAL Message 1: +init=epsg: syntax is deprecated. It might return a CRS with a non-EPSG compliant axis order. so it feels like a combination of how PROJ has been installed, and how it has been compiled into the R packages. On 9/11/19 9:38 PM, James Sample wrote: > Dear all, > > I am trying to setup R and rgdal within a JupyterLab environment alongside > my usual Python tools. I realise this is probably an unfamiliar setup for > many, but I hope someone might be able to help nevertheless. > > I'm using a Docker container based on Ubuntu 18.04 and derived from the > Jupyter Data Science Notebook ( > https://github.com/jupyter/docker-stacks/tree/master/datascience-notebook). > I have Python 3.7 and R 3.6 installed, and I'm using "conda" as my package > manager. > > I have successfully installed GDAL and PROJ, together with various Python > and R packages, including 'sp' and 'rgdal'. When I run > > require(rgdal) > > I see the following > > Loading required package: rgdal > Loading required package: sp > rgdal: version: 1.4-4, (SVN revision 833) > Geospatial Data Abstraction Library extensions to R successfully loaded > Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28 > Path to GDAL shared files: /opt/conda/share/gdal > GDAL binary built with GEOS: TRUE > Loaded PROJ.4 runtime: Rel. 6.1.0, May 15th, 2019, [PJ_VERSION: 610] > Path to PROJ.4 shared files: (autodetected) > Linking to sp version: 1.3-1 > > > which seems OK. Most things work as expected, but this > > CRS("+init=epsg:3035") > > gives an exception > > Error in CRS("+init=epsg:3035"): no arguments in initialization list > > The same code works fine in R-Studio, although I note that my R-Studio > installation has PROJ 4.9.2 (with the same versions of rgdal and sp as > listed above). Unfortunately I can't downgrade PROJ, since some of my > Python packages require version 6.1. > > I have read that rgdal is compatible with PROJ6 and I haven't been able to > find (m)any similar issues online, so I assume I'm doing something wrong. > As a workaround, this runs successfully > > showP4(showWKT("+init=epsg:3035")) > > But the original syntax is cleaner and I'd rather not refactor all my old > code if I can help it. Can anyone point me in the right direction, please? > Is the "+init=epsg" syntax supported with PROJ 6, or should I be using an > alternative? > > Thanks and best wishes, > > > James. > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How can I manipulate the parameters of st_interpolate_aw so it doesn't fail as a result of point and line intersections?
Thanks for the clear diagnosis; a small reprex would be library(sf) a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0 b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5, -0.5,-0.5,1,1 c<-st_sfc(a,b) aa = st_sf(value = 1, geom = c[1]) bb = st_sf(value = 2, geom = c[2]) cc = st_sf(value = 3, geom = c[1] + c(.5,.5)) st_interpolate_aw(aa, cc, extensive = TRUE) cc = st_sf(value = 3, geom = c[1] + c(1,1)) st_interpolate_aw(aa, cc, extensive = TRUE) st_interpolate_aw(aa, bb, extensive = TRUE) revealing two problems: 0 or 1 dimensionsal intersections, and intersections giving GEOMETRYCOLLECTION. This now should work when you install sf from github. On 9/4/19 8:32 AM, Chris Fowler wrote: > I have also asked this question here: > https://stackoverflow.com/questions/57767022/how-do-you-use-st-interpolate-aw-with-polygon-layers-that-legitimately-include-p?noredirect=1#comment101987833_57767022 > > but expect that the specialized knowledge on this list serve may achieve a > better result and that members of this list who are, like me, still making > the transition from sp to sf could benefit from the responses it generates. > > I am trying to use st_interpolate_aw() to assign values from one polygon > layer to another polygon layer (voting district vote totals assigned to > census tracts). The operation fails because st_interpolate_aw() calls > st_intersection() which finds point and line intersections that then get > treated as points, lines, or geometrycollections and break the subsequent > code casting geometries to multipolygons (I get the error "Error in > st_cast.POINT(x[1], to, ...) :cannot create MULTIPOLYGON from POINT") > > Within st_interpolate_aw the relevant code looks like this where 'x' is my > first polygon layer and 'to' is my second: > > i = st_intersection(st_geometry(x), st_geometry(to)) > idx = attr(i, "idx") > i = st_cast(i, "MULTIPOLYGON") > > When I use debug to step through the above lines of code in > st_interpolate_aw I can look at summary(i) after running st_intersection > and I can see multiple POINTS, GEOMETRYCOLLECTIONS and LINESTRINGS that > were created in addition to the desired POLYGONS and MULTIPOLYGONS. This is > not a problem of messy layers per se--we should expect point and line > intersections when comparing these kinds of administrative boundaries, but > they shouldn't be relevant for interpolation purposes and could be > reasonably dropped. The problem is that st_cast doesn't seem set up to make > allowances. > > I expect the answer to this question is to figure out what parameters I can > set via ... in st_interpolate_aw to alter the behavior of st_intersection() > or st_cast() to make this work. Alternatively, there may be a path through > tricks like st_set_precision() or st_snap() to pre-process my polygons so > this doesn't happen (see somewhat relevant back and forth here > <https://github.com/r-spatial/sf/issues/860>). I have already run > st_is_valid() on my polygons and they are fine. Also, I fundamentally think > that st_intersection is behaving appropriately with regards to the data, > and that it is st_cast not being willing to drop geometries with no area > that causes the problem. > > The example below doesn't get me all the way to the solution as it deals > with the code internal to st_interpolate_aw not the function itself, but it > is small, draws on polygons from the vignette > <https://cran.r-project.org/web/packages/sf/vignettes/sf1.html>, and > generates the same error as my code when running st_interpolate_aw(). If > the example below could be run with a parameter change to st_intersection > or st_cast so that st_cast doesn't return an error I > think I could get st_interpolate_aw to work > > a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0 > b <- st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,.5,0,0,0.5, > -0.5,-0.5,1,1 > c<-st_sfc(a,b) > plot(c) > i <- st_intersection(c[[1]],c[[2]]) > i=st_cast(i,"MULTIPOLYGON") > > Thanks, > Chris > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Error in netCDF file: cells are not equally spaced
depend a lot on very specific details. (angstroms >>> package has similar helpers for the approach but is unlikely to help here >>> without access to the file to try) >>> >>> To help us guess further you can use the "ncdump" output, raster will >>> print this with >>> >>> print(raster("yourfile")) >>> >>> and from that we could provide better guesses at variable names and >>> subsetting etc. >>> >>> But, as Edzer asked, nothing is as good as having the file - any chance >>> you can share it? (Personally I think the world rather needs more R >>> attention on climate model output.) >>> >>> Cheers, Mike. >>> >>> On Thu, Aug 22, 2019 at 1:53 PM Frederico Faleiro >>> wrote: >>>> >>>> Dear list, >>>> >>>> I am using some netCDF files from the CMIP5 climate models in raster >>>> package, but I am having some issues with one model. >>>> The netCDF file from the GFDL-ESM2G model (e.g. >>>> >>> http://aims3.llnl.gov/thredds/fileServer/css03_data/cmip5/output1/NOAA-GFDL/GFDL-ESM2G/rcp45/mon/atmos/Amon/r1i1p1/v20120412/pr/pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc >>>> ) >>>> have the error message as in bellow example. >>>> >>>> #Example >>>> library(raster) >>>> s1 <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc") >>>> Error in .rasterObjectFromCDF(x, type = objecttype, band = band, ...) : >>>> cells are not equally spaced; you should extract values as points >>>> >>>> # I check some solutions that recomend force read the file with the >>>> argument "stopIfNotEqualSpaced = F" as bellow. >>>> sf <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc", >>>> *stopIfNotEqualSpaced >>>> = F*) >>>> bf <- brick("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc", >>>> *stopIfNotEqualSpaced >>>> = F*) >>>> >>>> I performed some tests and only "works" with brick. However I did not >>> find >>>> any solution to check where is the problem and fix it in the file. >>>> >>>> How can I check if it is really an issue and fix it? >>>> >>>> sessionInfo() >>>> R version 3.5.1 (2018-07-02) >>>> Platform: x86_64-w64-mingw32/x64 (64-bit) >>>> Running under: Windows 10 x64 (build 18362) >>>> >>>> Matrix products: default >>>> >>>> locale: >>>> [1] LC_COLLATE=Portuguese_Brazil.1252 LC_CTYPE=Portuguese_Brazil.1252 >>>> [3] LC_MONETARY=Portuguese_Brazil.1252 LC_NUMERIC=C >>>> [5] LC_TIME=Portuguese_Brazil.1252 >>>> >>>> attached base packages: >>>> [1] stats graphics grDevices utils datasets methods base >>>> >>>> other attached packages: >>>> [1] raster_2.8-19 sp_1.3-1 >>>> >>>> loaded via a namespace (and not attached): >>>> [1] compiler_3.5.1 rgdal_1.4-3 Rcpp_1.0.1 codetools_0.2-15 >>>> ncdf4_1.16.1 >>>> [6] grid_3.5.1 lattice_0.20-35 >>>> >>>> Best regards, >>>> >>>> Frederico >>>> >>>> [[alternative HTML version deleted]] >>>> >>>> ___ >>>> R-sig-Geo mailing list >>>> R-sig-Geo@r-project.org >>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >>> >>> >>> -- >>> Michael Sumner >>> Software and Database Engineer >>> Australian Antarctic Division >>> Hobart, Australia >>> e-mail: mdsum...@gmail.com >>> >> > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Error in netCDF file: cells are not equally spaced
I would be happy to look into this, but the file you point to requires authentification. On 8/22/19 5:52 AM, Frederico Faleiro wrote: > Dear list, > > I am using some netCDF files from the CMIP5 climate models in raster > package, but I am having some issues with one model. > The netCDF file from the GFDL-ESM2G model (e.g. > http://aims3.llnl.gov/thredds/fileServer/css03_data/cmip5/output1/NOAA-GFDL/GFDL-ESM2G/rcp45/mon/atmos/Amon/r1i1p1/v20120412/pr/pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc > ) > have the error message as in bellow example. > > #Example > library(raster) > s1 <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc") > Error in .rasterObjectFromCDF(x, type = objecttype, band = band, ...) : > cells are not equally spaced; you should extract values as points > > # I check some solutions that recomend force read the file with the > argument "stopIfNotEqualSpaced = F" as bellow. > sf <- stack("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc", > *stopIfNotEqualSpaced > = F*) > bf <- brick("pr_Amon_GFDL-ESM2G_rcp45_r1i1p1_204101-204512.nc", > *stopIfNotEqualSpaced > = F*) > > I performed some tests and only "works" with brick. However I did not find > any solution to check where is the problem and fix it in the file. > > How can I check if it is really an issue and fix it? > > sessionInfo() > R version 3.5.1 (2018-07-02) > Platform: x86_64-w64-mingw32/x64 (64-bit) > Running under: Windows 10 x64 (build 18362) > > Matrix products: default > > locale: > [1] LC_COLLATE=Portuguese_Brazil.1252 LC_CTYPE=Portuguese_Brazil.1252 > [3] LC_MONETARY=Portuguese_Brazil.1252 LC_NUMERIC=C > [5] LC_TIME=Portuguese_Brazil.1252 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] raster_2.8-19 sp_1.3-1 > > loaded via a namespace (and not attached): > [1] compiler_3.5.1 rgdal_1.4-3 Rcpp_1.0.1 codetools_0.2-15 > ncdf4_1.16.1 > [6] grid_3.5.1 lattice_0.20-35 > > Best regards, > > Frederico > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] regression-kriging and co-kriging (Edzer Pebesma)
On 8/15/19 12:31 PM, Emanuele Barca wrote: > Dear Edzer, > > sorry for bothering you once more but I need to be sure about my R script. > > In summary, i'm comparing the performance of regression-kriging and > collocated co-kriging. > > Regression-kriging was based on an unique covariate, the elevation Z. > > I use Z as unique ancillary variable in the Co-kriging. > > As first attempt, the final raster maps were completely different. It > appeared that it was due > > to the fact that the dataset was non-stationary and only > regression-kriging overcomes this issue, while co-kriging not. > > But if I pass to universal co-kriging introducing Z as covariate, it > bacomes useless as ancillary variable! > > What is my mistake? You assume that someone else can be sure about your analysis by looking at your R script without having access to your data. And you are starting a personal dialogue on a list, essentially uninviting everyone else to get involved. > > emanuele > > > > > Il 2019-08-15 12:00 r-sig-geo-requ...@r-project.org ha scritto: >> Send R-sig-Geo mailing list submissions to >> r-sig-geo@r-project.org >> >> To subscribe or unsubscribe via the World Wide Web, visit >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> or, via email, send a message with subject or body 'help' to >> r-sig-geo-requ...@r-project.org >> >> You can reach the person managing the list at >> r-sig-geo-ow...@r-project.org >> >> When replying, please edit your Subject line so it is more specific >> than "Re: Contents of R-sig-Geo digest..." >> >> >> Today's Topics: >> >> 1. Error running codes (Enoch Gyamfi Ampadu) >> 2. Re: regression-kriging and co-kriging (Edzer Pebesma) >> >> -- >> >> Message: 1 >> Date: Thu, 15 Aug 2019 06:51:49 +0200 >> From: Enoch Gyamfi Ampadu >> To: r-sig-geo@r-project.org >> Subject: [R-sig-Geo] Error running codes >> Message-ID: >> >> Content-Type: text/plain; charset="utf-8" >> >> Dear List, >> >> Please I have been running a certain the codes below to read my image, >> shapeffile to enable me classify for cover with Random forest. I have >> gotten to the point of extracting the raster values and the raster >> information. However I keep getting errors. though I have tried different >> combination and other codes like readOGR for importing the training >> polygon. I will be glad if you could be of help as I am new to using R. >> >> #import clip image of study area >> >> TestImg3 <- brick("C:\\Users\\hp\\Desktop\\Data collection\\Nkandla >> Images\\Landsat\\2019\\08052019\\Corrected\\New >> Folder\\Image08052019_Clip.tif") >> >> #Assign name of bands >> names(TestImg3) <- c(paste0("B", 1:7, coll=""), "B9") >> plotRGB(TestImg3, r=4, g=3, b=2, stretch ="lin") >> >> #Import training shapefile >> sample <- shapefile("C:/Users/hp/Desktop/Data collection/Nkandla >> Images/Landsat/2019/08052019/Corrected/Train08052019_2.shp") >> >> responseCol <- "class" >> >> #(I have tried options of changing "class" to "classname" as reflect >> actaul >> name assigned in ArcMap) >> >> # Overlap the sample polygons on the image (combine the class information >> with extracted values). >> >> pixels = data.frame(matrix(vector(), nrow = 0, ncol = >> length(names(img)) + >> 1)) >> for (i in 1:length(unique(sample[[responseCol]]))){ >> category <- unique(sample[[responseCol]])[i] >> categorymap <- sample[sample[[responseCol]] == category,] >> dataSet <- extract(img, categorymap) >> dataSet <- dataSet[!unlist(lapply(dataSet, is.null))] >> if(is(sample, "SpatialPointsDataFrame")){ >> dataSet <- cbind(dataSet, class = as.numeric(category)) >> pixeles <- rbind(pixeles, dataSet) >> } >> if(is(sample, "SpatialPolygonsDataFrame")){ >> dataSet <- lapply(dataSet, function(x){cbind(x, class = >> as.numeric(rep(category, >> >> nrow(x}) >> df <- do.call("rbind", dataSet) >> pixels <- rbind(pixeles, df) >> } >> >> } >> >> THIS IS THE ERROR I GET FROM RUNNING THE ABOVE CODES >> >>> for (i in 1:length(unique(samples[[responseCol]]))){ >> + category <- unique(samples[[responseCol]])[
Re: [R-sig-Geo] regression-kriging and co-kriging
On 8/12/19 8:21 PM, Emanuele Barca wrote: > Dear Edzer, > > maybe I found the solution. I found this in the predict function help: > "When a non-stationary (i.e., non-constant) mean is used, both for > simulation and prediction purposes the variogram model defined should be > that of the residual process, and not that of the raw observations" > Since my data were, actually, non-stationary, I applied the universal > co-kriging instead usual co-kriging. > now the maps of regression-kring and co-kriging are actually similar s > expected. > did I understand correctly the quoted sentence? I think so, but hard to be sure given the information you provide. > > regards > > emanuele barca > -- >> >> Message: 2 >> Date: Sat, 10 Aug 2019 10:41:38 +0200 >> From: Edzer Pebesma >> To: r-sig-geo@r-project.org >> Subject: Re: [R-sig-Geo] regression-kriging and co-kriging >> Message-ID: >> Content-Type: text/plain; charset="utf-8" >> >> Hard to tell from your script. Maybe give a reproducible example? >> >> On 8/6/19 1:07 PM, Emanuele Barca wrote: >>> Dear r-sig-geo friends, >>> >>> I produced two maps garnered in the following way: >>> >>> # for regression-kriging >>> Piezo.map <-autoKrige(LivStat ~ Z, input_data = mydata.sp, new_data >>> = covariates, model = "Ste") >>> >>> Piezork.pred <- Piezo.map$krige_output$var1.pred >>> Piezork.coords <- Piezo.map$krige_output@coords >>> Piezork.out <- as.data.frame(cbind(Piezork.coords, Piezork.pred)) >>> colnames(Piezork.out)[1:2] <- c("X", "Y") >>> coordinates(Piezork.out) = ~ X + Y >>> gridded(Piezork.out) <- TRUE >>> >>> spplot(Piezork.out, main = list(label = "R-k Hydraulic head", cex = >>> 1.5)) >>> >>> #for co-kriging >>> g <- gstat(id = "Piezo", formula = LivStat ~ 1, data = mydata.sp, set >>> = list(nocheck = 1)) >>> g <- gstat(g, id = "Z", formula = Z ~ 1, data = mydata.sp, set = >>> list(nocheck = 1)) >>> >>> v.g <- variogram(g) >>> >>> #g <- gstat(g, id = "Piezo", model = vgm(150, "Mat", 1350, 0.0, kappa >>> = 1.9), fill.all = T)# >>> g <- gstat(g, id = "Piezo", model = vgm(0.7, "Ste", 1300, 18, kappa = >>> 1.9), fill.all = T)# >>> g.fit <- fit.lmc(v.g, g, fit.lmc = TRUE, correct.diagonal = 1.01) # >>> fit multivariable variogram model , fit.lmc = TRUE, correct.diagonal >>> = 1.01 >>> g.fit >>> plot(v.g, model = g.fit, main = "Fitted Variogram Models - Raw Data")# >>> #gridded(covariates) <- TRUE >>> g.cok <- predict(g.fit, newdata = covariates)#grid >>> >>> g.cok.pred <- g.cok@data$Piezo.pred >>> <- na.omit(g.cok.pred) >>> g.cok.coords <- g.cok@coords >>> g.cok.out <- as.data.frame(cbind(g.cok.coords, g.cok.pred)) >>> colnames(g.cok.out)[1:2] <- c("X", "Y") >>> coordinates(g.cok.out) = ~ X + Y >>> gridded(g.cok.out) <- TRUE >>> spplot(g.cok.out, main = list(label = "Hydraulic head with >>> Co-kriging", cex = 1.5)) >>> >>> ### >>> >>> >>> I am unable to understand why the first map appears as a raster and >>> the second not, notwithstanding the fact that they are both computed >>> on the same "covariates" DEM??? >>> >>> where is the mistake??? >>> >>> regards >>> >>> emanuele >>> >>> >>> Emanuele Barca Researcher >>> Water Research Institute (IRSA-CNR) >>> Via De Blasio, 5 70123 Bari (Italy) >>> Phone +39 080 5820535 Fax +39 080 5313365 >>> Mobile +39 340 3420689 >>> _ >>> >>> >>> >>> --- >>> Questa e-mail è stata controllata per individuare virus con Avast >>> antivirus. >>> https://www.avast.com/antivirus >>> >>> [[alternative HTML version deleted]] >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo
Re: [R-sig-Geo] regression-kriging and co-kriging
Hard to tell from your script. Maybe give a reproducible example? On 8/6/19 1:07 PM, Emanuele Barca wrote: > Dear r-sig-geo friends, > > I produced two maps garnered in the following way: > > # for regression-kriging > Piezo.map <-autoKrige(LivStat ~ Z, input_data = mydata.sp, new_data > = covariates, model = "Ste") > > Piezork.pred <- Piezo.map$krige_output$var1.pred > Piezork.coords <- Piezo.map$krige_output@coords > Piezork.out <- as.data.frame(cbind(Piezork.coords, Piezork.pred)) > colnames(Piezork.out)[1:2] <- c("X", "Y") > coordinates(Piezork.out) = ~ X + Y > gridded(Piezork.out) <- TRUE > > spplot(Piezork.out, main = list(label = "R-k Hydraulic head", cex = 1.5)) > > #for co-kriging > g <- gstat(id = "Piezo", formula = LivStat ~ 1, data = mydata.sp, set > = list(nocheck = 1)) > g <- gstat(g, id = "Z", formula = Z ~ 1, data = mydata.sp, set = > list(nocheck = 1)) > > v.g <- variogram(g) > > #g <- gstat(g, id = "Piezo", model = vgm(150, "Mat", 1350, 0.0, kappa > = 1.9), fill.all = T)# > g <- gstat(g, id = "Piezo", model = vgm(0.7, "Ste", 1300, 18, kappa = > 1.9), fill.all = T)# > g.fit <- fit.lmc(v.g, g, fit.lmc = TRUE, correct.diagonal = 1.01) # > fit multivariable variogram model , fit.lmc = TRUE, correct.diagonal = 1.01 > g.fit > plot(v.g, model = g.fit, main = "Fitted Variogram Models - Raw Data")# > #gridded(covariates) <- TRUE > g.cok <- predict(g.fit, newdata = covariates)#grid > > g.cok.pred <- g.cok@data$Piezo.pred > <- na.omit(g.cok.pred) > g.cok.coords <- g.cok@coords > g.cok.out <- as.data.frame(cbind(g.cok.coords, g.cok.pred)) > colnames(g.cok.out)[1:2] <- c("X", "Y") > coordinates(g.cok.out) = ~ X + Y > gridded(g.cok.out) <- TRUE > spplot(g.cok.out, main = list(label = "Hydraulic head with > Co-kriging", cex = 1.5)) > > ### > > I am unable to understand why the first map appears as a raster and > the second not, notwithstanding the fact that they are both computed > on the same "covariates" DEM??? > > where is the mistake??? > > regards > > emanuele > > > Emanuele Barca Researcher > Water Research Institute (IRSA-CNR) > Via De Blasio, 5 70123 Bari (Italy) > Phone +39 080 5820535 Fax +39 080 5313365 > Mobile +39 340 3420689 > _ > > > > --- > Questa e-mail è stata controllata per individuare virus con Avast antivirus. > https://www.avast.com/antivirus > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] st_sample with raster
sf::st_sample has a size argument, which can be specified as a vector to operate per feature. If you give that a value equal to st_area(obj) * obj$required_density and maybe set exact = FALSE, you should get the densities you want. On 7/30/19 12:25 AM, Tyler Frazier via R-sig-Geo wrote: > Hi All, > > Spatstat::rpoint has an argument f that accepts a pixel image object (class > .im) to generate points based on the proportionate probability density of > that raster. Just wondering, does sf::st_sample or perhaps another function > have a similar available method? > > Thanks so much! > Ty > > Tyler Frazier, Ph.D. > Lecturer of Interdisciplinary Studies > Data Science Program > College of William and Mary > Webpage: http://tjfrazier.people.wm.edu > Phone: +01 (757) 386-1269 > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to reshape a LINESTRING sfc into equal length segments
You may want to look into sf::st_segmentize. On 7/3/19 10:58 AM, Andrea Gilardi wrote: > Hi everyone. This is my first mail to this mailing list so excuse me if I'm > doing something wrong. I was wondering if it's possible to reshape a > LINESTRING sfc of all connected segments into equal length segments. > > I tried to explain a little bit on why I want to do that here: > https://gis.stackexchange.com/questions/327376/how-to-divide-linestring-spatial-network-into-equal-length-segments-using-r-sf > > I'm pretty sure that I'm overcomplicating something simple and, in that > case, could you point me in the right direction? > > Thank you in advance > Andrea > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] testing a crs string
On 6/20/19 4:18 PM, Roy Mendelssohn - NOAA Federal via R-sig-Geo wrote: > Hi All: > > I am extending a package I have to allow users to provide a crs string, which > will eventually through some other packages be sent to sf. I like to try and > make at least minimal tests that inputs are valid. Is there some existing > code that will check if a string is a valid crs string? > inherits(try(sf::st_crs("+proj=invalid"), silent = TRUE), "try-error") [1] TRUE > inherits(try(sf::st_crs("+proj=longlat"), silent = TRUE), "try-error") [1] FALSE > > Thanks, > > -Roy > > > ** > "The contents of this message do not reflect any position of the U.S. > Government or NOAA." > ** > Roy Mendelssohn > Supervisory Operations Research Analyst > NOAA/NMFS > Environmental Research Division > Southwest Fisheries Science Center > ***Note new street address*** > 110 McAllister Way > Santa Cruz, CA 95060 > Phone: (831)-420-3666 > Fax: (831) 420-3980 > e-mail: roy.mendelss...@noaa.gov www: http://www.pfeg.noaa.gov/ > > "Old age and treachery will overcome youth and skill." > "From those who have been given much, much will be expected" > "the arc of the moral universe is long, but it bends toward justice" -MLK Jr. > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] install rgdal from source macOS
D_PROJ_API_H > -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sp/include" > -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk > -I/usr/local/include -fPIC -Wall -g -O2 -c ogr_proj.cpp -o ogr_proj.o > clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" > -DNDEBUG -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H > -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sp/include" > -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk > -I/usr/local/include -fPIC -Wall -g -O2 -c ogrdrivers.cpp -o ogrdrivers.o > clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" > -DNDEBUG -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H > -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sp/include" > -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk > -I/usr/local/include -fPIC -Wall -g -O2 -c ogrsource.cpp -o ogrsource.o > clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" > -DNDEBUG -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H > -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sp/include" > -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk > -I/usr/local/include -fPIC -Wall -g -O2 -c proj_info6.cpp -o proj_info6.o > clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" > -DNDEBUG -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -I/Users/dosc3612/Applications/miniconda3/envs/rgdal/include > -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H > -I"/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sp/include" > -isysroot /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk > -I/usr/local/include -fPIC -Wall -g -O2 -c projectit.cpp -o projectit.o > clang++ -std=gnu++11 -dynamiclib -Wl,-headerpad_max_install_names > -undefined dynamic_lookup -single_module -multiply_defined suppress > -L/Library/Frameworks/R.framework/Resources/lib -L/usr/local/lib -o > rgdal.so OGR_write.o gdal-bindings.o init.o inverser.o local_stubs.o > ogr_geom.o ogr_polygons.o ogr_proj.o ogrdrivers.o ogrsource.o proj_info6.o > projectit.o -L/Users/dosc3612/Applications/miniconda3/envs/rgdal/lib -lgdal > -L/Users/dosc3612/Applications/miniconda3/envs/rgdal/lib -lproj > -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework > -Wl,CoreFoundation > installing to > /Library/Frameworks/R.framework/Versions/3.6/Resources/library/00LOCK-rgdal/00new/rgdal/libs > ** R > ** data > ** inst > ** byte-compile and prepare package for lazy loading > ** help > *** installing help indices > ** building package indices > ** installing vignettes > ** testing if installed package can be loaded from temporary location > Error: package or namespace load failed for ‘rgdal’ in dyn.load(file, > DLLpath = DLLpath, ...): > unable to load shared object > '/Library/Frameworks/R.framework/Versions/3.6/Resources/library/00LOCK-rgdal/00new/rgdal/libs/rgdal.so': > > dlopen(/Library/Frameworks/R.framework/Versions/3.6/Resources/library/00LOCK-rgdal/00new/rgdal/libs/rgdal.so, > 6): Library not loaded: @rpath/libgdal.20.dylib > Referenced from: > /Library/Frameworks/R.framework/Versions/3.6/Resources/library/00LOCK-rgdal/00new/rgdal/libs/rgdal.so > Reason: image not found > Error: loading failed > Execution halted > ERROR: loading failed > * removing > ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal’ > > The downloaded source packages are in > ‘/private/var/folders/bb/13z2kq0j01jg5q_0ypjc76grgn/T/RtmplsyMiu/downloaded_packages’ > Warning message: > In install.packages("rgdal", type = "source", configure.args = > c("--with-proj-include=/Users/dosc3612/Applications/miniconda3/envs/rgdal/include", > : > installation of package ‘rgdal’ had non-zero exit status > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spacetime and sf
Dear ƒacu, spacetime will remain the way it is, just like sp, too much stuff depends on how it is now. For STFDF objects in spacetime, a follow-up package is stars: it supports raster and vector data cubes, and can deal with more cube dimensions than STFDF ever could; it also gets better dealing with e.g. large satellite images. For the STIDF (irregular: s/t point patterns) there is currently no follow-up package that builds on top of sf. As you note, its functionality would be pretty much limited to registering which column denotes time in an attribute, and add a class label. So far, I haven't considered that worth the effort: you may have to add lots of method instances to have that information survive all the current methods available for sf objects. Have you thought about giving the time variable a fixed name? Or is there a need for maintaining multiple time variables in sf objects? Best, On 6/4/19 5:01 PM, Facundo Muñoz wrote: > Dear Edzer, > > How do you conceive the future of the package `spacetime` and more > generally, the representation of spatio-temporal data, in light of the > transition from `sp` to `sf`? > > For context, I'm working on a package that handles spatio-temporal point > patterns. For my analyses I have relied on `sf` with some additional > attribute for the dates. However, for packaging I think it is safer to > rely on explicit classes for spatio-temporal data (so I don't have to > worry all the time about which is the temporal variable, in which > format, etc.). But it seems that only `spacetime` provides such classes > generically, which in turn relies on `sp` for the spatial dimension. Is > this going to continue to be supported or are there other plans? > > Thanks in advance > > ƒacu.- > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] GDAL 3.0.0 and rgdal?
The development version installs with GDAL 3, try install.packages("rgdal", repos="http://R-Forge.R-project.org;) On 5/27/19 1:43 AM, Thiago V. dos Santos via R-sig-Geo wrote: > Dear all, > Yesterday I replaced GDAL v2.4.1 with the latest v3.0.0 on my Linux server. > Then, when I tried to reinstall rgdal, I got an error due to wrong system > requirements: rgdal requires GDAL between versions 1.11.4 and 2.5.0. > Therefore, I am wondering: are there plans for rgdal to support GDAL 3.0.0 in > the next few days? I like to be up-to-date with GIS software, but I don't > really need to use GDAL 3.0.0 at all. If it is going to be a little while > until it is supported by rgdal, I can happily revert back to v2.4.1. > Thanks, -- Thiago V. dos Santos > Postdoctoral Research FellowDepartment of Climate and Space Science and > EngineeringUniversity of Michigan > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Mark variograms: 3D(xyt) points + 1D marks
On 5/11/19 6:11 PM, Tim Pollington wrote: > Thank you Edzer for your suggestion however I'm trying to avoid > geostatistical methods because my data is disease cases recorded at > household locations. Therefore strictly speaking there is not a random > field that pervades the space between the households and so I must treat > it as a spatiotemporal marked point process instead. I attemping to > extend Stoyan et al's mark variograms for ST > processes https://www.sciencedirect.com/science/article/pii/S2211675317300696 > to > my data. Regrettably (an ironically, as an associate editor) I have no longer access to that journal. Did you contact the authors of the paper whether they have software they're willing to share? > > Hopefully I'll have gotten a bit further by July so perhaps I can show > you more if you're attending Spatial Statistics in Sitges. Looking forward to that! > > Kind regards, > > Tim. > -------- > *From:* Edzer Pebesma > *Sent:* 10 May 2019 13:18 > *To:* r-sig-geo@r-project.org > *Subject:* Re: [R-sig-Geo] Mark variograms: 3D(xyt) points + 1D marks > > No, but you can model them as a function of space and time, see e.g. > https://journal.r-project.org/archive/2016/RJ-2016-014/index.html > > On 5/10/19 12:53 PM, Tim Pollington wrote: >> On reflection please ignore my question. One can't construct an xyt >> variogram with additional marks as one can't define a proper distance in an >> xyt space as different units are involved. >> >> [[alternative HTML version deleted]] >> >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > -- > Edzer Pebesma > Institute for Geoinformatics > Heisenbergstrasse 2, 48151 Muenster, Germany > Phone: +49 251 8333081 -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Semi-join?
On 5/7/19 8:29 PM, Tim Keitt wrote: > I keep finding things that are trivial in sql are not so in R, so perhaps > you can help me out. > > What would be the way in the 'sf' package to retain only polygons from one > coverage that intersect points in another coverage? I think that's a > semi-join. > > Roughly: 'select a.* from polys a, points b where st_intersects(a.geom, > b.geom)' a[b,] # or st_join(a, b, left = FALSE) # no left join gives you inner join > > Thanks. > > THK > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial join/intersection in R
On 5/6/19 12:23 PM, Edzer Pebesma wrote: > # raster does the same as stars: That was too quick: raster::extract returns cell values, stars::st_intersects returns cell indexes. If we fill the raster with values 100:1 rather than 1:100, we see library(stars) # Loading required package: abind # Loading required package: sf # Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0 # create 10 x 10 raster st_as_stars(matrix(100:1, 10, 10)) %>% st_set_dimensions(names = c("x", "y")) %>% st_set_dimensions("x", values = 0:9) %>% st_set_dimensions("y", values = 9:0) -> r # three points: at node, edge, and inside cell: p3 = st_sfc(st_point(c(1,1)), st_point(c(1, 1.5)), st_point(c(1.1, 1.8))) st_intersects(r, p3, transpose = TRUE) # Sparse geometry binary predicate list of length 3, where the predicate was `intersects' # 1: 82 # 2: 72 # 3: 72 library(raster) # Loading required package: sp extract(as(r, "Raster"), as(p3, "Spatial")) # [1] 19 29 29 -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Spatial join/intersection in R
y way > to avoid this? > > Thanks in advance. > Roozbeh > > > image.png > > > > -- > *Roozbeh Valavi* > PhD Candidate > The Quantitative & Applied Ecology Group <http://qaeco.com/> > School of BioSciences | Faculty of Science > The University of Melbourne, VIC 3010, Australia > Mobile: +61 423 283 238 > ___________ > 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 > > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 Rplots.pdf Description: Adobe PDF document pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] R-sig-Geo Digest, Vol 189, Issue 3
i, > PhD Forest Science - Ecological Mathematics > Skype ID: maurizioxyz > linux user 552742 > #EUFGIS national Focal point for Italy > #Annals of Silvicultural Research Associated Editor > https://scholar.google.it/citations?hl=en=_2X6fu8J > <https://scholar.google.it/citations?hl=en=_2X6fu8J*> > > > > > Il giorno sab 4 mag 2019 alle ore 12:10 <mailto:r-sig-geo-requ...@r-project.org>> ha scritto: > > Send R-sig-Geo mailing list submissions to > r-sig-geo@r-project.org <mailto:r-sig-geo@r-project.org> > > To subscribe or unsubscribe via the World Wide Web, visit > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > or, via email, send a message with subject or body 'help' to > r-sig-geo-requ...@r-project.org > <mailto:r-sig-geo-requ...@r-project.org> > > You can reach the person managing the list at > r-sig-geo-ow...@r-project.org > <mailto:r-sig-geo-ow...@r-project.org> > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of R-sig-Geo digest..." > > > Today's Topics: > > 1. Re: Strange spatial reference system netCDF (Edzer Pebesma) > > -- > > Message: 1 > Date: Fri, 3 May 2019 16:58:01 +0200 > From: Edzer Pebesma <mailto:edzer.pebe...@uni-muenster.de>> > To: r-sig-geo@r-project.org <mailto:r-sig-geo@r-project.org> > Subject: Re: [R-sig-Geo] Strange spatial reference system netCDF > Message-ID: <82193e8b-8fcc-3308-7874-a22ae6985...@uni-muenster.de > <mailto:82193e8b-8fcc-3308-7874-a22ae6985...@uni-muenster.de>> > Content-Type: text/plain; charset="utf-8" > > Dear Maurizio, > > an attempt to do something along this line is: > > f = "tas_rcp85_land-rcm_eur_12km_01_mon_198012-208011.nc > <http://tas_rcp85_land-rcm_eur_12km_01_mon_198012-208011.nc>" > library(stars) > r = read_ncdf(f, ncsub = cbind(start = c(1, 1, 1, 1), count = c(418, > 406, 3, 1)), eps=1e-3) > rx = read_stars(f, proxy = TRUE) # only for the crs! > st_crs(r) = st_crs(rx) > r0 = stars:::st_transform_proj.stars(r, 4326) > png("x.png") > plot(r0[,,,1], border = NA, axes = TRUE, reset = FALSE) > library(rnaturalearth) > plot(ne_coastline(returnclass = "sf"), add = TRUE, col = 'orange') > > Note that the output object time stamps respect the 360-day calendar > (using pkg PCICt). > > You'll find some discussion, and the output image, here: > https://github.com/r-spatial/stars/issues/175 > > Feel free to post follow-up questions here, or to github. > > On 5/2/19 6:49 PM, Maurizio Marchi wrote: > > Dear list, > > I'm working with large netCDF files from the UK Climate > Projections portal ( > > https://www.metoffice.gov.uk/research/collaboration/ukcp). More in > detail > > I'm reading with the ncdf4 library this > > > > <http://data.ceda.ac.uk/badc/ukcp18/data/land-rcm/eur/12km/rcp85/01/tas/mon/latest/> > > file. Once loaded O can easily handle it but it is a strange reference > > system as follow: > > rotated_latitude_longitude[] > > grid_mapping_name: rotated_latitude_longitude > > longitude_of_prime_meridian: 0 > > earth_radius: *6371229* > > grid_north_pole_latitude: 39.25 > > grid_north_pole_longitude: 198 > > north_pole_grid_longitude: 0 > > If I load the .nc file in QGIS I see that the SR is "+proj=longlat > > +a=6371229 +b=6371229 +no_defs" where the values associated to *a* and > > *b* variables > > are the earth radius. Even if QGIS can read it, the raw file isn't > > projected properly (it seems tha QGIS is not able to handle such > > projection). > > Finally, using the ncdf4+raster libraries, I can easily generate e > raster > > image but then I'm not able to understand I could I re project > this raster > > in WGS84 reference system. > > Any helps? > > Cheers > > > > -- > > Maurizio Marchi, > > PhD Forest Science - Ecological Mathematics > > Skype ID: maurizioxyz > > linux user 552742 > > #EUFGIS national Focal point for Italy > > #Annals of Silvicultural Research Associated Editor > > https
Re: [R-sig-Geo] Strange spatial reference system netCDF
Dear Maurizio, an attempt to do something along this line is: f = "tas_rcp85_land-rcm_eur_12km_01_mon_198012-208011.nc" library(stars) r = read_ncdf(f, ncsub = cbind(start = c(1, 1, 1, 1), count = c(418, 406, 3, 1)), eps=1e-3) rx = read_stars(f, proxy = TRUE) # only for the crs! st_crs(r) = st_crs(rx) r0 = stars:::st_transform_proj.stars(r, 4326) png("x.png") plot(r0[,,,1], border = NA, axes = TRUE, reset = FALSE) library(rnaturalearth) plot(ne_coastline(returnclass = "sf"), add = TRUE, col = 'orange') Note that the output object time stamps respect the 360-day calendar (using pkg PCICt). You'll find some discussion, and the output image, here: https://github.com/r-spatial/stars/issues/175 Feel free to post follow-up questions here, or to github. On 5/2/19 6:49 PM, Maurizio Marchi wrote: > Dear list, > I'm working with large netCDF files from the UK Climate Projections portal ( > https://www.metoffice.gov.uk/research/collaboration/ukcp). More in detail > I'm reading with the ncdf4 library this > <http://data.ceda.ac.uk/badc/ukcp18/data/land-rcm/eur/12km/rcp85/01/tas/mon/latest/> > file. Once loaded O can easily handle it but it is a strange reference > system as follow: > rotated_latitude_longitude[] > grid_mapping_name: rotated_latitude_longitude > longitude_of_prime_meridian: 0 > earth_radius: *6371229* > grid_north_pole_latitude: 39.25 > grid_north_pole_longitude: 198 > north_pole_grid_longitude: 0 > If I load the .nc file in QGIS I see that the SR is "+proj=longlat > +a=6371229 +b=6371229 +no_defs" where the values associated to *a* and > *b* variables > are the earth radius. Even if QGIS can read it, the raw file isn't > projected properly (it seems tha QGIS is not able to handle such > projection). > Finally, using the ncdf4+raster libraries, I can easily generate e raster > image but then I'm not able to understand I could I re project this raster > in WGS84 reference system. > Any helps? > Cheers > > -- > Maurizio Marchi, > PhD Forest Science - Ecological Mathematics > Skype ID: maurizioxyz > linux user 552742 > #EUFGIS national Focal point for Italy > #Annals of Silvicultural Research Associated Editor > https://scholar.google.it/citations?hl=en=_2X6fu8J > <https://scholar.google.it/citations?hl=en=_2X6fu8J*> > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] raster - units assignment error in reading NetCDF4
It looks like this is a curvilinear grid; you could read & plot it with library(stars) r = read_stars("L8_OLI_2017_03_18_00_09_19_093086_L2R.nc", curvilinear = c("lon", "lat")) plot(r[1], border = NA) Printing r shows the units of each of the variables. The plotting takes a while because the data grid is not aligned with the pixel grid of the plotting device. On 4/30/19 12:43 PM, Alexandre Castagna wrote: > Hi group, > > I'm trying to read a NetCDF4 file, but I get the following error message: > > fl <- 'L8_OLI_2017_03_18_00_09_19_093086_L2R.nc' > r <- raster(fl, varname = 'rhos_443') > > Error in (function (cl, name, valueClass) : > assignment of an object of class “numeric” is not valid for @‘unit’ in an > object of class “.SingleLayerData”; is(value, "character") is not TRUE > > This data is OLI Landsat 8 atmospheric 'corrected' imagery with the ACOLITE > software (link <https://github.com/acolite/acolite>). > In a previous version of the software, exploring the file after nc_open > function revealed that units was numeric ('num' qualifier was shown in R > object print); but since, the developer has committed a change for all > units to be written as characters. The unit in question is '1' (unitless, > for bihemispherical reflectance). A new exploration after the update to > check the units type does not show the qualifier 'num' anymore, but the > error keep the same. Maybe R is translating '1' to 1 automatically? > > I can go around and create a raster directly from data read with ncdf4 > package, like: > > ncfl <- nc_open(fl) > r <- raster(t(ncvar_get(ncfl, 'rhos_443'))) > > But for a number of reasons that is less ideal. A temporary link to > download the example data file is: https://we.tl/t-nENAV7tmBg > > Kind regards, > > Alexandre Castagna Mourão e Lima > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] tpk files
Thanks; not sure this is what is intended, but it seemed to work: > library(sf) Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0 > r = st_read("EMU.gpkg", query = "select * from EMU_Master where depth_lvl = 1") Reading layer `EMU_Master' from data source `/home/edzer/Downloads/EMU.gpkg' using driver `GPKG' Simple feature collection with 677109 features and 33 fields geometry type: POINT dimension: XYZ bbox: xmin: -179.875 ymin: -78.375 xmax: 179.875 ymax: 89.875 epsg (SRID):4326 proj4string:+proj=longlat +datum=WGS84 +no_defs > object.size(r) 438817848 bytes On 4/12/19 12:48 AM, Barry Rowlingson wrote: > QGIS could "read" it pretty smartly, only loading in the bits in the > current view extent, and doing the loading in a separate thread so the > GUI was still active. > > Marta has told me that she only needs the surface layer points - the X Y > and Z coordinates are duplicated as attributes so I think a selection of > Z=0 (or something) might make a small enough subset to read into R > directly. This might be doable via the `query` option of `st_read`, or > failing that a "select" in SQLite3 and dumping the results to a CSV > file. There seems to be hundreds of data points at each location so > taking the surface points only could thin it out considerably. It may > still take some time but if its a one-off... > > Barry > > > On Thu, Apr 11, 2019 at 9:57 PM Edzer Pebesma > mailto:edzer.pebe...@uni-muenster.de>> > wrote: > > It's a 30 Gb 3D point file: > > $ ogrinfo EMU.gpkg > INFO: Open of `EMU.gpkg' > using driver `GPKG' successful. > 1: EMU_Master (3D Point) > > @Shaun: homebrew seems to be supported neither by apple, nor by CRAN, so > you are a bit on your own there. Have you tried the CRAN binary packages > using GDAL? > > In any case, the windows binary does (should) support gpkg, see > https://github.com/rwinlib/gdal2 > > Trying to read this file into R with sf::st_read will require a lot of > RAM, or some strategy to read in parts only. > > On 4/11/19 5:52 PM, Shaun Walbridge wrote: > > I think the issue is, most GDAL installations don't have the > Geopackage raster driver [1] installed by default, which lists > "Needs libsqlite3 (and any or all of PNG, JPEG, WEBP drivers)" for > it to be available. At least on my Homebrew installation of GDAL, > this driver wasn't built out of the box. If you rebuild GDAL with > this additional driver, or find a prebuilt binary which has it, it > should be able to open. A simple test is if `gdalinfo EMU.gpkg` > returns information about the dataset outside of R. > > > > 1. > > https://urldefense.proofpoint.com/v2/url?u=https-3A__www.gdal.org_drv-5Fgeopackage-5Fraster.html=DwIGaQ=n6-cguzQvX_tUIrZOS_4Og=fCPRb7QX-vd5bnO9gIJHCiX852SVUtyYX--xtCKtpfk=p5ULiF5de1gKZBP-IzWbMO9Pe5LFzv9uaZ5VJYnWw1Y=d6xaKGlN0jpd8mBdjKXAhzst7N3Bgo43BvJlLnDSngk= > > > > On 4/11/19, 11:41 AM, "Barry Rowlingson" <mailto:b.rowling...@gmail.com>> wrote: > > > > What did you try? The instructions at the top say: > > > > "Download 3.3GB tile package and rename extension from .tpk to > .zip. > > Extract to get EMU.gpkg" > > > > If that's a valid GeoPackage then `sf` should be able to read > it. Not sure > > what might be in the geopackage though, "tile package" sounds > like rasters, > > but GeoPackages are generally vector... > > > > I'll try in five minutes when all 3.3Gb have downloaded > > > > On Thu, Apr 11, 2019 at 3:37 PM Marta Rufino > mailto:marta.m.ruf...@gmail.com>> > > wrote: > > > > > Hi, > > > > > > I would like to open (and use) a 'tpk' file from arcgis in r. > > > For example: > > > > > > > > https://esri.maps.arcgis.com/home/item.html?id=24885cd6bd9544f5a8e15d0bf40f67d6 > > > > > > I tried raster and sf package, but no luck. > > > > > > Any ideia if we can do this in r? > > > > > > Thank you very much in advance, > > > > > > Best wishes, > > > M. > > > > > > [[alternative HTML version deleted]] > > > > > > ___ > > > R-sig-G
Re: [R-sig-Geo] tpk files
It's a 30 Gb 3D point file: $ ogrinfo EMU.gpkg INFO: Open of `EMU.gpkg' using driver `GPKG' successful. 1: EMU_Master (3D Point) @Shaun: homebrew seems to be supported neither by apple, nor by CRAN, so you are a bit on your own there. Have you tried the CRAN binary packages using GDAL? In any case, the windows binary does (should) support gpkg, see https://github.com/rwinlib/gdal2 Trying to read this file into R with sf::st_read will require a lot of RAM, or some strategy to read in parts only. On 4/11/19 5:52 PM, Shaun Walbridge wrote: > I think the issue is, most GDAL installations don't have the Geopackage > raster driver [1] installed by default, which lists "Needs libsqlite3 (and > any or all of PNG, JPEG, WEBP drivers)" for it to be available. At least on > my Homebrew installation of GDAL, this driver wasn't built out of the box. If > you rebuild GDAL with this additional driver, or find a prebuilt binary which > has it, it should be able to open. A simple test is if `gdalinfo EMU.gpkg` > returns information about the dataset outside of R. > > 1. > https://urldefense.proofpoint.com/v2/url?u=https-3A__www.gdal.org_drv-5Fgeopackage-5Fraster.html=DwIGaQ=n6-cguzQvX_tUIrZOS_4Og=fCPRb7QX-vd5bnO9gIJHCiX852SVUtyYX--xtCKtpfk=p5ULiF5de1gKZBP-IzWbMO9Pe5LFzv9uaZ5VJYnWw1Y=d6xaKGlN0jpd8mBdjKXAhzst7N3Bgo43BvJlLnDSngk= > > On 4/11/19, 11:41 AM, "Barry Rowlingson" wrote: > > What did you try? The instructions at the top say: > > "Download 3.3GB tile package and rename extension from .tpk to .zip. > Extract to get EMU.gpkg" > > If that's a valid GeoPackage then `sf` should be able to read it. Not sure > what might be in the geopackage though, "tile package" sounds like > rasters, > but GeoPackages are generally vector... > > I'll try in five minutes when all 3.3Gb have downloaded > > On Thu, Apr 11, 2019 at 3:37 PM Marta Rufino > wrote: > > > Hi, > > > > I would like to open (and use) a 'tpk' file from arcgis in r. > > For example: > > > > > https://esri.maps.arcgis.com/home/item.html?id=24885cd6bd9544f5a8e15d0bf40f67d6 > > > > I tried raster and sf package, but no luck. > > > > Any ideia if we can do this in r? > > > > Thank you very much in advance, > > > > Best wishes, > > M. > > > > [[alternative HTML version deleted]] > > > > ___ > > R-sig-Geo mailing list > > R-sig-Geo@r-project.org > > > https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo=DwICAg=n6-cguzQvX_tUIrZOS_4Og=YFaRLkcUCdDkLrpTbNOUV9J1CwYBCTMwgm5tdQkRSm4=nytIIxO936Ls0xrd3zZkBd1WNjQB3DwlOK88GErq19M=uahwJjXmsZFnUVQXvkICr3EfAbNjOuaXl6iwziIexTM= > > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > > https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo=DwICAg=n6-cguzQvX_tUIrZOS_4Og=YFaRLkcUCdDkLrpTbNOUV9J1CwYBCTMwgm5tdQkRSm4=nytIIxO936Ls0xrd3zZkBd1WNjQB3DwlOK88GErq19M=uahwJjXmsZFnUVQXvkICr3EfAbNjOuaXl6iwziIexTM= > > > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Reverse y-axis with geom_sf
On 3/19/19 2:55 PM, Kent Johnson wrote: > Hi, > > I am working with data representing cells in a tissue sample with lines and > polygonal regions defined in the same space. The cells and polygons are > represented and manipulated as simple features objects and visualized using > ggplot2 and geom_sf. Mostly this works very well. The problem is that the > coordinate system of my data has the origin at the top left corner. For > data with no CRS, geom_sf puts the origin at the bottom left so all my > plots are inverted. > > Is there a way to invert the y axis while staying within sf (or possibly > sp?) The only answers I have found involve extracting the coordinates from > the sf objects and inverting or plotting from the raw data. Maybe by > creating a CRS with the correct orientation? > > For a very simple example with just a few vectors - the code below plots an > arrow pointing down; I would like to invert the y-axis so it points up. > library(sf) > library(ggplot2) > s1 <- rbind(c(9, 11), c(10, 10)) > s2 <- rbind(c(11, 11), c(10, 10)) > s3 <- rbind(c(10,14), c(10, 12), c(10,10)) > mls <- st_multilinestring(list(s1,s2,s3)) > > ggplot(mls) + geom_sf() You mean, like ggplot(mls * c(1, -1)) + geom_sf() ? > > Thank you for any help, > Kent Johnson > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] create a raster time series with negative years (years Before Present)
See ?as.Date: "Years before 1CE (aka 1AD) will probably not be handled correctly." On 2/22/19 3:59 PM, Jackson Rodrigues wrote: > Hey guys, > > I need some help to set a time serie (ts) to raster stack, however my > raster serie has negative ages (years Before Present (yrBP)). > > I have tried several codes and packages unsuccessfully. > > If succeed, I would like to apply seasonal and other calculations to my > new raster time serie. > > Below are some codes from package rts. I adapted years to negative value > only as a example. > > # >> path <- system.file("external", package="rts") # location of files >> lst <- list.files(path=path,pattern='.asc$',full.names=TRUE) >> lst # list of raster files >> r <- stack(lst) # creating a RasterStack object >> r >> d <- c("-2000-02-01","-2000-03-01","-2000-04-01","-2000-05-01") # > corresponding dates to 4 rasters >> d <- as.Date(d) # or d <- as.POSIXct(d) > Error in charToDate(x) : > character string is not in a standard unambiguous format >> rt <- rts(r,d) # creating a RasterStackTS object >> rt >> plot(rt) > # > > best, > > Jackson Rodrigues > -- > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Reverse color in stplot
I think that it accepts an argument col.regions where you specify a color palette; if you have a particular one, called pal, and want to revert it, use rev(pal). On 2/21/19 7:02 AM, Mahsa Sameti wrote: > Dear all > > I want to reverse color in my spatio-temporal maps that are created with > stplot function. How can i solve this problem? > > Thanks alot > Mahsa Sameti > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal segmentation fault during installation
"b") + nchar(sm[1L], type = > "b")if (w > LONG) prefix <- paste0(prefix, "\n ")} >else prefix <- "Error : "msg <- paste0(prefix, conditionMessage(e), > "\n").Internal(seterrmessage(msg[1L]))if (!silent && > isTRUE(getOption("show.error.messages"))) {cat(msg, file = outFile) >.Internal(printDeferredWarnings())}invisible(structure(msg, > class = "try-error", condition = e))}) > 20: try(suppressPackageStartupMessages(library(pkg_name, lib.loc = lib, > character.only = TRUE, logical.return = TRUE))) > 21: tools:::.test_load_package("rgdal", > "/home/ecor/R/x86_64-pc-linux-gnu-library/3.5") > An irrecoverable exception occurred. R is aborting now ... > Segmentation fault (core dumped) > ERROR: loading failed > * removing ‘/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/rgdal’ > > The downloaded source packages are in > ‘/tmp/Rtmpq9SibK/downloaded_packages’ > Warning message: > In install.packages("rgdal") : > installation of package ‘rgdal’ had non-zero exit status >> sessionInfo() > R version 3.5.1 (2018-07-02) > Platform: x86_64-pc-linux-gnu (64-bit) > Running under: Ubuntu 16.04.5 LTS > > Matrix products: default > BLAS: /usr/lib/libblas/libblas.so.3.6.0 > LAPACK: /usr/lib/lapack/liblapack.so.3.6.0 > > locale: > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > [3] LC_TIME=en_IE.UTF-8LC_COLLATE=en_US.UTF-8 > [5] LC_MONETARY=en_IE.UTF-8LC_MESSAGES=en_US.UTF-8 > [7] LC_PAPER=en_IE.UTF-8 LC_NAME=C > [9] LC_ADDRESS=C LC_TELEPHONE=C > [11] LC_MEASUREMENT=en_IE.UTF-8 LC_IDENTIFICATION=C > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > loaded via a namespace (and not attached): > [1] compiler_3.5.1 tools_3.5.1 >> > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] stars::RasterIO using extent info?
On 11/13/18 4:10 PM, Howard, Tim G (DEC) via R-sig-Geo wrote: > Dear list, > > I am exploring the different options for reading parts of large imagery > object in stars, as discussed here: > > https://r-spatial.github.io/stars/articles/proxy.html > > My ultimate goal is to read into RAM only a clipped portion of a large raster > (well, actually a raster stack, but taking baby steps here). > > My immediate question: the `RasterIO` option of read_stars defines cell > offsets and cell counts (*Size). Is there a straightforward way to calculate > these values given extent information? > > Reproducible example (mostly taken from here: > https://www.r-spatial.org/r/2018/03/22/stars2.html): > > library(stars) > tif <- system.file("tif/L7_ETMs.tif", package = "stars") > x <- read_stars(tif) # read entire tif into ram > x <- x[,,,1] #get just one layer for now > # calculate a circular polygon at the center of the raster > pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500) > plot(x) > # interestingly, I don't think the circle is in the right place when plotted > plot(st_geometry(pol), add = TRUE, border = "red") > # this is what I'd like to be able to restrict to what is read in memory: > plot(x[pol]) > > ## read only portion of tif using proxy object > x <- read_stars(tif, proxy = TRUE) > x <- x[,,,1] > y <- st_as_stars(x[pol]) > plot(y) # this is cropped to the extent (but not the circle - let's not worry > about that right now) > > Question: can I do the equivalent with the RasterIO options in stars? Said > another way, instead of setting up the proxy, can I map my extent object (or > bounding box) directly to the cell count values needed for RasterIO? stars can do the math, and so can you; it is explained here: https://r-spatial.github.io/stars/articles/data_model.html stars uses some functions directly from GDAL which it doesn't expose to the user, but there is no magic going on here. > > > Thanks in advance for any tips. > Tim > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Varying measurement error in Kriging predictions
Please read p 74, "Kriging data with known measurement errors", of http://gstat.org/gstat.pdf . It refers to a method published by Delhomme (1978). The software was written long before the Christensen paper you mention, to which I don't have access. On 11/6/18 3:55 PM, Antonis Alexiadis wrote: > Hello, my question is how can I implement known varying measurement error > in my > Kriging predictions? > > I have two separate datasets: the first one is a 2-D dataset that includes > the observations and > the second is the a 2-D dataset that includes the measurement error in each > specific location. > The measurement error dataset follows a structure, it's not random (which > can be characterized via the use of the experimental variogram). > > I have searched a lot online on how to incorporate the measurement > uncertainty in > Kriging predictions but it seems to still be an open question in the forums > (a paper solving this issue has been developed by William F. Christensen > titled as: Filtered Kriging for Spatial Data with Heterogeneous Measurement > Error Variances). I have tried to use gstat and incorporate the variances > using the weights functionality but after I do the kriging and visualize > the predicted variance field, even though it qualitatively resembles the > defined one, the values of the variances are magnitudes lower than the ones > proposed. > > Does anyone have an idea on how to solve it, or aware of some > software-package that > has already implemented this functionality? (It's quite tough to understand > the stated paper already, rather having to program > its contents) . > > > Thank you. > > Regards, > Antonios. > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Reordering geographical coordinates clockwise to make a polygon
Sorry for my previous answer, I didn't get the question fully. I'm afraid I still don't get the question fully, but functions that might help are sf::st_line_merge (creates a LINESTRING from line pieces) and sf::st_polygonize (creates a polygon from a LINESTRING that forms a closed ring) On 11/5/18 5:35 PM, Patrick Giraudoux wrote: > > Apologize to answer to myself: the way described below would work only > if there are no concave "peninsula" towards north or south inside the > polygon. Even thinking about an alternate solution, e.g. the Traveling > Salesman Problem (TSP), ie. the shortest route linking all points, one > could get wrong since points of a narrowpeninsula could have points > closest to the opposite border than from the next point on the border... > > Suppose we are stuck, and should redraw polygons by hand > > > Le 05/11/2018 à 14:02, Patrick Giraudoux a écrit : >> >> Dear listers, >> >> There is an interesting post here: >> https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order >> dealing on the issue. However, I would like to know if a function has >> been already developped in a R package. >> >> I am coping with a young colleague's shapefile where borders of >> polygons have been drawn as lines quite inconsistently. To change this >> SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is >> to select the points bordering each polygon (delete the others), >> define the point set as a Polygon, then rebuild step by step a >> SpatialPolygonsDataFrame with all its (~25) polygons. It would be much >> quicker than to redraw one by one the 5160 data points (two times on >> borders shared by two polygons). >> >> The problem is that the points must be reordered clockwise (the way >> lines making the borders is all except clockwise) before making them >> a polygon. >> >> Any function already developped for that ? >> >> Cheers, >> >> Patrick >> >> >> > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Reordering geographical coordinates clockwise to make a polygon
Look for argument check_ring_dir in the documentation of sf::st_read. On 11/5/18 2:02 PM, Patrick Giraudoux wrote: > Dear listers, > > There is an interesting post here: > https://stackoverflow.com/questions/6989100/sort-points-in-clockwise-order > dealing on the issue. However, I would like to know if a function has > been already developped in a R package. > > I am coping with a young colleague's shapefile where borders of polygons > have been drawn as lines quite inconsistently. To change this > SpatialLinesDataFrame into a SpatialPoints object is easy. The idea is > to select the points bordering each polygon (delete the others), > define the point set as a Polygon, then rebuild step by step a > SpatialPolygonsDataFrame with all its (~25) polygons. It would be much > quicker than to redraw one by one the 5160 data points (two times on > borders shared by two polygons). > > The problem is that the points must be reordered clockwise (the way > lines making the borders is all except clockwise) before making them a > polygon. > > Any function already developped for that ? > > Cheers, > > Patrick > > > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 pEpkey.asc Description: application/pgp-keys ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] www.asdar-book.org site issues
The website seems to be up again. On 08/11/2018 01:36 PM, Roger Bivand wrote: > The ASDAR book site is not being shown for viewing, but the content is > still there. Anyone needing content (several inquiries off-list) may try: > > ASDAR_BOOK <- "http://www.asdar-book.org/bundles2ed; > chapters <- c("hello", "cm", "vis", "die", "cm2", "std", "sppa", "geos", > "lat", "dismap") > for (i in chapters) {fn <- paste(i, "bundle.zip", sep="_"); > download.file(paste(ASDAR_BOOK, fn, sep = "/"), fn)} > > to download the chapter bundles. Summer maintenance takes time, and we > hope to fix the underlying problem in a week or so. > > Roger > -- Edzer Pebesma Institute for Geoinformatics Heisenbergstrasse 2, 48151 Muenster, Germany Phone: +49 251 8333081 ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] rgdal compile confusion
please report the entire output you see when trying to install rgdal. On 08/19/2018 08:14 PM, Rich Shepard wrote: > On Sat, 11 Aug 2018, Rich Shepard wrote: > > When I update.packages() I'm told that the attempt ended with a > non-zero error message: > > trying URL 'https://ftp.osuosl.org/pub/cran/src/contrib/rgdal_1.3-4.tar.gz' > Content type 'application/x-gzip' length 1664774 bytes (1.6 MB) > == > downloaded 1.6 MB > > Warning message: > In install.packages(update[instlib == l, "Package"], l, contriburl = > contriburl, : > installation of package ‘rgdal’ had non-zero exit status > > The update to version 1.3-3 last May worked. What might be the reason the > upgrade to 1.3-4 fails? > > R produces contradictory messages: > > configure: C++11 support available > configure: rgdal: 1.3-4 > ... > checking for gdal-config... /usr/bin/gdal-config > checking gdal-config usability... yes > configure: GDAL: 2.3.0 > checking C++11 support for GDAL >= 2.3.0... yes > checking GDAL version >= 1.11.4... yes > > Looks OK, yes? But, immediately following is this: > > In file included from /usr/include/gdal.h:45:0, > from gdal_test.cc:1: > /usr/include/cpl_port.h:187:6: error: #error Must have C++11 or newer. > # error Must have C++11 or newer. > ^ > ... > configure: Install failure: compilation and/or linkage problems. > configure: error: GDALAllRegister not found in libgdal. > ERROR: configuration failed for package ‘rgdal’ > * removing ‘/usr/lib/R/library/rgdal’ > * restoring previous ‘/usr/lib/R/library/rgdal’ > > Yet configure found C++11 support for gdal and the version of the > latter is > greater than 1.11.4. > > TIA, > > Rich > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] How to use gstat 'variogram' function for nonlinear formulae?
On 07/04/2018 06:14 PM, Joelle k. Akram wrote: > Dear all, > > > For a linear model I use the 'variogram' function in gstat as follows: > > > coordinates(data) <- c('x','y') > myvario <- variogram(dep~var1+var2+var3, data) > > Now, I have the following Nonlinear formula that I want to use in variogram. > > dep~ var1*var2^3 This is still a linear model. > > > My 2 questions are: > > i) Can the variogram function in gstat handle nonlinear models as this > formula? > > ii) What is the syntax to use this nonlinear formulae within the variogram > function? is it simply : > > myvario <- variogram(dep~ var1*var2^3, data) I think that like in other areas of R you'll have to use myvario <- variogram(dep~ var1*I(var2)^3, data) Note that the * will be interpreted as indicating you want main effects and an interaction, not simply the product of var1*var3^3; see ?formula > > Any advice is appeciated. > > > cheers > > Chris Akram > > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Difference in area calculation between QGIS and R
On 06/27/2018 05:29 PM, Michael Sumner wrote: > Raster uses a (discretized) cosine-of-latitude approximati (popular amongst > longlat map makers). As far as I can tell, raster branches into the libgeographic code (copied into the package sources, although Karney is not mentioned as one of the package copyright holders), using close to exact geodesic / ellipsoidal computation. > > QGIS uses a project to local equal area projection method or maybe some > other approach. > > There's lots and f options, all that matters is what your work needs. > > Cheers, Mike > > On Wed, 27 Jun 2018, 22:14 Suncus Etruscus, > wrote: > >> Dear List, >> I am trying to use R ("sp" and "raster" packages) to calculate the area of >> several polygons (CRS of the shapefile EPSG: 4326 - WGS84). >> >> I used this line of code: >> >> shapefile_name$area_km2 <- area(shapefile_name)/100 >> >> However, when I used QGIS to calculate the area of the same polygons >> (through the "$area" function), I found there was a slight difference (in >> every polygon). >> >> For example, for a polygon of about 30 000 km2, the area calculated in R >> was 50 km2 smaller. >> >> What could be the cause? >> >> Thanks in advance, >> N. >> >> [[alternative HTML version deleted]] >> >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] What is the Coordinate Reference System?
Mb >>> wrf <- nc_open(a) >>> paste("The file has",wrf$nvars,"variables") >>> paste("The file has",names(wrf$var)) >>> lat <- ncvar_get( wrf, "XLAT" ) >>> lon <- ncvar_get( wrf, "XLONG" ) >>> t2 <- ncvar_get( wrf, "T2" ) - 273.15 >>> dim(t2) #73 60 13 >>> #image >>> image(t2[, , 12]) >>> # raster >>> rt2 <- raster(t(t2[1:dim(t2)[1], >>> dim(t2)[2]:1, >>> 12]), >>> xmn = min(lon), >>> xmx = max(lon), >>> ymn = min(lat), >>> ymx = max(lat), >>> * crs="+init=epsg:4326") # ???* >>> spplot(rt2) >>> >>> >>> >>> -- >>> Sergio Ibarra Espinosa >>> Post Doc >>> PhD in Atmospheric Sciences >>> Instituto de Astronomia, Geofísica e Ciências Atmosféricas >>> Universidade de São Paulo >>> Rua do Matão, 1226 >>> São Paulo-SP - Brasil - >>> 05508-090 >>> +55-11-97425-3791 >>> Skype: sergio_ibarra1 >>> >>>[[alternative HTML version deleted]] >>> >>> ___ >>> R-sig-Geo mailing list >>> R-sig-Geo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >>> >> -- >> Dr. Michael Sumner >> Software and Database Engineer >> Australian Antarctic Division >> 203 Channel Highway >> Kingston Tasmania 7050 Australia >> >> [[alternative HTML version deleted]] >> >> ___ >> R-sig-Geo mailing list >> R-sig-Geo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo >> > > Ben Tupper > Bigelow Laboratory for Ocean Sciences > 60 Bigelow Drive, P.O. Box 380 > East Boothbay, Maine 04544 > http://www.bigelow.org > > Ecological Forecasting: https://eco.bigelow.org/ > > > > > > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > 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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo