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: 1 hour 59 min ago

Re: Problem with initGRASS

Fri, 11/23/2018 - 02:44
On Fri, 23 Nov 2018, Jaime Burbano Girón wrote:

> Hi everyone,
>
> I'm trying to use GRASS in R, but I can't get to initialize it. I'm using
> Ubuntu 18.04 and GRASS 7.4.2. This is the system information:
>
> Versión de GRASS: 7.4.2
>
> Revisión SVN de GRASS: exported
>
> Fecha de compilación: 2018-10-28
>
> Construir plataforma: x86_64-pc-linux-gnu
>
> GDAL: 2.2.3
>
> PROJ.4: 4.9.3
>
> GEOS: 3.6.2
>
> SQLite: 3.22.0
>
> Python: 2.7.15rc1
>
> wxPython: 3.0.2.0
>
> Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic
>
>
> I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
> And this is the error: No, ?initGRASS says:

gisBase: The directory path to GRASS binaries and libraries

and grass74/bin only contains (for me):

$ ls bin
grass74

which is a Python script, ASCII text executable. Trying

grass74/grass-7.4.2

(no idea where your installer puts this):

$ ls grass-7.4.2
AUTHORS           contributors_extra.csv  fonts    REQUIREMENTS.html
bin               COPYING                 GPL.TXT  scripts
CHANGES           demolocation            gui      share
CITING            docs                    include  tools
config.status     driver                  INSTALL  translators.csv
contributors.csv  etc                     lib

which contains the binaries in bin and the libraries in lib. The required
string is what is set in the shell variable GISBASE when running GRASS
interactively:

echo $GISBASE

hence the name of the argument. initGRASS() cannot discover this as GRASS
is not running, and has not set its environment variables. I have added a
test looking for a directory called bin in the directory given as the
gisBase argument, check on https://r-forge.r-project.org/R/?group_id=2020 
to see when rgrass7 0.1-13 (Rev: 68) has build status current, try that
with your current setting, the correct setting and report back.

Hope this clarifies,

Roger


>
> sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
> not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
> g.version: not foundError in system(paste("g.version", get("addEXE",
> envir = .GRASS_CACHE),  :
>  error in running commandAdemás: Warning messages:1: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command2: In system(paste(paste("g.gisenv",
> get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command4: In system(paste(paste("g.gisenv",
> get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
> system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
> error in running command
>
>
> I verified that g.gisenv and g.version are located in /usr/lib/grass74/bin.
> Any suggestion?
>
> Thanks in advance.
>
> Best,
>
> Jaime.
> --
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

Problem with initGRASS

Thu, 11/22/2018 - 21:42
Hi everyone,

I'm trying to use GRASS in R, but I can't get to initialize it. I'm using
Ubuntu 18.04 and GRASS 7.4.2. This is the system information:

Versión de GRASS: 7.4.2

Revisión SVN de GRASS: exported

Fecha de compilación: 2018-10-28

Construir plataforma: x86_64-pc-linux-gnu

GDAL: 2.2.3

PROJ.4: 4.9.3

GEOS: 3.6.2

SQLite: 3.22.0

Python: 2.7.15rc1

wxPython: 3.0.2.0

Plataforma: Linux-4.15.0-38-generic-x86_64-with-Ubuntu-18.04-bionic


I'm running: > initGRASS(gisBase="/usr/lib/grass74/bin",override=T)
And this is the error:

sh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv:
not foundsh: 1: g.gisenv: not foundsh: 1: g.gisenv: not foundsh: 1:
g.version: not foundError in system(paste("g.version", get("addEXE",
envir = .GRASS_CACHE),  :
  error in running commandAdemás: Warning messages:1: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command2: In system(paste(paste("g.gisenv",
get("addEXE", envir = .GRASS_CACHE),  :  error in running command3: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command4: In system(paste(paste("g.gisenv",
get("addEXE", envir = .GRASS_CACHE),  :  error in running command5: In
system(paste(paste("g.gisenv", get("addEXE", envir = .GRASS_CACHE),  :
 error in running command


I verified that g.gisenv and g.version are located in /usr/lib/grass74/bin.
Any suggestion?

Thanks in advance.

Best,

Jaime.
--
Jaime Burbano Girón
Estudiante Doctoral
Doctorado en Estudios Ambientales y Rurales
Pontificia Universidad Javeriana
Bogotá, Colombia
(57-1) 3208320 Ext. 4844

        [[alternative HTML version deleted]]

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

Re: Doubt: add a boderline map

Wed, 11/21/2018 - 10:42
Yes, it worked.

Thank you very much!

Lara

Hugo Costa <[hidden email]> escreveu no dia terça, 20/11/2018 à(s)
18:15:

> Hi Lara,
>
> have you tried to run the following?
>
> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
> plot(nc.sids, asp=1, main="", add=TRUE)
>
> Hugo
>
> Lara Silva <[hidden email]> escreveu no dia terça, 20/11/2018
> à(s) 16:04:
>
>> Hello,
>>
>> I got a map through the plot function.
>>
>> Code:
>>
>> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
>> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>>
>> However I needed to add a border line in my study area to become more
>> perceptible (Figure 1).
>>
>> I tried the border(), rect () functions  but it gave error.
>>
>> So, I tried to create a border line
>>
>> Code:
>> nc.sids <- readShapePoly("PoligonoSM.shp")
>>
>> plot(nc.sids, asp=1, main="")
>>
>> nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
>> plot(nc.border, asp=1, main="")
>>
>>
>> But I am having a lot of difficulties in building the script.
>>
>> How can cross the information of the 2 plots?
>>
>> https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
>> Thanks.
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> 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

job position on geostatistics at the London School of Hygiene & Tropical Medicine

Wed, 11/21/2018 - 06:57
Dear community,

We are looking for a geostatistician interested in spatial modelling of infectious diseases to join our group (London Applied and Spatial Epidemiology Research group) at the London School of Hygiene & Tropical Medicine. More details in the link below. Also, further information in our research projects in these links: www.thiswormyworld.org<http://www.thiswormyworld.org/> / http://laser.lshtm.ac.uk/

__________________________________________________________________________

Assistant Professor in Spatial Statistics and Epidemiology
Dept of Disease Control, LSHTM
Salary: �45,878 to �52,520 per annum, inclusive.
Closing Date:  Wednesday 05 December 2018

The successful applicant will be conducting research on the epidemiology and control of neglected tropical diseases (NTDs) in sub-Saharan Africa. This research is a component of the Expanded Special Project for the Elimination of NTDs within the Africa Region of the WHO. The post-holder will be responsible for leading the development and implementation of appropriate statistical methods to quantify the spatial distribution of NTDs.
Job Description and application: https://jobs.lshtm.ac.uk/vacancy.aspx?ref=ITD-DCD-2018-24


Thanks,

________________________________________________________
Jorge Cano | Assistant Professor | Department of Disease Control
London Applied & Spatial Epidemiology Research Group (LASER)
www.thiswormyworld.org<http://www.thiswormyworld.org/> / http://laser.lshtm.ac.uk/
T: +44 (0) 20 7927 2584



        [[alternative HTML version deleted]]


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

Re: Doubt: add a boderline map

Tue, 11/20/2018 - 13:15
Hi Lara,

have you tried to run the following?

map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)

plot(nc.sids, asp=1, main="", add=TRUE)

Hugo

Lara Silva <[hidden email]> escreveu no dia terça, 20/11/2018
à(s) 16:04:

> Hello,
>
> I got a map through the plot function.
>
> Code:
>
> map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
> main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)
>
> However I needed to add a border line in my study area to become more
> perceptible (Figure 1).
>
> I tried the border(), rect () functions  but it gave error.
>
> So, I tried to create a border line
>
> Code:
> nc.sids <- readShapePoly("PoligonoSM.shp")
>
> plot(nc.sids, asp=1, main="")
>
> nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
> plot(nc.border, asp=1, main="")
>
>
> But I am having a lot of difficulties in building the script.
>
> How can cross the information of the 2 plots?
>
> https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
> Thanks.
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> 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

Re: Unable to use factors in raster::predict

Tue, 11/20/2018 - 13:11
This is not of much help, but I run the code, and the line
pc <- predict(logo, m, OOB=TRUE, factors=f)
did work fine. I got a Rasterlayer, which is the R logo with pixels ranging
from 0 to 1.
Hugo

Gonzalez-Mirelis, Genoveva via R-sig-Geo <[hidden email]> escreveu
no dia terça, 20/11/2018 à(s) 10:23:

> 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
>
        [[alternative HTML version deleted]]

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

Doubt: add a boderline map

Tue, 11/20/2018 - 09:04
Hello,

I got a map through the plot function.

Code:

map <-plot(proj$A_EMmeanByROC_mergedAlgo_mergedRun_mergedData,
main= "", xlab ="x.coords", ylab="y.coords", cex.axis=caxis)

However I needed to add a border line in my study area to become more
perceptible (Figure 1).

I tried the border(), rect () functions  but it gave error.

So, I tried to create a border line

Code:
nc.sids <- readShapePoly("PoligonoSM.shp")

plot(nc.sids, asp=1, main="")

nc.border <- unionSpatialPolygons(nc.sids, rep(1, nrow(nc.sids)))
plot(nc.border, asp=1, main="")


But I am having a lot of difficulties in building the script.

How can cross the information of the 2 plots?

https://www.dropbox.com/sh/y5by1sxjcvdl6rh/AAAIF4srkChW2FxeEfFG4d_ta?dl=0
Thanks.

        [[alternative HTML version deleted]]

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

Unable to use factors in raster::predict

Tue, 11/20/2018 - 03:22
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

Pages