Re: [R-sig-Geo] as.im in spatstat::ppm not working as before

2024-05-11 Thread Edzer Pebesma



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

2024-05-08 Thread Edzer Pebesma
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

2024-02-28 Thread Edzer Pebesma




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

2023-07-13 Thread Edzer Pebesma

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

2023-02-10 Thread Edzer Pebesma
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?

2023-01-31 Thread Edzer Pebesma
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?

2023-01-31 Thread Edzer Pebesma




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?

2023-01-25 Thread Edzer Pebesma




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?

2023-01-25 Thread Edzer Pebesma
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

2023-01-25 Thread Edzer Pebesma




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

2022-12-14 Thread Edzer Pebesma
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

2022-11-29 Thread Edzer Pebesma

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

2022-10-22 Thread Edzer Pebesma

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

2022-10-12 Thread Edzer Pebesma

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

2022-10-07 Thread Edzer Pebesma
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

2022-04-29 Thread Edzer Pebesma

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

2021-11-14 Thread Edzer Pebesma



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

2021-11-14 Thread Edzer Pebesma
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

2021-11-14 Thread Edzer Pebesma
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

2021-11-02 Thread Edzer Pebesma

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

2021-10-14 Thread Edzer Pebesma
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

2021-09-14 Thread Edzer Pebesma
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

2021-09-10 Thread Edzer Pebesma
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

2021-07-22 Thread Edzer Pebesma
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

2021-06-21 Thread Edzer Pebesma

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

2021-06-07 Thread Edzer Pebesma




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

2021-06-01 Thread Edzer Pebesma





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

2021-06-01 Thread Edzer Pebesma




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

2021-04-01 Thread Edzer Pebesma




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)

2021-03-18 Thread Edzer Pebesma

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)

2021-03-15 Thread Edzer Pebesma

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

2021-02-09 Thread Edzer Pebesma
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

2021-02-09 Thread Edzer Pebesma

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

2021-02-04 Thread Edzer Pebesma




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

2021-02-03 Thread Edzer Pebesma
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

2020-12-29 Thread Edzer Pebesma
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

2020-12-01 Thread Edzer Pebesma
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.

2020-10-31 Thread Edzer Pebesma




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

2020-09-14 Thread Edzer Pebesma
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

2020-09-09 Thread Edzer Pebesma
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

2020-09-08 Thread Edzer Pebesma




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

2020-06-20 Thread Edzer Pebesma


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

2020-06-06 Thread Edzer Pebesma


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

2020-06-06 Thread Edzer Pebesma
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

2020-06-03 Thread Edzer Pebesma
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

2020-04-29 Thread Edzer Pebesma
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)

2020-03-24 Thread Edzer Pebesma
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

2020-03-24 Thread Edzer Pebesma
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.

2020-03-19 Thread Edzer Pebesma
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

2020-03-18 Thread Edzer Pebesma
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.

2020-03-18 Thread Edzer Pebesma
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

2020-02-19 Thread Edzer Pebesma
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()?

2020-01-21 Thread Edzer Pebesma
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

2019-12-12 Thread Edzer Pebesma
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

2019-12-06 Thread Edzer Pebesma
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

2019-10-29 Thread Edzer Pebesma
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

2019-10-28 Thread Edzer Pebesma
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

2019-10-26 Thread Edzer Pebesma


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

2019-10-25 Thread Edzer Pebesma
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

2019-10-15 Thread Edzer Pebesma
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

2019-10-15 Thread Edzer Pebesma
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

2019-10-07 Thread Edzer Pebesma
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

2019-10-01 Thread Edzer Pebesma
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

2019-09-12 Thread Edzer Pebesma
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

2019-09-12 Thread Edzer Pebesma


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

2019-09-11 Thread Edzer Pebesma
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?

2019-09-04 Thread Edzer Pebesma
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

2019-08-22 Thread Edzer Pebesma
 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

2019-08-21 Thread Edzer Pebesma
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)

2019-08-15 Thread 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

2019-08-15 Thread Edzer Pebesma


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

2019-08-10 Thread Edzer Pebesma
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

2019-08-10 Thread Edzer Pebesma
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

2019-07-04 Thread Edzer Pebesma
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

2019-06-20 Thread Edzer Pebesma


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

2019-06-11 Thread Edzer Pebesma
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

2019-06-04 Thread Edzer Pebesma
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?

2019-05-27 Thread Edzer Pebesma
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

2019-05-12 Thread Edzer Pebesma


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?

2019-05-07 Thread Edzer Pebesma


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

2019-05-06 Thread Edzer Pebesma


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

2019-05-06 Thread Edzer Pebesma
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

2019-05-04 Thread Edzer Pebesma
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

2019-05-03 Thread Edzer Pebesma
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

2019-04-30 Thread Edzer Pebesma
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

2019-04-12 Thread Edzer Pebesma
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

2019-04-11 Thread Edzer Pebesma
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

2019-03-19 Thread Edzer Pebesma


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)

2019-02-22 Thread Edzer Pebesma
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

2019-02-20 Thread Edzer Pebesma
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

2018-11-17 Thread Edzer Pebesma
 "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?

2018-11-13 Thread Edzer Pebesma


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

2018-11-07 Thread Edzer Pebesma
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

2018-11-05 Thread Edzer Pebesma
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

2018-11-05 Thread Edzer Pebesma
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

2018-08-21 Thread Edzer Pebesma
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

2018-08-19 Thread Edzer Pebesma
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?

2018-07-04 Thread Edzer Pebesma



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

2018-06-27 Thread Edzer Pebesma


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?

2018-06-18 Thread Edzer Pebesma
 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


  1   2   3   4   5   6   7   8   >