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

Unable to use factors in raster::predict

11 hours 43 min ago
Dear list,

I would like to use the 'predict' function in the 'raster' package in an implementation of species distribution modelling with a couple of factor variables. My case can be set up exactly as the cforest example listed in the help file. Unfortunately, I cannot get the example to work: it throws an error because of the unused argument 'factors', which I need. Here is the code I am referring to:

# create a RasterStack or RasterBrick with with a set of predictor layers
logo <- brick(system.file("external/rlogo.grd", package="raster"))
names(logo)

# known presence and absence points
p <- matrix(c(48, 48, 48, 53, 50, 46, 54, 70, 84, 85, 74, 84, 95, 85,
              66, 42, 26, 4, 19, 17, 7, 14, 26, 29, 39, 45, 51, 56, 46, 38, 31,
              22, 34, 60, 70, 73, 63, 46, 43, 28), ncol=2)

a <- matrix(c(22, 33, 64, 85, 92, 94, 59, 27, 30, 64, 60, 33, 31, 9,
              99, 67, 15, 5, 4, 30, 8, 37, 42, 27, 19, 69, 60, 73, 3, 5, 21,
              37, 52, 70, 74, 9, 13, 4, 17, 47), ncol=2)

# extract values for points
xy <- rbind(cbind(1, p), cbind(0, a))
v <- data.frame(cbind(pa=xy[,1], extract(logo, xy[,2:3])))

# cforest (other Random Forest implementation) example with factors argument

v$red <- as.factor(round(v$red/100))
logo$red <- round(logo[[1]]/100)

library(party)
m <- cforest(pa~., control=cforest_unbiased(mtry=3), data=v)
f <- list(levels(v$red))
names(f) <- 'red'
pc <- predict( m,logo, OOB=TRUE, factors=f)

Error in RET@prediction_weights(newdata = newdata, mincriterion = mincriterion,  :
  unused argument (factors = f)


Note that I needed to change the order in the arguments in the last line. If I run it the way it was in the example, like this:

pc <- predict(logo, m, OOB=TRUE, factors=f)


then the error is the following:


Error in v[cells, ] <- predv :

  number of items to replace is not a multiple of replacement length

Also note that if I run the line without the 'factors' argument the resulting value (r) is a matrix, and not a raster object.

I am using Package 'raster' version 2.6-7 on Windows 7 (where this worked fine in the past).

Many thanks for any help,

Genoveva

Genoveva Gonzalez Mirelis, Scientist
Institute of Marine Research
Nordnesgaten 50
5005 Bergen, Norway
Phone number +47 55238510


        [[alternative HTML version deleted]]

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

Dealing with non-nested geographic units using the RQGIS and tigris R packages

Mon, 11/19/2018 - 09:34
I recently came across this tutorial suggesting that the RQGIS and tigris R
packages could be used to overcome issues with non-nested geographies:

https://rpubs.com/corey_sparks/256942

I'm trying to apply this to some data I have where I would like to look at
Census Block Groups nested within Zip Code Tabulation Areas. Specifically,
I'd like to compute a segregation index using Census Block Groups as the
lower units and Zip Code Tabulation Areas as the higher level units.

I am confused about how these R packages deal with this issue of non-nested
units. Here's the explanation from the tutorial:

"In order to get around this problem [non-nested units], we can use the
QGIS application to perform a geographic intersection between two spatial
polygon layers that we want to use. This will allow us to merge the
population data at the census tract level to a higher level, so we can
calculate a segregation index."

Can anyone say more about how this is possible? I can run the code in the
tutorial but I don't understand what it's doing.

Thanks!

        [[alternative HTML version deleted]]

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

set up months to a time serie

Mon, 11/19/2018 - 05:26
Dear all,

I am trying to set a sequence of  48000 months to a simulation between
10200 - 6201  before present.

However I always get the following message
Error in charToDate(TS) : #character string is not in a standard
unambiguous format
Could anyone help me?
####
library(zoo)
TS1<-brick("TS1.nc")
# dimensions  : 12, 12, 144, 48000  (nrow, ncol, ncell, nlayers)
TS1 <- as.yearmon(10200 - seq(0, 48000)/12)

####

Moreover, I would like to detect changes inside my rasterbrick serie, I
read that package bfast can detect changes in such of time series and
files, anyone experienced with bfast could help me with it as well?

thank you all,

Best

Jackson

        [[alternative HTML version deleted]]

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

Re: Pixel calculation using extract() in geoTiff raster: faster and free of computer lock in R

Sun, 11/18/2018 - 10:02
Hi Alexandre,

it seems to me this is a job for function 'focal' rather then 'extract'.
Even if so, I'm not sure about the speed it would take.

Cheers
Hugo

ASANTOS via R-sig-Geo <[hidden email]> escreveu no dia domingo,
18/11/2018 à(s) 14:10:

> Dear R-sig-geo Members,
>
>      I've like to use R and raster package (extract() function) for
> pixel calculations (SD, Skewness and Kurtosis in 1  unit radius of a
> buffer) in a geoTiff raster, but my original geoTiff image has
> 6244(nrow), 8721(ncol) and 54453924 (ncell) dimensions. This image size
> causes slow computation process or computer lock, for example in my code
> below:
>
> ### <code r>
> #Packages
> library(raster)
> library(moments) #Measures of Skewness and Kurtosis
>
> ## Create artificial raster for my geoTiff simulation - dimensions  of
> my original geoTiff: 6244, 8721, 54453924  (nrow, ncol, ncell)
> r <- raster(nc=8721, nr=6244)
> r <- setValues(r, round(runif(ncell(r))* 255))
>
> #Extract all pixel coordinates in raster
> coord_r<-coordinates(r)
>
> #Extract standard deviation, skewness and kurtosis
> Buffer<-1
> SD<-function (x, na.rm = TRUE)
> {
> if (is.matrix(x))
>      apply(x, 2, sd, na.rm = na.rm)
> else if (is.vector(x))
>      sqrt(var(x, na.rm = na.rm))
> else if (is.data.frame(x))
>      sapply(x, sd, na.rm = na.rm)
> else sqrt(var(as.vector(x), na.rm = na.rm))
> }
> desv_pad_R<-extract(r, coord_r, buffer = Buffer, fun = SD)
> str(desv_pad_R)
> sk_R <-extract(r,coord_r,buffer=Buffer, fun=skewness, na.rm = TRUE)
> str(sk_R)
> k_R <-extract(r,coord_r,buffer=Buffer, fun=kurtosis, na.rm = TRUE)
> str(k_R)
> # <END code>
>
>      There are different approaches (eg. integration with SAGA GIS or
> GRASS, convert image to ASCII) for my problem?
>
> Thanks in advanced,
>
> Alexandre
>
> --
> ======================================================================
> Alexandre dos Santos
> Proteção Florestal
> IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
> Campus Cáceres
> Caixa Postal 244
> Avenida dos Ramires, s/n
> Bairro: Distrito Industrial
> Cáceres - MT                      CEP: 78.200-000
> Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)
>
>          [hidden email]
> Lattes: http://lattes.cnpq.br/1360403201088680
> OrcID: orcid.org/0000-0001-8232-6722
> Researchgate: www.researchgate.net/profile/Alexandre_Santos10
> LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
> Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
        [[alternative HTML version deleted]]

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

Pixel calculation using extract() in geoTiff raster: faster and free of computer lock in R

Sun, 11/18/2018 - 07:10
Dear R-sig-geo Members,

     I've like to use R and raster package (extract() function) for
pixel calculations (SD, Skewness and Kurtosis in 1  unit radius of a
buffer) in a geoTiff raster, but my original geoTiff image has
6244(nrow), 8721(ncol) and 54453924 (ncell) dimensions. This image size
causes slow computation process or computer lock, for example in my code
below:

### <code r>
#Packages
library(raster)
library(moments) #Measures of Skewness and Kurtosis

## Create artificial raster for my geoTiff simulation - dimensions  of
my original geoTiff: 6244, 8721, 54453924  (nrow, ncol, ncell)
r <- raster(nc=8721, nr=6244)
r <- setValues(r, round(runif(ncell(r))* 255))

#Extract all pixel coordinates in raster
coord_r<-coordinates(r)

#Extract standard deviation, skewness and kurtosis
Buffer<-1
SD<-function (x, na.rm = TRUE)
{
if (is.matrix(x))
     apply(x, 2, sd, na.rm = na.rm)
else if (is.vector(x))
     sqrt(var(x, na.rm = na.rm))
else if (is.data.frame(x))
     sapply(x, sd, na.rm = na.rm)
else sqrt(var(as.vector(x), na.rm = na.rm))
}
desv_pad_R<-extract(r, coord_r, buffer = Buffer, fun = SD)
str(desv_pad_R)
sk_R <-extract(r,coord_r,buffer=Buffer, fun=skewness, na.rm = TRUE)
str(sk_R)
k_R <-extract(r,coord_r,buffer=Buffer, fun=kurtosis, na.rm = TRUE)
str(k_R)
# <END code>

     There are different approaches (eg. integration with SAGA GIS or
GRASS, convert image to ASCII) for my problem?

Thanks in advanced,

Alexandre

--
======================================================================
Alexandre dos Santos
Proteção Florestal
IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
Campus Cáceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
Cáceres - MT                      CEP: 78.200-000
Fone: (+55) 65 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

         [hidden email]
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/

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

Re: logLik value of the spml spatial panel model

Sat, 11/17/2018 - 13:59
Dear Dr. Piras: This is great. Thank you very much for the response.

Best,

Danlin


On 2018/11/14 16:32, Gianfranco Piras wrote:
> Sorry for the very slow reply. I was  looking back at the email exchange and this problem sounded familiar.
>
> I have had another user contacting me because he was experiencing the same issue with the pooled model. He also compared the value of the likelihood from R with the value that he was obtaining with codes in Matlab (in the spatial econometrics toolbox).  Anyway, Roger was right when he said that the log-likelihood  from the "fix” was calculated omitting constants. I fixed this specific problem a long time ago, but since the email that was sent to me by the other user had several  other questions (which I didn’t had time to address),  I never committed the new version on R-forge.
>
> Anyway, here is a script and the updated package. I will try to get to the other questions asap and finally upload to CRAN.
>
> Sorry  for any inconvenience related to my delay!
>
> Best,
> G.
>
>
>
>
>      
>
>
>
>> On Oct 25, 2018, at 10:56 AM, Millo Giovanni <[hidden email]> wrote:
>>
>> Correct. It should not be used for comparing FE to anything else. In particular, comparisons between FE and RE specifications have to concentrate on the treatment of the individual effects.
>>
>> Best, Giovanni
>>
>> -----Messaggio originale-----
>> Da: Danlin Yu [mailto:[hidden email]]
>> Inviato: giovedì 25 ottobre 2018 16:54
>> A: Millo Giovanni; Jeremie Juste; Roger Bivand
>> Cc: [hidden email]; Gianfranco Piras
>> Oggetto: Re: R: [R-sig-Geo] logLik value of the spml spatial panel model
>>
>> This is great response. Thank you all, Giovanni, Roger and Jeremie.
>> Based on the response, it seems that the LL element for fixed effect
>> model is only useful for optimization, but not cross model comparison?
>> The fixed effect spatial lag panel model, however, does produce the LL
>> value. Does that mean this LL value shall not be used for model
>> comparison purpose as well?
>>
>> Thank you all.
>>
>> Best,
>>
>> Danlin
>>
>>
>> On 2018/10/25 8:46, Millo Giovanni wrote:
>>> Hello. Thanks for keeping us posted.
>>>
>>> As was aptly observed, while RE and pooled models do actually have a logLik element, FE ones don't but there's a reason. While it is not difficult to output the  LL for the FE models too, like Jeremie did (or perhaps even setting quiet=FALSE and taking the last loglikelihood value from optimization), we never managed to decide whether outputting this is appropriate or not. In fact the LL element from the procedure is relative to the "pooled" model on demeaned data, so it is sort of a crossbred likelihood; most people asking for it want it to compare to the LL of the random effects model or similar uses, which I would not advise.
>>>
>>> Gianfranco (author of the FE stuff) do please step in if appropriate.
>>> Best,
>>> Giovanni
>>>
>>> -----Messaggio originale-----
>>> Da: Jeremie Juste [mailto:[hidden email]]
>>> Inviato: giovedì 25 ottobre 2018 11:21
>>> A: Roger Bivand
>>> Cc: Danlin Yu; [hidden email]; Millo Giovanni; Gianfranco Piras
>>> Oggetto: Re: [R-sig-Geo] logLik value of the spml spatial panel model
>>>
>>>
>>> Hello,
>>>
>>> Thanks for the feedback.
>>>
>>> Here is the diff. It is my first one so I'm not sure I'm doing everything right. :-)
>>>
>>>
>>> modified   R/likelihoodsFE.R
>>> @@ -459,7 +459,7 @@ if(Hess) asyv <- NULL  else asyv <- asyv
>>>
>>>
>>> -       return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se, rho.se = rho.se, s2.se = s2.se, asyvar1=asyvar1, residuals = r, asyv = asyv)
>>> +    return<-list(coeff=betas, rho = rho, s2 = s2, rest.se = rest.se,
>>> + rho.se = rho.se, s2.se = s2.se, ll=LL, asyvar1=asyvar1, residuals = r,
>>> + asyv = asyv)
>>>   }
>>>
>>> I have compressed splm with the modification with the .git inside. I have installed the modified package again and checked that the modification worked.
>>>
>>>> fespaterr <- spml(fm, data = Produc, listw = mat2listw(usaww),
>>>> model="within", spatial.error="b", Hess = FALSE) fespaterr$logLik
>>>> It may be,however, that there was a reason for omitting it, so checks
>>>> to verify that it actually is a logLik value on the same basis as
>>>> other logLik values (not from an abbreviated log-likelihood function
>>>> used for fitting by omitting constants) would be helpful. It would be
>>>> good to provide logLik methods for the ML fitted model objects.
>>> I have checked that it is the logLik that correspond to the estimation of the rho.
>>>
>>> Best regards,
>>>
>>> Jeremie
>>>
>>>
>>>
>>>
>>> Ai sensi del Reg. UE 2016/679 e normativa vigente si precisa che le informazioni contenute in questo messaggio sono riservate ed a uso esclusivo del destinatario. Qualora il messaggio in parola Le fosse pervenuto per errore, La invitiamo ad eliminarlo senza copiarlo e a non inoltrarlo a terzi, dandocene gentilmente comunicazione. Grazie.
>>>
>>> Pursuant to Italian Law and EU Reg. 2016/679, you are hereby informed that this message contains confidential information intended only for the use of the addressee. If you are not the addressee, and have received this message by mistake, please delete it and immediately notify us. You may not copy or disseminate this message to anyone. Thank you.
>> --
>> ___________________________________________
>> Danlin Yu, Ph.D.
>> Professor of GIS and Urban Geography
>> Department of Earth & Environmental Studies
>> Montclair State University
>> Montclair, NJ, 07043
>> Tel: 973-655-4313
>> Fax: 973-655-4072
>> Office: CELS 314
>> Email: [hidden email]
>> webpage: csam.montclair.edu/~yu
>>
--
___________________________________________
Danlin Yu, Ph.D.
Professor of GIS and Urban Geography
Department of Earth & Environmental Studies
Montclair State University
Montclair, NJ, 07043
Tel: 973-655-4313
Fax: 973-655-4072
Office: CELS 314
Email: [hidden email]
webpage: csam.montclair.edu/~yu


        [[alternative HTML version deleted]]

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

Re: rgdal segmentation fault during installation

Sat, 11/17/2018 - 07:04
It looks like you have multiple versions of GDAL and/or PROJ installed
on your system; see https://github.com/r-spatial/sf the section on
"Multiple GDAL, GEOS and/or PROJ versions on your system" and there the
link to https://github.com/r-spatial/sf/issues/844 which contains a very
ugly solution (setting LD_LIBRARY_PATH from within R).

On 11/17/18 1:57 PM, Emanuele Cordano wrote:
> Dear list,
>
> I got the following segmentation fault error when I try to install rgdal on
> my local Ubuntu 16.04 LTS machine. Did anyone encounter something similar?
> Please below the error message with traceback and R session info.
> Similar problems appears on
> https://stackoverflow.com/questions/50064673/cannot-install-rgdal-package-in-local-server-an-irrecoverable-exception-occurre
> which remains unsolved . In my case (see below)  all checks, PROJ4
> included,  were OK!
>
> Thank you for any kind of feedback
> Best
> Emanuele Cordano
>
> user@machine:~$ R
>
> R version 3.5.1 (2018-07-02) -- "Feather Spray"
> Copyright (C) 2018 The R Foundation for Statistical Computing
> Platform: x86_64-pc-linux-gnu (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.
>
> [Previously saved workspace restored]
>
>> install.packages("rgdal")
> Installing package into ‘/home/ecor/R/x86_64-pc-linux-gnu-library/3.5’
> (as ‘lib’ is unspecified)
> trying URL 'http://cran.r-project.org/src/contrib/rgdal_1.3-6.tar.gz'
> Content type 'application/x-gzip' length 1666975 bytes (1.6 MB)
> ==================================================
> downloaded 1.6 MB
>
> * installing *source* package ‘rgdal’ ...
> ** package ‘rgdal’ successfully unpacked and MD5 sums checked
> configure: R_HOME: /usr/lib/R
> configure: CC: gcc -std=gnu99
> configure: CXX: g++
> configure: C++11 support available
> configure: rgdal: 1.3-6
> checking for /usr/bin/svnversion... yes
> configure: svn revision: 773
> checking for gdal-config... /usr/bin/gdal-config
> checking gdal-config usability... yes
> configure: GDAL: 2.1.3
> checking GDAL version >= 1.11.4... yes
> checking gdal: linking with --libs only... yes
> checking GDAL: /usr/share/gdal/2.1/pcs.csv readable... yes
> configure: pkg-config proj exists, will use it
> configure: PROJ version: 4.9.2
> checking proj_api.h presence and usability... yes
> checking PROJ version >= 4.8.0... yes
> checking projects.h presence and usability... yes
> checking PROJ.4: epsg found and readable... yes
> checking PROJ.4: conus found and readable... yes
> configure: Package CPP flags:  -I/usr/include/gdal
> configure: Package LIBS:  -L/usr/lib -lgdal -lproj
> configure: creating ./config.status
> config.status: creating src/Makevars
> ** libs
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c OGR_write.cpp -o OGR_write.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c gdal-bindings.cpp -o gdal-bindings.o
> gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g  -c init.c -o init.o
> gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g  -c inverser.c -o inverser.o
> gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g  -c local_stubs.c -o local_stubs.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c ogr_geom.cpp -o ogr_geom.o
> gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g  -c ogr_polygons.c -o ogr_polygons.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c ogr_proj.cpp -o ogr_proj.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c ogrdrivers.cpp -o ogrdrivers.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c ogrsource.cpp -o ogrsource.o
> g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
> -I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
> -O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
> -D_FORTIFY_SOURCE=2 -g -c projectit.cpp -o projectit.o
> g++ -std=gnu++11 -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions
> -Wl,-z,relro -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
> projectit.o -L/usr/lib -lgdal -lproj -L/usr/lib/R/lib -lR
> installing to /home/ecor/R/x86_64-pc-linux-gnu-library/3.5/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
>
>  *** caught segfault ***
> address 0x20, cause 'memory not mapped'
>
> Traceback:
>  1: assign("has_proj_def.dat", .Call("PROJ4_proj_def_dat_Installed",
> PACKAGE = "rgdal"), envir = .RGDAL_CACHE)
>  2: fun(libname, pkgname)
>  3: doTryCatch(return(expr), name, parentenv, handler)
>  4: tryCatchOne(expr, names, parentenv, handlers[[1L]])
>  5: tryCatchList(expr, classes, parentenv, handlers)
>  6: tryCatch(fun(libname, pkgname), error = identity)
>  7: runHook(".onLoad", env, package.lib, package)
>  8: loadNamespace(package, lib.loc)
>  9: doTryCatch(return(expr), name, parentenv, handler)
> 10: tryCatchOne(expr, names, parentenv, handlers[[1L]])
> 11: tryCatchList(expr, classes, parentenv, handlers)
> 12: tryCatch({    attr(package, "LibPath") <- which.lib.loc    ns <-
> loadNamespace(package, lib.loc)    env <- attachNamespace(ns, pos = pos,
> deps)}, error = function(e) {    P <- if (!is.null(cc <- conditionCall(e)))
>         paste(" in", deparse(cc)[1L])    else ""    msg <-
> gettextf("package or namespace load failed for %s%s:\n %s",
> sQuote(package), P, conditionMessage(e))    if (logical.return)
> message(paste("Error:", msg), domain = NA)    else stop(msg, call. = FALSE,
> domain = NA)})
> 13: library(pkg_name, lib.loc = lib, character.only = TRUE, logical.return
> = TRUE)
> 14: withCallingHandlers(expr, packageStartupMessage = function(c)
> invokeRestart("muffleMessage"))
> 15: suppressPackageStartupMessages(library(pkg_name, lib.loc = lib,
> character.only = TRUE, logical.return = TRUE))
> 16: doTryCatch(return(expr), name, parentenv, handler)
> 17: tryCatchOne(expr, names, parentenv, handlers[[1L]])
> 18: tryCatchList(expr, classes, parentenv, handlers)
> 19: tryCatch(expr, error = function(e) {    call <- conditionCall(e)    if
> (!is.null(call)) {        if (identical(call[[1L]], quote(doTryCatch)))
>         call <- sys.call(-4L)        dcall <- deparse(call)[1L]
>  prefix <- paste("Error in", dcall, ": ")        LONG <- 75L        sm <-
> strsplit(conditionMessage(e), "\n")[[1L]]        w <- 14L + nchar(dcall,
> type = "w") + nchar(sm[1L], type = "w")        if (is.na(w))             w
> <- 14L + nchar(dcall, type = "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-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_IE.UTF-8    LC_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

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

pEpkey.asc (2K) Download Attachment

rgdal segmentation fault during installation

Sat, 11/17/2018 - 06:57
Dear list,

I got the following segmentation fault error when I try to install rgdal on
my local Ubuntu 16.04 LTS machine. Did anyone encounter something similar?
Please below the error message with traceback and R session info.
Similar problems appears on
https://stackoverflow.com/questions/50064673/cannot-install-rgdal-package-in-local-server-an-irrecoverable-exception-occurre
which remains unsolved . In my case (see below)  all checks, PROJ4
included,  were OK!

Thank you for any kind of feedback
Best
Emanuele Cordano

user@machine:~$ R

R version 3.5.1 (2018-07-02) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (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.

[Previously saved workspace restored]

> install.packages("rgdal")
Installing package into ‘/home/ecor/R/x86_64-pc-linux-gnu-library/3.5’
(as ‘lib’ is unspecified)
trying URL 'http://cran.r-project.org/src/contrib/rgdal_1.3-6.tar.gz'
Content type 'application/x-gzip' length 1666975 bytes (1.6 MB)
==================================================
downloaded 1.6 MB

* installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
configure: R_HOME: /usr/lib/R
configure: CC: gcc -std=gnu99
configure: CXX: g++
configure: C++11 support available
configure: rgdal: 1.3-6
checking for /usr/bin/svnversion... yes
configure: svn revision: 773
checking for gdal-config... /usr/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 2.1.3
checking GDAL version >= 1.11.4... yes
checking gdal: linking with --libs only... yes
checking GDAL: /usr/share/gdal/2.1/pcs.csv readable... yes
configure: pkg-config proj exists, will use it
configure: PROJ version: 4.9.2
checking proj_api.h presence and usability... yes
checking PROJ version >= 4.8.0... yes
checking projects.h presence and usability... yes
checking PROJ.4: epsg found and readable... yes
checking PROJ.4: conus found and readable... yes
configure: Package CPP flags:  -I/usr/include/gdal
configure: Package LIBS:  -L/usr/lib -lgdal -lproj
configure: creating ./config.status
config.status: creating src/Makevars
** libs
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c OGR_write.cpp -o OGR_write.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c gdal-bindings.cpp -o gdal-bindings.o
gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g  -c init.c -o init.o
gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g  -c inverser.c -o inverser.o
gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g  -c local_stubs.c -o local_stubs.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c ogr_geom.cpp -o ogr_geom.o
gcc -std=gnu99 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g  -c ogr_polygons.c -o ogr_polygons.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c ogr_proj.cpp -o ogr_proj.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c ogrdrivers.cpp -o ogrdrivers.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c ogrsource.cpp -o ogrsource.o
g++ -std=gnu++11 -I"/usr/share/R/include" -DNDEBUG -I/usr/include/gdal
-I"/home/ecor/R/x86_64-pc-linux-gnu-library/3.5/sp/include"    -fpic  -g
-O2 -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time
-D_FORTIFY_SOURCE=2 -g -c projectit.cpp -o projectit.o
g++ -std=gnu++11 -shared -L/usr/lib/R/lib -Wl,-Bsymbolic-functions
-Wl,-z,relro -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
projectit.o -L/usr/lib -lgdal -lproj -L/usr/lib/R/lib -lR
installing to /home/ecor/R/x86_64-pc-linux-gnu-library/3.5/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

 *** caught segfault ***
address 0x20, cause 'memory not mapped'

Traceback:
 1: assign("has_proj_def.dat", .Call("PROJ4_proj_def_dat_Installed",
PACKAGE = "rgdal"), envir = .RGDAL_CACHE)
 2: fun(libname, pkgname)
 3: doTryCatch(return(expr), name, parentenv, handler)
 4: tryCatchOne(expr, names, parentenv, handlers[[1L]])
 5: tryCatchList(expr, classes, parentenv, handlers)
 6: tryCatch(fun(libname, pkgname), error = identity)
 7: runHook(".onLoad", env, package.lib, package)
 8: loadNamespace(package, lib.loc)
 9: doTryCatch(return(expr), name, parentenv, handler)
10: tryCatchOne(expr, names, parentenv, handlers[[1L]])
11: tryCatchList(expr, classes, parentenv, handlers)
12: tryCatch({    attr(package, "LibPath") <- which.lib.loc    ns <-
loadNamespace(package, lib.loc)    env <- attachNamespace(ns, pos = pos,
deps)}, error = function(e) {    P <- if (!is.null(cc <- conditionCall(e)))
        paste(" in", deparse(cc)[1L])    else ""    msg <-
gettextf("package or namespace load failed for %s%s:\n %s",
sQuote(package), P, conditionMessage(e))    if (logical.return)
message(paste("Error:", msg), domain = NA)    else stop(msg, call. = FALSE,
domain = NA)})
13: library(pkg_name, lib.loc = lib, character.only = TRUE, logical.return
= TRUE)
14: withCallingHandlers(expr, packageStartupMessage = function(c)
invokeRestart("muffleMessage"))
15: suppressPackageStartupMessages(library(pkg_name, lib.loc = lib,
character.only = TRUE, logical.return = TRUE))
16: doTryCatch(return(expr), name, parentenv, handler)
17: tryCatchOne(expr, names, parentenv, handlers[[1L]])
18: tryCatchList(expr, classes, parentenv, handlers)
19: tryCatch(expr, error = function(e) {    call <- conditionCall(e)    if
(!is.null(call)) {        if (identical(call[[1L]], quote(doTryCatch)))
        call <- sys.call(-4L)        dcall <- deparse(call)[1L]
 prefix <- paste("Error in", dcall, ": ")        LONG <- 75L        sm <-
strsplit(conditionMessage(e), "\n")[[1L]]        w <- 14L + nchar(dcall,
type = "w") + nchar(sm[1L], type = "w")        if (is.na(w))             w
<- 14L + nchar(dcall, type = "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-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_IE.UTF-8    LC_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
>

--
Emanuele Cordano, PhD
Environmental Engineer / Ingegnere per l' Ambiente e il territorio nr.
3587 (Albo A - Provincia di Trento)
cell: +39 3282818564
email: [hidden email],[hidden email],
[hidden email]
PEC: [hidden email]
URL: www.rendena100.eu
LinkedIn: https://www.linkedin.com/in/emanuele-cordano-31995333
GitHub: https://github.com/ecor

        [[alternative HTML version deleted]]

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

Re: stars::RasterIO using extent info?

Thu, 11/15/2018 - 06:31
Mike,
Thanks for the feedback - I'll check it out.
Best,
Tim

From: Michael Sumner [mailto:[hidden email]]
Sent: Wednesday, November 14, 2018 5:55 PM
To: Howard, Tim G (DEC) <[hidden email]>
Cc: [hidden email]
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?

FWIW I have a raster-delivering-front-end for RasterIO in this dev package: 

https://github.com/hypertidy/lazyraster

It uses the more obvious extent() idioms and will even use an open plot if nothing else is specified. 

(It uses an independent binding to GDAL in the vapour package). 

That might help, or not. It's on my list to find a sensible way for stars to leverage this obvious ease-of-use.  rgdal::readGDAL always had an interface to RasterIO, but only raster ever made use of that, and it's indirect - raster doesn't allow a mix of extent and resolution in its approach. 

Cheers, Mike. 

On Thu, 15 Nov 2018 at 02:27 Howard, Tim G (DEC) via R-sig-Geo <mailto:[hidden email]> wrote:
Ok, fair enough that there's no magic involved. I've worked through the details with the small example as follows. The result is only a couple of cells different in each direction.

library(stars)
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
x <- x[,,,1]
# calculate a circular polygon at the center of the raster
pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500)

bb_pol <- st_bbox(pol)

xoff <- st_dimensions(x)$x$offset
xdelt <- st_dimensions(x)$x$delta
yoff <- st_dimensions(x)$y$offset
ydelt <- st_dimensions(x)$y$delta
cropXoff <- (bb_pol$xmin - xoff + xdelt)/xdelt
cropXsize <- (bb_pol$xmax - bb_pol$xmin)/xdelt
cropYoff <- (bb_pol$ymin - yoff + ydelt)/ydelt
cropYsize <- (bb_pol$ymax - bb_pol$ymin)/ydelt

# if ydelt is negative, get abs of ysize and move yoffset to the top
if(cropYsize < 0){
  cropYsize <- abs(cropYsize)
  cropYoff <- cropYoff - cropYsize 
}

rasterio <- list(nXOff = cropXoff, nYOff = cropYoff, nXSize = cropXsize, nYSize = cropYsize, bands = c(1))
(z = read_stars(tif, RasterIO = rasterio))
plot(z)

# crop it if desired
plot(z[pol])

## compare to proxy method

x <- read_stars(tif, proxy = TRUE)
x <- x[,,,1]
y <- st_as_stars(x[pol])
plot(y)

# only a couple of cells off!
z
y




------------------------------
Date: Tue, 13 Nov 2018 17:26:16 +0100
From: Edzer Pebesma <mailto:[hidden email]>
To: mailto:[hidden email]
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <mailto:[hidden email]>
Content-Type: text/plain; charset="utf-8"



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
> mailto:[hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: tel:+49%20251%208333081

_______________________________________________
R-sig-Geo mailing list
mailto:[hidden email]
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
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: stars::RasterIO using extent info?

Wed, 11/14/2018 - 16:55
FWIW I have a raster-delivering-front-end for RasterIO in this dev package:

https://github.com/hypertidy/lazyraster

It uses the more obvious extent() idioms and will even use an open plot if
nothing else is specified.

(It uses an independent binding to GDAL in the vapour package).

That might help, or not. It's on my list to find a sensible way for stars
to leverage this obvious ease-of-use.  rgdal::readGDAL always had an
interface to RasterIO, but only raster ever made use of that, and it's
indirect - raster doesn't allow a mix of extent and resolution in its
approach.

Cheers, Mike.

On Thu, 15 Nov 2018 at 02:27 Howard, Tim G (DEC) via R-sig-Geo <
[hidden email]> wrote:

> Ok, fair enough that there's no magic involved. I've worked through the
> details with the small example as follows. The result is only a couple of
> cells different in each direction.
>
> library(stars)
> tif <- system.file("tif/L7_ETMs.tif", package = "stars")
> x <- read_stars(tif)
> x <- x[,,,1]
> # calculate a circular polygon at the center of the raster
> pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500)
>
> bb_pol <- st_bbox(pol)
>
> xoff <- st_dimensions(x)$x$offset
> xdelt <- st_dimensions(x)$x$delta
> yoff <- st_dimensions(x)$y$offset
> ydelt <- st_dimensions(x)$y$delta
> cropXoff <- (bb_pol$xmin - xoff + xdelt)/xdelt
> cropXsize <- (bb_pol$xmax - bb_pol$xmin)/xdelt
> cropYoff <- (bb_pol$ymin - yoff + ydelt)/ydelt
> cropYsize <- (bb_pol$ymax - bb_pol$ymin)/ydelt
>
> # if ydelt is negative, get abs of ysize and move yoffset to the top
> if(cropYsize < 0){
>   cropYsize <- abs(cropYsize)
>   cropYoff <- cropYoff - cropYsize
> }
>
> rasterio <- list(nXOff = cropXoff, nYOff = cropYoff, nXSize = cropXsize,
> nYSize = cropYsize, bands = c(1))
> (z = read_stars(tif, RasterIO = rasterio))
> plot(z)
>
> # crop it if desired
> plot(z[pol])
>
> ## compare to proxy method
>
> x <- read_stars(tif, proxy = TRUE)
> x <- x[,,,1]
> y <- st_as_stars(x[pol])
> plot(y)
>
> # only a couple of cells off!
> z
> y
>
>
>
>
> ------------------------------
> Date: Tue, 13 Nov 2018 17:26:16 +0100
> From: Edzer Pebesma <[hidden email]>
> To: [hidden email]
> Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="utf-8"
>
>
>
> 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
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081 <+49%20251%208333081>
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> 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
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: stars::RasterIO using extent info?

Wed, 11/14/2018 - 09:26
Ok, fair enough that there's no magic involved. I've worked through the details with the small example as follows. The result is only a couple of cells different in each direction.

library(stars)
tif <- system.file("tif/L7_ETMs.tif", package = "stars")
x <- read_stars(tif)
x <- x[,,,1]
# calculate a circular polygon at the center of the raster
pol <- x %>% st_bbox() %>% st_as_sfc() %>% st_centroid() %>% st_buffer(500)

bb_pol <- st_bbox(pol)

xoff <- st_dimensions(x)$x$offset
xdelt <- st_dimensions(x)$x$delta
yoff <- st_dimensions(x)$y$offset
ydelt <- st_dimensions(x)$y$delta
cropXoff <- (bb_pol$xmin - xoff + xdelt)/xdelt
cropXsize <- (bb_pol$xmax - bb_pol$xmin)/xdelt
cropYoff <- (bb_pol$ymin - yoff + ydelt)/ydelt
cropYsize <- (bb_pol$ymax - bb_pol$ymin)/ydelt

# if ydelt is negative, get abs of ysize and move yoffset to the top
if(cropYsize < 0){
  cropYsize <- abs(cropYsize)
  cropYoff <- cropYoff - cropYsize  
}

rasterio <- list(nXOff = cropXoff, nYOff = cropYoff, nXSize = cropXsize, nYSize = cropYsize, bands = c(1))
(z = read_stars(tif, RasterIO = rasterio))
plot(z)

# crop it if desired
plot(z[pol])

## compare to proxy method

x <- read_stars(tif, proxy = TRUE)
x <- x[,,,1]
y <- st_as_stars(x[pol])
plot(y)

# only a couple of cells off!
z
y




------------------------------
Date: Tue, 13 Nov 2018 17:26:16 +0100
From: Edzer Pebesma <[hidden email]>
To: [hidden email]
Subject: Re: [R-sig-Geo] stars::RasterIO using extent info?
Message-ID: <[hidden email]>
Content-Type: text/plain; charset="utf-8"



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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

st_transform problem

Wed, 11/14/2018 - 04:33

I have a shapefile from CBS/Kadaster in EPSG:28992 (Amersfoort / RD New) projection and I use st_transform to reproject to EPSG:4326. 
Small R script example:


setwd("c:/testTransform") library("sf") polyPurmerend <- st_read(dsn = "c:/testTransform", layer = "Purmerend",stringsAsFactors = FALSE) polyPurmerend_WGS84 <- st_transform(polyPurmerend,4326) st_write(polyPurmerend_WGS84,dsn=paste0("","Purmerend_WGS84",".gpkg"),"polygon", layer_options="OVERWRITE=YES")
Unfortunately the resulting layer is at a different position.


I cannot figure out why.
Rinus Deurloo University of Amsterdam.





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

Re: citing the spatio-temporal kriging vignette

Wed, 11/14/2018 - 02:54
Awesome!

Thanks so much!
Erin

Erin Hodgess, PhD
mailto: [hidden email]


On Wed, Nov 14, 2018 at 1:13 AM Roger Bivand <[hidden email]> wrote:

> On Wed, 14 Nov 2018, Erin Hodgess wrote:
>
> > Hello!
> >
> > What is the correct way to cite the (very fine) vignette on
> spatio-temporal
> > kriging, please?
>
> Unless it has changed a lot since publication, print(citation("gstat"),
> bibtex=TRUE) gives among others:
>
> @Article{,
>      title = {Spatio-Temporal Interpolation using gstat},
>      author = {Benedikt Gräler and Edzer Pebesma and Gerard Heuvelink},
>      year = {2016},
>      journal = {The R Journal},
>      volume = {8},
>      issue = {1},
>      pages = {204-218},
>      url = {
> https://journal.r-project.org/archive/2016-1/na-pebesma-heuvelink.pdf}
> }
>
> but should currently be (different URL):
>
> @article{RJ-2016-014,
>    author = {Benedikt Gräler and Edzer Pebesma and Gerard Heuvelink},
>    title = {{Spatio-Temporal Interpolation using gstat}},
>    year = {2016},
>    journal = {{The R Journal}},
>    url = {
> https://journal.r-project.org/archive/2016/RJ-2016-014/index.html},
>    pages = {204--218},
>    volume = {8},
>    number = {1}
> }
>
> Roger
>
> >
> > Thanks,
> > Erin
> >
> > Erin Hodgess, PhD
> > mailto: [hidden email]
> >
> >       [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > [hidden email]
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: [hidden email]
> http://orcid.org/0000-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
        [[alternative HTML version deleted]]

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

Re: citing the spatio-temporal kriging vignette

Wed, 11/14/2018 - 02:13
On Wed, 14 Nov 2018, Erin Hodgess wrote:

> Hello!
>
> What is the correct way to cite the (very fine) vignette on spatio-temporal
> kriging, please?

Unless it has changed a lot since publication, print(citation("gstat"),
bibtex=TRUE) gives among others:

@Article{,
     title = {Spatio-Temporal Interpolation using gstat},
     author = {Benedikt Gräler and Edzer Pebesma and Gerard Heuvelink},
     year = {2016},
     journal = {The R Journal},
     volume = {8},
     issue = {1},
     pages = {204-218},
     url = {https://journal.r-project.org/archive/2016-1/na-pebesma-heuvelink.pdf}
}

but should currently be (different URL):

@article{RJ-2016-014,
   author = {Benedikt Gräler and Edzer Pebesma and Gerard Heuvelink},
   title = {{Spatio-Temporal Interpolation using gstat}},
   year = {2016},
   journal = {{The R Journal}},
   url = {https://journal.r-project.org/archive/2016/RJ-2016-014/index.html},
   pages = {204--218},
   volume = {8},
   number = {1}
}

Roger

>
> Thanks,
> Erin
>
> Erin Hodgess, PhD
> mailto: [hidden email]
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: [hidden email]
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway

citing the spatio-temporal kriging vignette

Wed, 11/14/2018 - 00:13
Hello!

What is the correct way to cite the (very fine) vignette on spatio-temporal
kriging, please?

Thanks,
Erin

Erin Hodgess, PhD
mailto: [hidden email]

        [[alternative HTML version deleted]]

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

Re: stars::RasterIO using extent info?

Tue, 11/13/2018 - 10:26


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
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> --
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081

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

pEpkey.asc (2K) Download Attachment

stars::RasterIO using extent info?

Tue, 11/13/2018 - 09:10
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?


Thanks in advance for any tips.
Tim

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

How to extract values in an sf Points object

Sun, 11/11/2018 - 22:14
I have a a **tibble** containing objects of class **sfc_POINT** such as

    # A tibble: 1 x 1
              geometry
    *          <POINT [°]>
    1 (-86.82857 32.78871)

where, of course, the two dbl represent longitude and latitude that I wish
to combine to another data frame for plotting as the centroids of
geographic units that can be (1) adjusted manually for position (it's a
more or less permanent basemap, so I only have to do this once) and (2)
optionally located on a thematic map, in my case, the US. (I apologize for
our current political situation.)

Although I can do this by subsetting

    str(foo[[1]][2])
    num 32.8

to get the latitude of a point (which, blessedly is only displayed, not
truncated to one decimal) and by making the second subindex [[1]] for
latitude, what I want is to do this at a pass for a tibble column of 50.

This is possible, of course, with some fancy footwork and, if I had to do
it 50 times, I would. Before that, however, I'd like to be sure I haven't
missed a recipe.

My plan is to use this with the **usmap** package to create a standardized,
projected base map on which to create thematic maps. The package returns a
*ggplot* friendly object of a converted shapefile of US states and
optionally counties and subsets of states. **And** it insets Alaska and
Hawaii and discards  the farflung territories. Since the boundaries (for
now) of the states are stable, once created, it can be saved as an Rda
object along with its centroid data and dyplr::inner_joint on the state
fips.

Again, the object is to decompose a POINT object -- (-86.82857 32.78871)
into a list c(-86.82857, 32.78871) or some other structure amenable to
tidyverse treatment.

Thanks in advance.

Richard Careaga
Unaffiliated (so much more a gracious term than permanently unemployable)
@technocrat
https://technocrat.rbind.io

        [[alternative HTML version deleted]]

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

Re: R code for inverse distance to a power

Thu, 11/08/2018 - 17:42
On 11/8/18 05:50, sudheesh pai wrote:
> Hi,
>
> I am looking for an R code to replace the inverse distance to a power
> algorithm in QGIS, which has the input as point shape file. the parameters
> of algorithm are radius 1 (R1), radius 2 (R2) and the z parameter.
> Also, my CRS is 3006 which means that i have to specify R1 and R2 in
> meters.
>
> If some one has an example, kindly let me know.
>
> Thank you
>
Sounds like you are talking about Inverse Distance Weighted functions
(IDW). Yes the gstat package has a function for that.

Example:
https://mgimond.github.io/Spatial/interpolation-in-r.html

Enjoy,
Alex

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

R code for inverse distance to a power

Thu, 11/08/2018 - 07:50
Hi,

I am looking for an R code to replace the inverse distance to a power
algorithm in QGIS, which has the input as point shape file. the parameters
of algorithm are radius 1 (R1), radius 2 (R2) and the z parameter.
Also, my CRS is 3006 which means that i have to specify R1 and R2 in
meters.

If some one has an example, kindly let me know.

Thank you

--
SUDHEESH.P.G
[hidden email]

        [[alternative HTML version deleted]]

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

Pages