Re: [R-sig-Geo] rgdal release candidate 1.5-9 rev. 1000 ready for testing

[Thiago V. dos Santos via R-sig-Geo]

Hi Roy,

I use and recommend MacPorts as a package manager for macOS computers. It is 
quite comprehensive in terms of programs/libraries available, almost always 
up-to-date and pretty simple to use. It does not always provide binaries 
though, it needs to compile a few packages from source (gdal is an example) - 
which can take a little bit, depending on your machine configuration.

I have some notes on how to set it up - maybe it's useful for you: 
https://thiagodossantos.com/post/1-mac-science-software/  

Greetings,
 -- Thiago V. dos Santos

ThiagoDosSantos.com
MudancasClimaticasBrasil.com






On Friday, June 5, 2020, 10:54:55 AM GMT-3, Roy Mendelssohn - NOAA Federal via 
R-sig-Geo  wrote: 





Hi All:

Is there a source for binaries of PROJ >= 6 and GDAL > 3  for the Mac, other 
than the Anaconda ones which cause problems because they use @rpart?

Thanks,

-Roy


> On Jun 5, 2020, at 4:29 AM, Roger Bivand  wrote:
> 
> The release candidate of rgdal_1.5-9 is ready for testing on R-forge:
> 
> https://r-forge.r-project.org/R/?group_id=884
> 
> Those insisting on installing on PROJ >= 6 and GDAL < 3 must use configure 
> argument --with-proj_api="proj_api.h"; with this used, this version builds 
> with --no-build-vignettes and installs and checks OK on PROJ 7.0.1 and GDAL 
> 2.2.4 with --with-proj_api="proj_api.h".
> 
> Otherwise checked OK with PROJ 4.8.0, 4.9.2, 4.9.3 and 5.2.0 with GDAL 
> 1.11.4; with PROJ 5.2.0 and GDAL 2.2.4, 2.3.2 and 2.4.2; with PROJ 6.3.2 and 
> GDAL 3.0.4; with PROJ 7.0.1 and GDAL 3.0.4 and 3.1.0.
> 
> All who have indicated issues with source installs are asked to try the 
> release candidate and to report back here by midnight CEST Monday 8 June. If 
> no indications are forthcoming, I'll assume that problems with 1.5-8 are 
> resolved.
> 
> Roger
> 
> -- 
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; e-mail: roger.biv...@nhh.no
> https://orcid.org/-0003-2392-6140
> https://scholar.google.no/citations?user=AWeghB0J=en
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

**
"The contents of this message do not reflect any position of the U.S. 
Government or NOAA."
**
Roy Mendelssohn
Supervisory Operations Research Analyst
NOAA/NMFS
Environmental Research Division
Southwest Fisheries Science Center
***Note new street address***
110 McAllister Way
Santa Cruz, CA 95060
Phone: (831)-420-3666
Fax: (831) 420-3980
e-mail: roy.mendelss...@noaa.gov www: https://www.pfeg.noaa.gov/

"Old age and treachery will overcome youth and skill."
"From those who have been given much, much will be expected" 
"the arc of the moral universe is long, but it bends toward justice" -MLK Jr.


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] rgdal 1.5-8 released on CRAN

[Thiago V. dos Santos via R-sig-Geo]

Hi Roger, 

Thanks for all the hints. However, I am still struggling with this.

I have PKG_CONFIG_PATH set to /opt/local/lib/proj6/lib/pkgconfig:
cat $PKG_CONFIG_PATH/proj.pc
prefix=/opt/local/lib/proj6
exec_prefix=${prefix}
libdir=${exec_prefix}/lib
includedir=${prefix}/include
datadir=${prefix}/share/proj

Name: proj
Description: Cartographic Projections Library.
Requires:
Version: 6.3.2
Libs: -L${libdir} -lproj
Libs.Private: -L/opt/local/lib -lsqlite3 -lstdc++
Cflags: -I${includedir}



Note that the paths described within proj.pc are all good:
ls /opt/local/lib/proj6/lib/
libproj.15.dylib* libproj.a         libproj.dylib@    pkgconfig/

ls /opt/local/lib/proj6/include/
geodesic.h            proj/                 proj_api.h            
proj_experimental.h
org_proj4_PJ.h        proj.h                proj_constants.h      
proj_symbol_rename.h

ls /opt/local/lib/proj6/share/proj/
BETA2007.gsb          ITRF2008              WO                    nad.lst       
        null                  prvi
CH                    ITRF2014              alaska                nad27         
        nzgd2kgrid0005.gsb    stgeorge
FL                    MD                    conus                 nad83         
        other.extra           stlrnc
GL27                  TN                    egm96_15.gtx          ntf_r93.gsb   
        proj.db               stpaul
ITRF2000              WI                    hawaii                ntv1_can.dat  
        projjson.schema.json  world



It is also recognized by the system:
pkg-config --list-all
...
netcdf                              netCDF - NetCDF Client Library for C
cairo-xml                           cairo-xml - xml surface backend for cairo 
graphics library
proj                                proj - Cartographic Projections Library.
damageproto                         DamageProto - Damage extension headers
libraw_r                            libraw - Raw image decoder library 
(thread-safe)
...

pkg-config proj --libs
-L/opt/local/lib/proj6/lib -lproj



I am sure that my system's gdal was built against proj6: 
gdal-config --dep-libs
-L/opt/local/lib -lIlmImf -lImath -lHalf -lIex -lIexMath -lIlmThread -lcrypto 
-lqhull -lqhull -L/opt/local/lib -lgeos_c -lwebp -L/opt/local/lib -lsqlite3 
-L/opt/local/lib -lexpat -L/opt/local/lib -lopenjp2 -L/opt/local/include 
-L/opt/local/include/lib -ljasper -L/opt/local/lib -lnetcdf -L/opt/local/lib 
-lhdf5 -L/opt/local/lib -lmfhdf -ldf -lsz -lgif -lCharLS -lpng -lzstd 
-L/opt/local/lib/proj6/lib -lproj -lsqlite3 -lz -lpthread -ldl -L/usr/local/lib 
-L/opt/local/lib -lspatialite -lpcre -L/opt/local/lib -lcurl -L/opt/local/lib 
-lxml2



and yet, this is what I get when I try to install rgdal relying only on 
pkg-config:

R
install.packages('rgdal', type='source')

Installing package into ‘/Users/thiago/Documents/r-packages’
(as ‘lib’ is unspecified)
trying URL 'https://cran.rstudio.com/src/contrib/rgdal_1.5-8.tar.gz'
Content type 'application/x-gzip' length 2299235 bytes (2.2 MB)
==
downloaded 2.2 MB


* installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
** using staged installation
configure: R_HOME: /Library/Frameworks/R.framework/Resources
configure: CC: /opt/local/bin/clang
configure: CXX: /opt/local/bin/clang++
configure: C++11 support available
configure: rgdal: 1.5-8
checking for /usr/bin/svnversion... yes
configure: svn revision: 990
checking for gdal-config... /opt/local/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 3.1.0
checking GDAL version >= 1.11.4... yes
checking GDAL version <= 2.5 or >= 3.0... yes
checking gdal: linking with --libs only... yes
checking GDAL: gdal-config data directory readable... yes
checking GDAL: /opt/local/share/gdal/stateplane.csv readable... yes
configure: pkg-config proj exists, will use it
configure: PROJ version: 6.3.2
configure: proj CPP flags: -DPROJ_H_API -I/opt/local/lib/proj6/include
configure: PROJ LIBS: -L/opt/local/lib/proj6/lib -lproj
checking PROJ header API:... proj.h
checking for gcc... /opt/local/bin/clang
checking whether the C compiler works... yes
checking for C compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether /opt/local/bin/clang accepts -g... yes
checking for /opt/local/bin/clang option to accept ISO C89... none needed
checking how to run the C preprocessor... /opt/local/bin/clang -E
checking for grep that handles long lines and -e... /usr/bin/grep
checking for egrep... /usr/bin/grep -E
checking for ANSI C header files... rm: conftest.dSYM: is a directory
rm: conftest.dSYM: is a directory
yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes

Re: [R-sig-Geo] rgdal 1.5-8 released on CRAN

[Thiago V. dos Santos via R-sig-Geo]

Dear Roger,

Many thanks for the effort keeping rgdal up-to-date with proj6.

I'd like to report that I am unable to install rgdal 1.5.8 on my macOS system. 
I am reporting this error here on the list because I thought it would be the 
best channel in terms of reaching future users experiencing the same error. 
Please apologize if my rationale is not right, and ignore this message.

GDAL was installed via MacPorts and is at its latest version, 3.1.0 (at the 
time of this writing). Proj6 has also been installed via Macports.

This is how I am trying to install it:

export PKG_CONFIG_PATH=/opt/local/lib/proj6/lib/pkgconfig
R
install.packages('rgdal', type="source", configure.args=c(
     '--with-proj-include=/opt/local/lib/proj6/include',
     '--with-proj-lib=/opt/local/lib/proj6/lib'))

and this is the output I get:

R version 4.0.0 (2020-04-24) -- "Arbor Day"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)


R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.


  Natural language support but running in an English locale


R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.


Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.


> install.packages('rgdal', type="source", configure.args=c(
+      '--with-proj-include=/opt/local/lib/proj6/include',
+      '--with-proj-lib=/opt/local/lib/proj6/lib'))
Installing package into ‘/Users/thiago/Documents/R-packages’
(as ‘lib’ is unspecified)


trying URL 'https://cran.r-project.org/src/contrib/rgdal_1.5-8.tar.gz'
Content type 'application/x-gzip' length 2299235 bytes (2.2 MB)
==
downloaded 2.2 MB


* installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
** using staged installation
configure: R_HOME: /Library/Frameworks/R.framework/Resources
configure: CC: /opt/local/bin/gcc
configure: CXX: /opt/local/bin/g++
configure: C++11 support available
configure: rgdal: 1.5-8
checking for /usr/bin/svnversion... yes
configure: svn revision: 990
checking for gdal-config... /opt/local/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 3.1.0
checking GDAL version >= 1.11.4... yes
checking GDAL version <= 2.5 or >= 3.0... yes
checking gdal: linking with --libs only... yes
checking GDAL: gdal-config data directory readable... yes
checking GDAL: /opt/local/share/gdal/stateplane.csv readable... yes
configure: pkg-config proj exists, will use it
configure: PROJ version: 6.3.2
configure: proj CPP flags: -DPROJ_H_API -I/opt/local/lib/proj6/include
configure: PROJ LIBS: -L/opt/local/lib/proj6/lib
checking PROJ header API:... proj.h
checking for gcc... /opt/local/bin/gcc
checking whether the C compiler works... yes
checking for C compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C compiler... yes
checking whether /opt/local/bin/gcc accepts -g... yes
checking for /opt/local/bin/gcc option to accept ISO C89... none needed
checking how to run the C preprocessor... /opt/local/bin/gcc -E
checking for grep that handles long lines and -e... /usr/bin/grep
checking for egrep... /usr/bin/grep -E
checking for ANSI C header files... rm: conftest.dSYM: is a directory
rm: conftest.dSYM: is a directory
yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking proj.h usability... yes
checking proj.h presence... yes
checking for proj.h... yes
checking for proj_context_create in -lproj... no
configure: error: proj_context_create not found in libproj.
ERROR: configuration failed for package ‘rgdal’
* removing ‘/Users/thiago/Documents/R-packages/rgdal’
* restoring previous ‘/Users/thiago/Documents/R-packages/rgdal’


The downloaded source packages are in

‘/private/var/folders/_z/01gg71zs19g816v6m2dddt8wgn/T/RtmpvZAChj/downloaded_packages’
Warning message:
In install.packages("rgdal", type = "source", configure.args = 
c("--with-proj-include=/opt/local/lib/proj6/include",  :
  installation of package ‘rgdal’ had non-zero exit status

I found this error to be quite mysterious, and could not find any previous 
discussion about it.

Are you familiar with it? Is there any other argument that I can pass to 
install.packages to solve it?

Greetings,
 -- Thiago V. dos Santos

ThiagoDosSantos.com
MudancasClimaticasBrasil.com






On Thursday, May 

[R-sig-Geo] How to find distance between grid cell and point that lies in specified radius?

[Thiago V. dos Santos via R-sig-Geo]

Hello all,

I am trying to implement a function in R to perform the Cressman scheme (which 
corrects the values of a gridded field, say precipitation, based on the closest 
observations available).

It looks like it hasn't been done yet, but please let me know if you are aware 
of any R package that implements the Cressman scheme.

I have a raster and a spatial points layer (which I will call "stations"). A 
toy example is provided below:

```
library(raster)


# create random raster
r <- raster(nrows=9, ncols=18, xmn=5153337, xmx=6570069, ymn=7462732, 
            ymx=8060416, crs = CRS("+proj=poly +lat_0=0 +lon_0=-54 +x_0=500 
+y_0=1000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"))
r[] <- runif(ncell(r))


# create spatial points
pts <- sampleRandom(r, 10, na.rm = TRUE, sp = TRUE)


# display everything
plot(r)
plot(pts, add=T, col="red", pch=16)
```

The starting point of my analysis would be:

1) for each grid cell in the raster, identify the station(s) that lie *within a 
radius of influence* (say 30 km);

2) for each station lying under the radius of influence, calculate a 
distance-weighted expression that will be used later to correct the raster 
values. Its formula is:

W = (R2 - r2)/(R2 + r2) 

where R = influence radius and r = distance between the station and the grid 
cell.

3) Based on the weighted correction factor W calculated above, update the value 
of the grid cell using the expression W * (gridpoint - station). The actual 
expression is a bit more sophisticated, but I want to get the general method 
here.

Based on some research, it looks like this could be done with sf. But since I 
don't have any experience with that package, would someone more familiar please 
point out some functions in sf or even other packages that could be helpful to 
solve my problem?

Thanks in advance,
 -- Thiago V. dos Santos

Postdoctoral Research Fellow
Department of Climate and Space Science and Engineering
University of Michigan

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] GDAL 3.0.0 and rgdal?

[Thiago V. dos Santos via R-sig-Geo]

Dear all,
Yesterday I replaced GDAL v2.4.1 with the latest v3.0.0 on my Linux server. 
Then, when I tried to reinstall rgdal, I got an error due to wrong system 
requirements: rgdal requires GDAL between versions 1.11.4 and 2.5.0.
Therefore, I am wondering: are there plans for rgdal to support GDAL 3.0.0 in 
the next few days? I like to be up-to-date with GIS software, but I don't 
really need to use GDAL 3.0.0 at all. If it is going to be a little while until 
it is supported by rgdal, I can happily revert back to v2.4.1.
Thanks, -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] rgdal installation on Fedora 28

[Thiago V. dos Santos via R-sig-Geo]

The following lines in your config.log:

/usr/bin/ld: /tmp/cc9pfZ1b.o: relocation R_X86_64_32 against .rodata' can not 
be used when making a shared object; recompile with -fPIC

make me wonder: maybe there is a -fPIC compilation flag lacking for you 
platform in the makefile?

Greetings, -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan 

On Friday, June 15, 2018, 3:27:44 PM EDT, Sarah Goslee 
 wrote:  
 
 Hi folks,

I updated all of my linux computers to Fedora 28, and now can’t install rgdal.

I can install sf with no problems, so it doesn't appear to be an issue
with GDAL or with proj.

I checked with Roger Bivand, the package maintainer, who asked me to
confirm whether I could install sf, to make sure it wasn’t a
dependency issue, and to post the logs online, then ask on the
R-sig-geo mailing list. I’m putting the full problem here, and the
abbreviated version on the mailing list.

I’d appreciate any thoughts: I rely very heavily on the incredibly
useful rgdal package.

> sessionInfo()
R version 3.5.0 (2018-04-23)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: Fedora 28 (Workstation Edition)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_US.UTF-8      LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8      LC_NAME=C
 [9] LC_ADDRESS=C              LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.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.0 tools_3.5.0


The install message is:

installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
configure: CC: gcc -m64
configure: CXX: g++ -m64
configure: rgdal: 1.3-2
checking for /usr/bin/svnversion... no
configure: svn revision: 755
checking whether the C++ compiler works... yes
checking for C++ compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... configure: error: in
`/tmp/RtmpfMGBY5/R.INSTALL66c3b8deddb/rgdal':
configure: error: cannot run C++ compiled programs.
If you meant to cross compile, use `--host'.
See `config.log' for more details
ERROR: configuration failed for package ‘rgdal’
* removing ‘/usr/lib64/R/library/rgdal’
* restoring previous ‘/usr/lib64/R/library/rgdal’


The complete information and config.log are posted online at:

http://numberwright.com/2018/06/stuck-on-rgdal/


Any ideas? Or any other information that would be useful?

Thanks!
Sarah

-- 
Sarah Goslee
http://www.functionaldiversity.org

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo  
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] CRAN releases of sp, rgdal and rgeos

[Thiago V. dos Santos via R-sig-Geo]

Hi Roger, thanks for the follow up. Answering to your comments:

a) I am using Macports's gdal version 2.3.0.20180523. I am pretty sure that its 
source code (i.e. the one that is compiled on my machine) is downloaded from 
OSGEO (at least according to 
https://github.com/macports/macports-ports/blob/master/gis/gdal/Portfile);

b) I do have more than one libproj installed on my system (as QGIS for macOS 
requires its own gdal and proj stuff), but only the Macports one is on my path;


c) I tried running:
tools::testInstalledPackage("rgdal", 
outDir=tempdir())list.files(tempdir())file.show(file.path(tempdir(), 
"rgdal-Ex.Rout"))

but it unexpectedly returned "Error in library("rgdal") : there is no package 
called ‘rgdal’. Execution halted".

It looks weird but so far I am able to normally load and use at least a couple 
of rgdal's functions.

I will keep monitoring it though, and report any other issue.

Cheers,
 -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan 

On Sunday, June 10, 2018, 4:46:04 AM EDT, Roger Bivand 
 wrote:  
 
 Thanks for reporting, comments inline below (note that I have no OSX 
access at all):

On Sat, 9 Jun 2018, Thiago V. dos Santos wrote:

> Dear Roger,
> Thank you very much for the excellent work done with those packages.
> Today I update both rgeos and rgdal on my system - macOS 10.13.5 with 
> all dependencies installed via MacPorts. I had updated sp a few days 
> earlier. The dependencies versions on my system are gdal 
> @2.3.0.20180523_0+grib+hdf4+hdf5+jasper+mpich+netcdf and proj @5.1.0_0. 
> Rgeos's update went flawlessly, but I got a few errors while updating 
> rgdal. It still compiled successfully, but I am concerned that some 
> functionality might be compromised due to the errors.
> This is what I got (I have to manually specify the location of proj or 
> rgdal won't find it):> install.packages('rgdal', type = "source", 
> configure.args=c(
> +    '--with-proj-include=/opt/local/lib/proj5/include',
> +    '--with-proj-lib=/opt/local/lib/proj5/lib'))
> Installing package into ‘/Users/thiago/Documents/R-packages’
> (as ‘lib’ is unspecified)
> trying URL 'https://cran.r-project.org/src/contrib/rgdal_1.3-2.tar.gz'
> Content type 'application/x-gzip' length 1667049 bytes (1.6 MB)
> ==
> downloaded 1.6 MB
>
> * installing *source* package ‘rgdal’ ...
> ** package ‘rgdal’ successfully unpacked and MD5 sums checked
> configure: CC: /usr/bin/clang
> configure: CXX: /usr/bin/clang++
> configure: rgdal: 1.3-2
> checking for /usr/bin/svnversion... yes
> configure: svn revision: 755
> checking whether the C++ compiler works... yes
> checking for C++ compiler default output file name... a.out
> checking for suffix of executables...
> checking whether we are cross compiling... no
> checking for suffix of object files... o
> checking whether we are using the GNU C++ compiler... yes
> checking whether /usr/bin/clang++ accepts -g... yes
> checking whether /usr/bin/clang++ supports C++11 features by default... no
> checking whether /usr/bin/clang++ supports C++11 features with 
> -std=gnu++11... yes
> configure: C++11 support available
> checking for gdal-config... /opt/local/bin/gdal-config
> checking gdal-config usability... yes
> configure: GDAL: 2.4.0

All OK up to the GDAL version returned by gdal-config - are you using the 
released GDAL 2.3.0 (probably not) or master?

> checking C++11 support for GDAL >= 2.3.0... yes
> checking GDAL version >= 1.11.4... yes
> checking gdal: linking with --libs only... yes
> checking GDAL: /opt/local/share/gdal/pcs.csv readable... yes
> checking proj_api.h presence and usability... yes
> ./configure: line 3395: test: =: unary operator expected

Will check, that line is:

if test ${PROJ_VERSION} = "" ; then

from configure.ac line 305. Possibly a shell dialect issue.

> checking PROJ version >= 4.8.0... yes
> checking projects.h presence and usability... yes

These relate to configure.ac lines 376-419, and the outcome: epsg found 
and readable is OK - could there be two libproj on your system (maybe for 
different architectures)?

> Undefined symbols for architecture x86_64:
>  "_pj_ctx_fclose", referenced from:
>      _main in proj_conf_test2-06fe7d.o
>  "_pj_get_default_ctx", referenced from:
>      _main in proj_conf_test2-06fe7d.o
>  "_pj_open_lib", referenced from:
>      _main in proj_conf_test2-06fe7d.o
> ld: symbol(s) not found for architecture x86_64
> clang: error: linker command failed with exit code 1 (use -v to see 
> invocation)
> ./configure: line 3511: ./proj_conf_test2: No such file or directory
> checking PROJ.4: epsg found and readable... yes

Same here for next block in configure.ac; conus found and readable.

> Undefined symbols for architecture x86_64:
>  "_pj_ctx_fclose", referenced from:
>      _main in proj_conf_test3-3b7aa2.o
>  "_pj_get_default_ctx", referenced from:
>      _main in 

Re: [R-sig-Geo] CRAN releases of sp, rgdal and rgeos

[Thiago V. dos Santos via R-sig-Geo]

Dear Roger,
Thank you very much for the excellent work done with those packages.
Today I update both rgeos and rgdal on my system - macOS 10.13.5 with all 
dependencies installed via MacPorts. I had updated sp a few days earlier.
The dependencies versions on my system are gdal 
@2.3.0.20180523_0+grib+hdf4+hdf5+jasper+mpich+netcdf and proj @5.1.0_0.
Rgeos's update went flawlessly, but I got a few errors while updating rgdal. It 
still compiled successfully, but I am concerned that some functionality might 
be compromised due to the errors.
This is what I got (I have to manually specify the location of proj or rgdal 
won't find it):> install.packages('rgdal', type = "source", configure.args=c(
+ '--with-proj-include=/opt/local/lib/proj5/include',
+ '--with-proj-lib=/opt/local/lib/proj5/lib'))
Installing package into ‘/Users/thiago/Documents/R-packages’
(as ‘lib’ is unspecified)
trying URL 'https://cran.r-project.org/src/contrib/rgdal_1.3-2.tar.gz'
Content type 'application/x-gzip' length 1667049 bytes (1.6 MB)
==
downloaded 1.6 MB

* installing *source* package ‘rgdal’ ...
** package ‘rgdal’ successfully unpacked and MD5 sums checked
configure: CC: /usr/bin/clang
configure: CXX: /usr/bin/clang++
configure: rgdal: 1.3-2
checking for /usr/bin/svnversion... yes
configure: svn revision: 755
checking whether the C++ compiler works... yes
checking for C++ compiler default output file name... a.out
checking for suffix of executables... 
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C++ compiler... yes
checking whether /usr/bin/clang++ accepts -g... yes
checking whether /usr/bin/clang++ supports C++11 features by default... no
checking whether /usr/bin/clang++ supports C++11 features with -std=gnu++11... 
yes
configure: C++11 support available
checking for gdal-config... /opt/local/bin/gdal-config
checking gdal-config usability... yes
configure: GDAL: 2.4.0
checking C++11 support for GDAL >= 2.3.0... yes
checking GDAL version >= 1.11.4... yes
checking gdal: linking with --libs only... yes
checking GDAL: /opt/local/share/gdal/pcs.csv readable... yes
checking proj_api.h presence and usability... yes
./configure: line 3395: test: =: unary operator expected
checking PROJ version >= 4.8.0... yes
checking projects.h presence and usability... yes
Undefined symbols for architecture x86_64:
  "_pj_ctx_fclose", referenced from:
  _main in proj_conf_test2-06fe7d.o
  "_pj_get_default_ctx", referenced from:
  _main in proj_conf_test2-06fe7d.o
  "_pj_open_lib", referenced from:
  _main in proj_conf_test2-06fe7d.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
./configure: line 3511: ./proj_conf_test2: No such file or directory
checking PROJ.4: epsg found and readable... yes
Undefined symbols for architecture x86_64:
  "_pj_ctx_fclose", referenced from:
  _main in proj_conf_test3-3b7aa2.o
  "_pj_get_default_ctx", referenced from:
  _main in proj_conf_test3-3b7aa2.o
  "_pj_open_lib", referenced from:
  _main in proj_conf_test3-3b7aa2.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)
./configure: line 3570: ./proj_conf_test3: No such file or directory
checking PROJ.4: conus found and readable... yes
configure: Package CPP flags:  -I/opt/local/include 
-I/opt/local/lib/proj5/include
configure: Package LIBS:  -L/opt/local/lib -lgdal -L/opt/local/lib/proj5/lib 
-lproj
configure: creating ./config.status
config.status: creating src/Makevars
** libs
clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" 
-DNDEBUG -I/opt/local/include -I/opt/local/lib/proj5/include 
-I"/Users/thiago/Documents/R-packages/sp/include" -I/usr/local/include   -fPIC  
-Wall -g -O2 -c OGR_write.cpp -o OGR_write.o
clang++ -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" 
-DNDEBUG -I/opt/local/include -I/opt/local/lib/proj5/include 
-I"/Users/thiago/Documents/R-packages/sp/include" -I/usr/local/include   -fPIC  
-Wall -g -O2 -c gdal-bindings.cpp -o gdal-bindings.o
/usr/bin/clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG 
-I/opt/local/include -I/opt/local/lib/proj5/include 
-I"/Users/thiago/Documents/R-packages/sp/include" -I/usr/local/include   -fPIC  
-Wall -g -O2  -c init.c -o init.o
/usr/bin/clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG 
-I/opt/local/include -I/opt/local/lib/proj5/include 
-I"/Users/thiago/Documents/R-packages/sp/include" -I/usr/local/include   -fPIC  
-Wall -g -O2  -c inverser.c -o inverser.o
/usr/bin/clang -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG 
-I/opt/local/include -I/opt/local/lib/proj5/include 
-I"/Users/thiago/Documents/R-packages/sp/include" -I/usr/local/include   -fPIC  
-Wall -g -O2  -c local_stubs.c -o local_stubs.o

[R-sig-Geo] Substitute raster stack values based on matrix index

[Thiago V. dos Santos via R-sig-Geo]

I have a large (7000x7000, 10-layered) raster stack whose values range from 0 
to 100.

I need to modify the raster values using the a "lookup table" consisted of a 
matrix which is 100 rows long by 10 cols wide, where the number of rows is 
associated with the 0-100 value range of the raster and the number of columns 
relates to the number of layers of the raster stack. 

The following code creates a toy example of my dataset:

# Create raster stackr <- raster(ncol=50, nrow=50)s <- stack(lapply(1:10, 
function(x) setValues(r, 1:ncell(r

# Create lookup tablem <- matrix(sample(100, 100*10, TRUE), 100, 10)

Therefore, I need to loop through each cell of the raster and check its value. 
For example, if the cell value is `20`, then the new value will be reassigned 
from the 20th line and 1st column in the matrix.

And so on for the 2nd raster stack layer and 2nd matrix column. And so on for 
the remaining raster stack layers and matrix columns.

Any ideas on how to do that?

I know this sounds a bit cumbersome, but that's how a few ISRIC-WISE's soil 
datasets are organized!
Greetings, -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Is the raster package still receiving updates?

[Thiago V. dos Santos via R-sig-Geo]

t; > > > Fra: Michael Sumner > Sendt: lørdag 21. oktober,
> 06.28 > Emne: Re: [R-sig-Geo] Is the raster package still receiving
> updates? > Til: Hodgess, Erin > Kopi: R-sig-geo Mailing List > > >
> Great! I'm keen, but haven't done the work needed to flesh this out
> properly: https://github.com/mdsumner/raster-rforge/issues I have not
> checked if there's been any commits to r-forge since I cloned this.
> Cheers, Mike On Sat, 21 Oct 2017, 12:58 Hodgess, Erin, wrote: > I could
> take a swat at it, if you wish. Please do send me the list of > fixes. I
> could start in about a week, if that would work. > Sincerely, > Erin > >
>> Erin M. Hodgess > Associate Professor > Department of Mathematics and
> Statistics > University of Houston - Downtown > mailto: hodge...@uhd.edu
> <mailto:hodge...@uhd.edu> > >  >
> From: R-sig-Geo on behalf of Michael > Sumner > Sent: Friday, October
> 20, 2017 8:05 PM > To: Thiago V. dos Santos > Cc: R-sig-geo Mailing List
>> Subject: Re: [R-sig-Geo] Is the raster package still receiving
> updates? > > No, it's effectively abandoned. No concrete plans known. >
>> I'm interested to keep it going and have kept a list of fixes needed,
> but > not sure when or if I'll get to it. > > I'd support any efforts in
> this, I'll rely on raster for many years still. > > Cheers, Mike > > On
> Sat, 21 Oct 2017, 02:05 Thiago V. dos Santos via R-sig-Geo, < >
> r-sig-geo@r-project.org <mailto:r-sig-geo@r-project.org>> wrote: > > >
> Dear list, > > > > I realized that the latest version of the raster
> package was released on > > CRAN over a year ago. > > > > I also noticed
> that Robert's participation in this list has become rather > > scarce. >
>> > > I was just wondering whether the raster package is still receiving
>> updates? > > > > Greetings, > > -- Thiago V. dos Santos > > > >
> Postdoctoral Research Fellow > > Department of Climate and Space Science
> and Engineering > > University of Michigan > > > >
> ___ > > R-sig-Geo mailing
> list > > R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org> > >
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > -- > Dr. Michael
> Sumner > Software and Database Engineer > Australian Antarctic Division
>> 203 Channel Highway > Kingston Tasmania 7050 Australia > >
> [[alternative HTML version deleted]] > >  >
> ___ > R-sig-Geo mailing list >
> R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org> >
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Dr. Michael Sumner
> Software and Database Engineer Australian Antarctic Division 203 Channel
> Highway Kingston Tasmania 7050 Australia [[alternative HTML version
> deleted]] ___ R-sig-Geo
> mailing list R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > [[alternative HTML
> version deleted]] > > ___ >
> R-sig-Geo mailing list > R-sig-Geo@r-project.org
> <mailto:R-sig-Geo@r-project.org> >
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo > -- Dr Christoph
> Reudenbach, Philipps-University of Marburg, Faculty of Geography, GIS
> and Environmental Modeling, Deutschhausstr. 10, D-35032 Marburg, fon:
> ++49.(0)6421.2824296, fax: ++49.(0)6421.2828950, web: gis-ma.org,
> giswerk.org, moc.environmentalinformatics-marburg.de
> ___ R-sig-Geo mailing list
> R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>>     [[alternative HTML version deleted]]
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
> 
> -- 
> Edzer Pebesma
> Institute for Geoinformatics
> Heisenbergstrasse 2, 48151 Muenster, Germany
> Phone: +49 251 8333081
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org <mailto:R-sig-Geo@r-project.org>
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

-- 
Edzer Pebesma
Institute for Geoinformatics
Heisenbergstrasse 2, 48151 Muenster, Germany
Phone: +49 251 8333081
  
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Is the raster package still receiving updates?

[Thiago V. dos Santos via R-sig-Geo]

Dear Edzer,
As a heavy user of both raster and sp (for remote sensing products and climate 
model outputs), I am really interested in the outcomes of your meeting with 
Robert and Etienne.
I was just wondering: eould it be possible/convenient at this point for you to 
share them with us?
Thanks!
 -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan 

On Saturday, October 21, 2017, 5:32:59 AM EDT, Edzer Pebesma 
<edzer.pebe...@uni-muenster.de> wrote:  
 
 The the NOTEs on CRAN and the issues raised on Mike's raster-rforge GH
repo to me suggest that raster has, for a while now, reached a mature
state. Of course, you'd wish that raster development, in terms of
extending its power and features, would continue the pace it did 10
years ago, but that's a different issue. Mike has write access to the
raster r-forge sources, it shouldn't be hard to support Robert with
current issues.

Right after rstudio::conf (Feb 4/5) I'm meeting with Robert Hijmans and
Etienne Racine, to discuss, and work on raster future and integration of
rasters and stars. Anyone interested, feel free to join.

Stars has the ambition to jump over raster, in the sense that it doesn't
want to have the local hard drive as a limitation. That calls for
rethinking quite a few things. So far, things are still all in main
memory but it already integrates with sf:

https://gist.github.com/edzer/381dac079ddcd5174be31209675b3822

and handles time stacks of different attributes in single objects. Feed
back welcome as always; will write up the first blog in a few weeks.


On 10/21/2017 10:41 AM, Roger Bivand wrote:
> Yes, this sounds fruitful. For now, I suggest that those wishing to enhance 
> raster do so in an additional package enhancing raster and depending on it. 
> Immediately, we should help Robert resolve the notes in the tests, they do 
> not affect functionality now, but deserve attention to bring the package into 
> full CRAN compliance. I can help with this.
> 
> Roger
> 
> Roger Bivand
> Norwegian School of Economics
> Bergen, Norway
> 
> 
> 
> Fra: Chris Reudenbach
> Sendt: lørdag 21. oktober, 09.01
> Emne: Re: [R-sig-Geo] Is the raster package still receiving updates?
> Til: r-sig-geo@r-project.org
> 
> 
> Many of us rely on the raster package as a crucial basis for their work. Big 
> thank to Robert.  Nevertheless to contribute needs to focus needs and 
> ressources.  Maybe an information exchange/discussion with Robert about 
> options such as if/how the interested community can support would be more 
> productive. cheers Chris Am 21.10.2017 um 06:52 schrieb Roger Bivand: > The 
> CRAN raster package has not been abandoned, if it had, it would have been 
> reclassified as orphaned. Its test results have notes but no warnings or 
> errors. It would be orphaned if the maintainer did not respond to CRAN 
> requests to resolve test errors. Consequently, you are suggesting a fork, and 
> should rename any package. Only the maintainer may update existing CRAN 
> packages. I would not think that forking an existing CRAN package is the most 
> productive way of contributing to the community. The package has many reverse 
> dependencies, so setting up a test framework to find backwards 
> incompatibilities would be crucial. > > Of course, all contributions are 
> welcome. > > Roger > > Roger Bivand > Norwegian School of Economics > Bergen, 
> Norway > > > > Fra: Michael Sumner > Sendt: lørdag 21. oktober, 06.28 > Emne: 
> Re: [R-sig-Geo] Is the raster package still receiving updates? > Til: 
> Hodgess, Erin > Kopi: R-sig-geo Mailing List > > > Great! I'm keen, but 
> haven't done the work needed to flesh this out properly: 
> https://github.com/mdsumner/raster-rforge/issues I have not checked if 
> there's been any commits to r-forge since I cloned this. Cheers, Mike On Sat, 
> 21 Oct 2017, 12:58 Hodgess, Erin, wrote: > I could take a swat at it, if you 
> wish. Please do send me the list of > fixes. I could start in about a week, 
> if that would work. > Sincerely, > Erin > > > Erin M. Hodgess > Associate 
> Professor > Department of Mathematics and Statistics > University of Houston 
> - Downtown > mailto: hodge...@uhd.edu > > 
>  > From: R-sig-Geo on behalf of 
> Michael > Sumner > Sent: Friday, October 20, 2017 8:05 PM > To: Thiago V. dos 
> Santos > Cc: R-sig-geo Mailing List > Subject: Re: [R-sig-Geo] Is the raster 
> package still receiving updates? > > No, it's effectively abandoned. No 
> concrete plans known. > > I'm interested to keep it going and have kept a 
> list of fixes needed, but > not sure when or if I'll get to it

[R-sig-Geo] Is the raster package still receiving updates?

[Thiago V. dos Santos via R-sig-Geo]

Dear list,

I realized that the latest version of the raster package was released on CRAN 
over a year ago.

I also noticed that Robert's participation in this list has become rather 
scarce.

I was just wondering whether the raster package is still receiving updates? 
 
Greetings,
 -- Thiago V. dos Santos

Postdoctoral Research Fellow
Department of Climate and Space Science and Engineering
University of Michigan

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Multi-file raster stack from netcdf specifying levels

[Thiago V. dos Santos via R-sig-Geo]

Yes, it works!
Thank you very much for the alternative! Greetings, -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan 

On Thursday, October 19, 2017 7:10 PM, Vijay Lulla <vijaylu...@gmail.com> 
wrote:
 

 Maybe you can try something like this?

h2osoil.l1 <- do.call(stack, lapply(files, function(file) raster(file, 
varname='H2OSOI', level=1)))
h2osoil.l2 <- do.call(stack, lapply(files, function(file) raster(file, 
varname='H2OSOI', level=2)))
...
h2osoil.l15 <- do.call(stack, lapply(files, function(file) raster(file, 
varname='H2OSOI', level=15)))


On Thu, Oct 19, 2017 at 6:33 PM, Thiago V. dos Santos via R-sig-Geo 
<r-sig-geo@r-project.org> wrote:

Dear all,

I am working with several netcdf files that contain soil moisture in 15 soil 
layers.

I am trying to create raster stacks of individual layers, but I am getting a 
strange error.

# Get list of files
files <- list.files('~/Downloads', pattern='.nc', full.names=T)

# Stack all the files
h2osoil.l1 <- stack(files, varname='H2OSOI', level=1) # this would be the first 
layer
h2osoil.l2 <- stack(files, varname='H2OSOI', level=2) # this would be the 
second layer

And this is the error I get:

Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘raster’ for signature 
‘"numeric"’
In addition: Warning messages:
1: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)
2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)
3: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)

Here are some files to reproduce this issue (3 MB each file):

https://www.dropbox.com/s/ 9z3fnadqdf1rzuz/cam5clm45_ctl. 
clm2.h0.1980-01.nc?dl=0
https://www.dropbox.com/s/ 0snz58au8ic196p/cam5clm45_ctl. 
clm2.h0.1980-02.nc?dl=0
https://www.dropbox.com/s/ v5ps5bkyrkflki7/cam5clm45_ctl. 
clm2.h0.1980-03.nc?dl=0


I know netcdf files can be quite tricky to work with, but is there any way to 
get this file read in the right way using the raster package? Any other option?


Greetings,
 -- Thiago V. dos Santos

Postdoctoral Research Fellow
Department of Climate and Space Science and Engineering
University of Michigan

__ _
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/ listinfo/r-sig-geo



-- 
Vijay Lulla, Ph.D.
Assistant Professor,Dept. of Geography,Indiana University Purdue University 
Indianapolis (IUPUI)425 University Blvd, CA-207C.Indianapolis, 
in-46202vlu...@iupui.edu

http://vijaylulla.com

   
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Multi-file raster stack from netcdf specifying levels

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I am working with several netcdf files that contain soil moisture in 15 soil 
layers.

I am trying to create raster stacks of individual layers, but I am getting a 
strange error.

# Get list of files
files <- list.files('~/Downloads', pattern='.nc', full.names=T)

# Stack all the files
h2osoil.l1 <- stack(files, varname='H2OSOI', level=1) # this would be the first 
layer
h2osoil.l2 <- stack(files, varname='H2OSOI', level=2) # this would be the 
second layer

And this is the error I get:

Error in (function (classes, fdef, mtable)  : 
unable to find an inherited method for function ‘raster’ for signature 
‘"numeric"’
In addition: Warning messages:
1: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)
2: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)
3: In .rasterObjectFromCDF(x, type = objecttype, band = band, ...) :
"level" set to 1 (there are 15 levels)

Here are some files to reproduce this issue (3 MB each file):

https://www.dropbox.com/s/9z3fnadqdf1rzuz/cam5clm45_ctl.clm2.h0.1980-01.nc?dl=0
https://www.dropbox.com/s/0snz58au8ic196p/cam5clm45_ctl.clm2.h0.1980-02.nc?dl=0
https://www.dropbox.com/s/v5ps5bkyrkflki7/cam5clm45_ctl.clm2.h0.1980-03.nc?dl=0


I know netcdf files can be quite tricky to work with, but is there any way to 
get this file read in the right way using the raster package? Any other option?

 
Greetings,
 -- Thiago V. dos Santos

Postdoctoral Research Fellow
Department of Climate and Space Science and Engineering
University of Michigan

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Optimized rasterOptions() for a (virtually) infinite RAM machine

[Thiago V. dos Santos via R-sig-Geo]

Dear all,
I am using the raster package to process a total of 32 daily climate files 
supplied as netcdf files. Each file is a raster brick with 100 rows x 95 cols x 
54750 time slices and weighs about 1.5 GB.
Essentially, all the processing I am performing on each netcdf file is:
a) to subset a specific date rangeb) to extract values using points
After that, I just convert the extracted data to data.tables and keep working 
in that format.
Since I extract data for about 450 points, and append all the data in a huge 
data.table, I need to use a computer with as much RAM as possible.
I ended up using a spot instance on Amazon EC2. Using an instance with 32 cores 
and 244GB of RAM will cost me around $0.30/hour.
Since I will be charged per hour, I need to optimize my code to get my results 
as fast as possible.
I don't even copy my data to the instance's hard disk; I send the files 
directly to the ram disk (/dev/shm). Even using 48GB of ram disk to store the 
files, I'll still have 196GB of RAM.
Under the scenario of having virtually infinite RAM memory, what would be the 
best rasterOptions() to make sure I am processing all my rasters in memory? Any 
other tips to benefit from such a large amount of RAM?
Thanks, -- Thiago V. dos Santos
Postdoctoral Research FellowDepartment of Climate and Space Science and 
EngineeringUniversity of Michigan
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to objectively subset cities by population

[Thiago V. dos Santos via R-sig-Geo]

Dear Dr. Gräler,
Thanks for your contribution. I very much enjoyed the clustering suggestion, 
and it seems to be available in R's leaflet through the "markerClusterOptions" 
command. It could solve my problem, so I will take a closer look at that.
Regarding your first suggestion, can you point me out some example that uses 
the overlaid grid approach?
Thanks, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Thursday, July 27, 2017, 3:00:55 AM CDT, Dr. Benedikt Gräler 
<b.grae...@52north.org> wrote:

Dear Thiago,

if you want them spatially evenly distributed, you could overlay a grid 
and select the largest per grid box - or maybe more intuitive, select 
the largest per predefined administrative areas (counties/postal 
codes/...). This could also change based on zoom-level. An alternative 
is to group sensors and expand and zoom in by clicking on the group (see 
e.g. [1]).

HTH,

  Ben

[1] http://sensorweb.demo.52north.org/client/#/map


On 27/07/2017 06:09, Thiago V. dos Santos via R-sig-Geo wrote:
> Dear all,
> 
> I have temperature records of nearly 1200 locations in southern Brazil.
> 
> I am writing a shiny app that will show an interactive map with the locations 
> plotted as circles, where the user can click a location to see its 
> temperature time series.
> 
> However, if I show all the locations in the map, it will look really bad, too 
> cramped.
> 
> Therefore, in an attempt to make the map look a bit cleaner, I am trying to 
> think of an objective way to subset the locations. My initial approach would 
> be to show only the "largest" locations, i.e. the ones with a population 
> above a certain threshold.
> 
> The problem is: the distribution of the population is so positively skewed 
> that I am having a hard time determining the optimal cutoff point.
> 
> Does anybody here know any tool or method, possibly spatial, that can assist 
> me with this analysis?
> 
> These are the locations I am working with:
> 
> #---
> # Download and summarize
> locs <- 
> read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1;)
> hist(locs$Population)
> summary(locs$Population)
> 
> # Convert to spatial points and plot
> require(sp)
> coordinates(locs) <- cbind(locs$Lon , locs$Lat)
> plot(locs)
> bubble(locs,"Population")
> #---
> 
> Thanks in advance,
>  -- Thiago V. dos Santos
> 
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>     [[alternative HTML version deleted]]
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Dr. Benedikt Gräler
52°North Initiative for Geospatial Open Source Software GmbH
Martin-Luther-King-Weg 24
48155 Muenster, Germany

E-Mail: b.grae...@52north.org
Fon: +49-(0)-251/396371-39
Fax: +49-(0)-251/396371-11

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

General Managers: Dr. Albert Remke, Dr. Andreas Wytzisk
Local Court Muenster HRB 10849
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to objectively subset cities by population

[Thiago V. dos Santos via R-sig-Geo]

Hi Kent,
Thank you very much for the response.
This is exactly the same approach I am using now, with a different radius 
formula though.
It is not at all terrible, but do you see how overlapped the largest cities 
are, and how difficult it is to click on some smaller cities?
This is why I am looking for an objective way to filter out the smallest 
cities, while at the same time keeping as much information on the map as 
possible (i.e. not leaving too many "blank" areas). 
Best, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Thursday, July 27, 2017, 7:29:44 AM CDT, Kent Johnson  
wrote:


Date: Thu, 27 Jul 2017 04:09:53 + (UTC)
From: "Thiago V. dos Santos" 
To: R-sig-geo Mailing List 
Subject: [R-sig-Geo] How to objectively subset cities by population
Message-ID: <1634627903.292923. 1501128593...@mail.yahoo.com>
Content-Type: text/plain; charset="UTF-8"

Dear all,

I have temperature records of nearly 1200 locations in southern Brazil.

I am writing a shiny app that will show an interactive map with the locations 
plotted as circles, where the user can click a location to see its temperature 
time series.

However, if I show all the locations in the map, it will look really bad, too 
cramped.


Have you considered using leaflet to make an interactive map? This is a good 
start:
library(leaflet)locs <- 
read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1;)
leaflet(locs) %>% addTiles() %>%   addCircles(radius=~20*sqrt(Population), 
label=~as.character(Geocode),              stroke=FALSE, fillOpacity=0.5)
You can configure popups to show HTML or get Shiny events on click, for example 
clicking on a city could display the time series in a separate panel.
Docs and many examples on the leaflet for R web site:
https://rstudio.github.io/leaflet/

Kent
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] How to objectively subset cities by population

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I have temperature records of nearly 1200 locations in southern Brazil.

I am writing a shiny app that will show an interactive map with the locations 
plotted as circles, where the user can click a location to see its temperature 
time series.

However, if I show all the locations in the map, it will look really bad, too 
cramped.

Therefore, in an attempt to make the map look a bit cleaner, I am trying to 
think of an objective way to subset the locations. My initial approach would be 
to show only the "largest" locations, i.e. the ones with a population above a 
certain threshold.

The problem is: the distribution of the population is so positively skewed that 
I am having a hard time determining the optimal cutoff point.

Does anybody here know any tool or method, possibly spatial, that can assist me 
with this analysis?

These are the locations I am working with:

#---
# Download and summarize
locs <- 
read.csv("https://www.dropbox.com/s/ykdd8x1mlc76klt/locations.csv?raw=1;)
hist(locs$Population)
summary(locs$Population)

# Convert to spatial points and plot
require(sp)
coordinates(locs) <- cbind(locs$Lon , locs$Lat)
plot(locs)
bubble(locs,"Population")
#---

Thanks in advance,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Adding text to rasterVis levelplot

[Thiago V. dos Santos via R-sig-Geo]

You can manually define the names that you want and pass them as arguments to 
'names.attr' (the name of each panel) and 'ylab' (the title of the y axis).
Following your example:
col.titles = c('col1','col2','','')
row.titles = c('row1','row2')

levelplot(s, layout=c(2,2),
  names.attr=col.titles,
  ylab=row.titles)

HTH,
Greetings, -- Thiago V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota


On Wednesday, July 19, 2017, 11:13:31 AM CDT, Mauricio Zambrano Bigiarini 
 wrote:

Dear all,

I would like to add custom text labels for the columns and rows of a
levelplot object created with he rasterVis package, similar to the
c("Fall", "Winter", "Spring", Summer") used to label the rows and the
c("a",..., "h") used to label the columns in the plot shown int he
following link:

https://stackoverflow.com/questions/43007241/how-to-add-text-to-a-specific-fixed-location-in-rastervis-levelplot

In particular, how can I add c("col1", "col2") and c("row1", "row2")
as column and row headers, respectively, to the following example
plot:

- START
library(rasterVis)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- stack(r, r+500, r-500, r+200)
levelplot(s, layout=c(2,2), names.att=rep("",4))
- END


Thanks in advance for your help,


Mauricio Zambrano-Bigiarini, PhD

=
Department of Civil Engineering
Faculty of Engineering and Sciences
Universidad de La Frontera, Temuco, Chile
=
"The ultimate inspiration is the deadline"
(Nolan Bushnell)
=
Linux user #454569 -- Linux Mint user

-- 
La información contenida en este correo electrónico y cualquier anexo o 
respuesta relacionada, puede contener datos e información confidencial y no 
puede ser usada o difundida por personas distintas a su(s) destinatario(s). 
Si usted no es el destinatario de esta comunicación, le informamos que 
cualquier divulgación, distribución o copia de esta información constituye 
un delito conforme a la ley chilena. Si lo ha recibido por error, por favor 
borre el mensaje y todos sus anexos y notifique al remitente.

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Specify color for "zero" raster values using levelplot

[Thiago V. dos Santos via R-sig-Geo]

Hi Mel,
Thank you ver much for the suggestion. I've reproduced your palette, which 
looks like exactly what I was looking for. However, the grey color is still not 
associated to zero values in the map.
After plotting my raster with your palette as an argument:
levelplot(annual.mask, cuts=14, col.regions=myPal(15), margin=F)

I ended up with a figure like this: http://i.imgur.com/aqZkCGZ.png, where I was 
hoping to have all white (i.e. zero values) in the map filled with grey.
Any other ideas to achieve that?
This is the file and code to reproduce the figure:
https://dl.dropboxusercontent.com/content_link/jnv58wx0QN8ObKaiqigVckXtDmanOVYlgBxXoiym4oty5MS93xHolItNZ8tJ5gVF/file?dl=1
library(rasterVis)
myPal <- read.table(sep=",", 
text="155,29,32,Brown238,36,37,Tomato238,77,34,Tomato250,143,35,Dark.orange254,215,24,Gold214,223,47,Bitter.lemon188,218,112,Sulu179,179,179,Dark.gray124,199,177,Keppel105,202,229,Viking72,156,211,Curious.blue69,101,173,Chetwode.blue58,86,164,Governor.bay58,72,155,Dark.slate.blue44,46,118,Blue.bell")
myPal <- colorRampPalette(rgb(myPal[, 1:3], names=as.character(myPal$V4), 
maxColorValue=255))
r.annual <- raster("Desktop/r.annual.tif")levelplot(r.annual, cuts=14, 
col.regions=myPal(15), margin=F) Greetings, -- Thiago 
V. dos Santos
PhD studentLand and Atmospheric ScienceUniversity of Minnesota 

On Thursday, May 18, 2017 3:06 AM, Melanie Bacou <m...@mbacou.com> wrote:
 

  Simply use the same palette as in the example: library(lattice)
 
 myPal <- read.table(sep="\t", text="
 155     29     32    Brown
 238     36     37    Tomato
 238 77 34    Tomato
 250    143 35    Dark orange
 254    215     24    Gold
 214    223 47    Bitter lemon
 188    218    112    Sulu
 179    179    179    Dark gray
 124    199    177    Keppel
 105    202    229    Viking
 72 156    211    Curious blue
 69 101    173    Chetwode blue
 58 86 164    Governor bay
 58 72 155    Dark slate blue
 44 46 118    Blue bell")
 
 myPal <- colorRampPalette(rgb(myPal[, 1:3], names=as.character(myPal$V4), 
maxColorValue=255))
 
 x <- seq(pi/4, 5 * pi, length.out = 100)
 y <- seq(pi/4, 5 * pi, length.out = 100)
 r <- as.vector(sqrt(outer(x^2, y^2, "+")))
 grid <- expand.grid(x=x, y=y)
 grid$z <- cos(r^2) * exp(-r/(pi^3))
 levelplot(z~x*y, grid, cuts=14, col.regions=myPal(15), margin=FALSE)
 
 --Mel.
 
 On 05/17/2017 01:25 PM, Thiago V. dos Santos via R-sig-Geo wrote:
  
 Dear all,

I am trying to change the color for zero values in a map produced using 
levelplot to plot a raster file. Specifically, I want to reproduce this figure: 
http://i.imgur.com/mjXxZhO.png, where a red to blue scale is used, but notice 
that zero values have been replaced by grey.

As an example, let's use an adapted version of the August irradiation code from 
the rasterVis webpage:


library(raster)
library(ncdf4)
library(rasterVis)

##Solar irradiation data from CMSAF 
old <- setwd(tempdir())
download.file('https://raw.github.com/oscarperpinan/spacetime-vis/master/data/SISmm2008_CMSAF.zip',
'SISmm2008_CMSAF.zip', method='wget')
unzip('SISmm2008_CMSAF.zip')

listFich <- dir(pattern='\\.nc')
stackSIS <- stack(listFich)
stackSIS <- stackSIS * 24 ##from irradiance (W/m2) to irradiation Wh/m2
idx <- seq(as.Date('2008-01-15'), as.Date('2008-12-15'), 'month')

SISmm <- setZ(stackSIS, idx)
names(SISmm) <- month.abb

setwd(old)

# Set color palette
myTheme=rasterTheme(region=brewer.pal('RdBu', n=11))

Aug <- raster(SISmm, 8)
meanAug <- cellStats(Aug, mean)
levelplot(Aug - meanAug, par.settings = myTheme, margin=FALSE)



In the example above, how can I replace the color of "zero values" with grey?

Thanks,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo




   
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Specify color for "zero" raster values using levelplot

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I am trying to change the color for zero values in a map produced using 
levelplot to plot a raster file. Specifically, I want to reproduce this figure: 
http://i.imgur.com/mjXxZhO.png, where a red to blue scale is used, but notice 
that zero values have been replaced by grey.

As an example, let's use an adapted version of the August irradiation code from 
the rasterVis webpage:


library(raster)
library(ncdf4)
library(rasterVis)

##Solar irradiation data from CMSAF 
old <- setwd(tempdir())
download.file('https://raw.github.com/oscarperpinan/spacetime-vis/master/data/SISmm2008_CMSAF.zip',
'SISmm2008_CMSAF.zip', method='wget')
unzip('SISmm2008_CMSAF.zip')

listFich <- dir(pattern='\\.nc')
stackSIS <- stack(listFich)
stackSIS <- stackSIS * 24 ##from irradiance (W/m2) to irradiation Wh/m2
idx <- seq(as.Date('2008-01-15'), as.Date('2008-12-15'), 'month')

SISmm <- setZ(stackSIS, idx)
names(SISmm) <- month.abb

setwd(old)

# Set color palette
myTheme=rasterTheme(region=brewer.pal('RdBu', n=11))

Aug <- raster(SISmm, 8)
meanAug <- cellStats(Aug, mean)
levelplot(Aug - meanAug, par.settings = myTheme, margin=FALSE)



In the example above, how can I replace the color of "zero values" with grey?
 
Thanks,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] xlsx or rJava issue

[Thiago V. dos Santos via R-sig-Geo]

Didier, I noticed you are using macOS. Are you by any chance using RStudio as 
well? If yes, besides Edzer's suggestion, the following command has also helped 
me in the past:

sudo ln -f -s $(/usr/libexec/java_home)/jre/lib/server/libjvm.dylib 
/usr/local/lib

Then, try re-installing rJava - *from source* - now. It should work.
 Hope this helps,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Friday, March 17, 2017 4:52 AM, Dr Didier G. Leibovici 
 wrote:
Hi,

I had some issues with java settings when using xlsx.

I had re-installed rJava as well.

> library(xlsx)
Error : .onLoad failed in loadNamespace() for 'xlsx', details:
   call: .jinit()
   error: JNI_GetCreatedJavaVMs returned -1

Error: package or namespace load failed for ‘xlsx’
JavaVM: requested Java version ((null)) not available. Using Java at "" 
instead.
JavaVM: Failed to load JVM: /bundle/Libraries/libserver.dylib
JavaVM FATAL: Failed to load the jvm library.

that file on my machine is in:
/Library/Java/JavaVirtualMachines/jdk1.8.0_25.jdk/Contents/Home/jre/lib/server/

is it a particular JAVA_HOME setting or a PATH issue?

I added this in the PATH, re-installed rJava, re-installed xlsx ... nope!
thanks,

Didier

on R version 3.3.2 (2016-10-31)

-- 
Dr Didier G. Leibovici
d-d'yeley-bow-v-c
Senior Research Fellow
Geocomputational Modelling & Geospatial Statistics
Nottingham Geospatial Institute
University of Nottingham, UK
+44 (0)115 84 13924
http://www.nottingham.ac.uk/ngi/people/didier.leibovici
Google+ didier.leibov...@gmail.com
Skype didierleibovici





This message and any attachment are intended solely for the addressee
and may contain confidential information. If you have received this
message in error, please send it back to me, and immediately delete it. 

Please do not use, copy or disclose the information contained in this
message or in any attachment.  Any views or opinions expressed by the
author of this email do not necessarily reflect the views of the
University of Nottingham.

This message has been checked for viruses but the contents of an
attachment may still contain software viruses which could damage your
computer system, you are advised to perform your own checks. Email
communications with the University of Nottingham may be monitored as
permitted by UK legislation.


[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo 

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Calculate anomalies on time-series rasters

[Thiago V. dos Santos via R-sig-Geo]

Thanks for the input, Michael and Loic.

It looks like remote's "anomalize" might work, but the overlay approach was 
more convenient to use based on the shape of my data.
 
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Sunday, February 12, 2017 5:24 PM, Loïc Dutrieux 
<loic.dutri...@conabio.gob.mx> wrote:


On 12/02/2017 03:14, Michael Sumner wrote:
> I believe the "remote" package has functions for doing exactly this.
> 
> HTH
> 
> On Sun, Feb 12, 2017, 17:42 Thiago V. dos Santos via R-sig-Geo <
> r-sig-geo@r-project.org> wrote:
> 
>> Dear all,
>>
>> I have a netcdf file with monthly temperatures values covering the period
>> of January 1961 to December 2010:
>>
>> library(raster)
>> # Create date sequence
>> idx = seq(as.Date("1961/1/1"), as.Date("2010/12/31"), by = "month")
>>
>> # Create raster stack and assign dates
>> r = raster(ncol=5, nrow=5)
>> s = stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
>> s = setZ(s, idx)
>>
>> First, let's consider only the period from 1961-1990 and take the global
>> monthly means (aka climatologies)
>>
>> # Split 1961-1990 period and take climatology
>> s61.90 = subset(s, which(getZ(s)>=as.Date('1961-01-01') &
>> getZ(s)<=as.Date('1990-12-31')))
>> s61.90.mon = zApply(s61.90, by=months, mean, name=month.abb[])
>>
>> s61.90.mon is a raster with 12 layers, one for each monthly mean
>> temperature.
>>
>> Now what I need to do is calculate the monthly 'anomaly' for the remainder
>> period (Jan 1991 to Dev 2010), i.e. the difference of a single month of the
>> period from the corresponding month in the climatological mean.
>>
>> # Here I split the remainder of the time series (1991 to 2010)
>> s91.2010 = subset(s, which(getZ(s)>=as.Date('1991-01-01') &
>> getZ(s)<=as.Date('2010-12-31')))
>>
>> In other words, I need to subtract 'jan' in the raster s.61.90.mon from
>> every 'jan' in the raster s91.2010, and so on for every other month.
>>
>>
>> I was wondering what is the best way to do that using raster (and possibly
>> zoo) functions?

If you're confident about the order of the layers in your rasterStacks
you could use overlay. See the note about recycling in the function help.

# lenght(x) needs to be a multiple of length(y), so that it can be recycled
# e.g. c(10, 10, 10, 10, 10, 10) - c(5, 7)
fun <- function(x, y) {
  x - y
}

annomalies <- overlay(x = s91.2010, y = s61.90.mon, fun = fun)

However, it looks like your climatologies are in alphabetic and not
chronological order, so you'll have to double check that.

Cheers,
Loïc

>>
>> Thank you very much in advance, -- Thiago V. dos Santos
>>
>> PhD candidate
>> Land and Atmospheric Science
>> University of Minnesota
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Calculate anomalies on time-series rasters

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I have a netcdf file with monthly temperatures values covering the period of 
January 1961 to December 2010:

library(raster)
# Create date sequence
idx = seq(as.Date("1961/1/1"), as.Date("2010/12/31"), by = "month")

# Create raster stack and assign dates
r = raster(ncol=5, nrow=5)
s = stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
s = setZ(s, idx)

First, let's consider only the period from 1961-1990 and take the global 
monthly means (aka climatologies)

# Split 1961-1990 period and take climatology
s61.90 = subset(s, which(getZ(s)>=as.Date('1961-01-01') & 
getZ(s)<=as.Date('1990-12-31')))
s61.90.mon = zApply(s61.90, by=months, mean, name=month.abb[])

s61.90.mon is a raster with 12 layers, one for each monthly mean temperature.

Now what I need to do is calculate the monthly 'anomaly' for the remainder 
period (Jan 1991 to Dev 2010), i.e. the difference of a single month of the 
period from the corresponding month in the climatological mean.

# Here I split the remainder of the time series (1991 to 2010)
s91.2010 = subset(s, which(getZ(s)>=as.Date('1991-01-01') & 
getZ(s)<=as.Date('2010-12-31')))

In other words, I need to subtract 'jan' in the raster s.61.90.mon from every 
'jan' in the raster s91.2010, and so on for every other month.


I was wondering what is the best way to do that using raster (and possibly zoo) 
functions?

Thank you very much in advance, -- Thiago V. dos Santos

PhD candidate
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Calculate the area of pixel in square kilometers from a raster in lat/long projection

[Thiago V. dos Santos via R-sig-Geo]

raster::area would calculate the area in square kilometers, and from there one 
can convert to hectare or square meters.

I was wondering: how do those calculations differ from raster::area's? 
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota


On Tuesday, November 22, 2016 8:15 PM, Isaque Daniel 
 wrote:



Hi there,
I found a example in stackoverflow of distance calculation and convertion, so I 
modificate that to compute the area by pixel in hectares using the raster in 
degrees.

#
# Function to compute the area  in hectares

measureAreahec <- function(imageRaster) {
  R <- 6378.137# radius of earth in Km
  areaImagetemp <- area(ceiImg)
  xresVal <- xres(areaImagetemp)
  yresVal <- yres(areaImagetemp)
  coords <- data.frame(coordinates(areaImagetemp))
  colnames(coords) <- c("lon","lat")
  lon1 <- coords$lon - (xresVal/2)
  lon2 <- coords$lon + (xresVal/2)
  lat1 <- coords$lat - (yresVal/2)
  lat2 <- coords$lat + (yresVal/2)

  # width
  dLat <- (lat1-lat1)*pi/180
  dLon <- (lon2-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlongitude <- R * c * 1000

  # heigth
  dLat <- (lat1-lat2)*pi/180
  dLon <- (lon1-lon1)*pi/180
  a <- sin((dLat/2))^2 + cos(lat1*pi/180)*cos(lat1*pi/180)*(sin(dLon/2))^2
  c <- 2 * atan2(sqrt(a), sqrt(1-a))
  dlatitude <- R * c * 1000
  return (dlatitude * dlongitude / 1)# distance 
in meters
}
##


It's a very simple one, but ,maybe can help somebody.

Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazlia University (UNB)



De: R-sig-Geo  em nome de Isaque Daniel 

Enviado: ter�-feira, 22 de novembro de 2016 17:48
Para: r-sig-geo
Assunto: [R-sig-Geo] Calculate the area of pixel in square kilometers from a 
raster in lat/long projection

Hi,


I'm looking for some advice to convert the area calculated by raster::area() in 
a raster image as the projection in lat/long.

[[elided Hotmail spam]]
Best
Isaque


--
Agronomist engineer
Master in Remote Sensing - National  Institute for Space Research (INPE) - 
Brazil
PHD Student in Transport - Brazilia University (UNB)

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
R-sig-Geo Info Page
stat.ethz.ch
R-sig-Geo -- R Special Interest Group on using Geographical data and Mapping 
About R-sig-Geo





[[alternative HTML version deleted]]


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] modis download authentication

[Thiago V. dos Santos via R-sig-Geo]

This is what worked for me, on a Mac:

myopts <- RCurl::curlOptions(netrc=TRUE, netrc.file=path.expand("~/.netrc"),
 cookiefile=path.expand("~/.urs_cookies"),
 followlocation=TRUE)


getBinaryURL(your_url, .opts=myopts)

where the .netrc file should look like this:

machine urs.earthdata.nasa.gov login your_login password your_pass
machine e4ftl01.cr.usgs.gov login your_login password your_pass

(note that both lines are required, and your_login and your_pass are your 
credentials to the EarthData system)

and the .urs_cookies is just an empty file where the access information will be 
saved to.

You can use a simple command that explicitly requires your credentials, but I 
don't think revealing your credentials on a R script is a good idea.

There are some helpful tips here as well:
http://disc.sci.gsfc.nasa.gov/registration/registration-for-data-access
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Thursday, September 8, 2016 6:19 AM, Luca Candeloro 
 wrote:
I used succesfully the function getBinaryURL from RCurl library to download
.hdf file, as in:

content <- getBinaryURL("
http://e4ftl01.cr.usgs.gov/MOLT/MOD11A2.005/2016.01.01/MOD11A2.A2016001.h18v04.005.2016012041335.hdf
")

Since http://e4ftl01.cr.usgs.gov/ has recently applied a login and pass
authentication to retrieve data, the previous code doesn't work.
It seems (from
https://lpdaac.usgs.gov/sites/default/files/public/get_data/docs/Command%20Line%20Access%20Tips%20for%20Utilizing%20Earthdata%20Login.docx
)
that curl options should be necessarily  set.
Does someone know how to correctly specify RCurl options in R, so that
getBinaryURL works again?
thanks,
Luca

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to calculate seasonal means in raster time series?

[Thiago V. dos Santos via R-sig-Geo]

Thanks Loïc,

I was almost there, but I was definitely missing "as.yearmon".
 
It works now.
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Wednesday, August 17, 2016 10:37 AM, Loïc Dutrieux <loic.dutri...@wur.nl> 
wrote:


On 17/08/2016 02:51, Bacou, Melanie wrote:
> Maybe create a function that takes a date as input and returns a
> meteorological season, and pass this function to `zapply(by=fun)`.
> --Mel.
>
> On 8/16/2016 8:38 PM, Thiago V. dos Santos via R-sig-Geo wrote:
>> Dear all,
>>
>> I am using the raster package to calculate seasonal averages of
>> climatic variables. I usually use zApply to perform temporal
>> calculations on gridded time series.
>>
>> However, this time I need to calculate seasonal averages according to
>> the meteorological nomenclature, with DJF (=winter in the southern
>> hemisphere: December, January, February), MAM, JJA, and SON. It is
>> similar to zoo's yearqtr, but with a offset of -1 month on each quarter.
>>
>> That means that December values comes from the previous year.
>> This is a typical raster I work with:
>>
>> library(raster)
>>
>> # Create date sequence
>> idx <- seq(as.Date("2001/1/1"), as.Date("2010/12/31"), by = "day")
>>
>> # Create raster stack and assign dates
>> r <- raster(ncol=5, nrow=5)
>> s <- stack(lapply(1:length(idx), function(x) setValues(r,
>> runif(ncell(r)
>> s <- setZ(s, idx)

Hi Thiago,

You can try to shift the aggregation indices.

library(zoo)
zApply(s, by = as.yearqtr(as.yearmon(idx) + 1/12), fun = mean)

Kind regards,
Loïc

>>
>>
>> On a raster stack like this, what would be the best strategy to
>> calculate "meteorological" seasonal means?
>>   Thanks in advance,
>>   -- Thiago V. dos Santos
>>
>> PhD student
>> Land and Atmospheric Science
>> University of Minnesota
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] How to calculate seasonal means in raster time series?

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I am using the raster package to calculate seasonal averages of climatic 
variables. I usually use zApply to perform temporal calculations on gridded 
time series.

However, this time I need to calculate seasonal averages according to the 
meteorological nomenclature, with DJF (=winter in the southern hemisphere: 
December, January, February), MAM, JJA, and SON. It is similar to zoo's 
yearqtr, but with a offset of -1 month on each quarter.

That means that December values comes from the previous year.
This is a typical raster I work with:

library(raster)

# Create date sequence
idx <- seq(as.Date("2001/1/1"), as.Date("2010/12/31"), by = "day")

# Create raster stack and assign dates
r <- raster(ncol=5, nrow=5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
s <- setZ(s, idx)


On a raster stack like this, what would be the best strategy to calculate 
"meteorological" seasonal means?
 Thanks in advance,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Incorrect month order in zApply function

[Thiago V. dos Santos via R-sig-Geo]

Thanks for the information, Oscar. Perhaps it would be a good idea to include 
it in the zApply description within the raster manual? 
Best, Thiago. 

Sent from Yahoo Mail on Android 
 
  On Thu, Aug 4, 2016 at 3:19 AM, Oscar Perpiñán 
Lamigueiro wrote:   Hello,

With the current code of zApply, you should use a function that returns an 
ordered vector. For example:

Month <- function(x)
{
    idx <- as.numeric(format(x, '%m'))
    factor(month.name[idx],
          levels = month.name,
          ordered = TRUE)
}


x2 <- zApply(s, by = Month, mean)

Best,

Oscar.
-- 
Oscar Perpiñán Lamigueiro
Dpto. Ing. Eléctrica, Electrónica, Automática y Física Aplicada
Escuela Técnica Superior de Ingeniería y Diseño Industrial
URL: http://oscarperpinan.github.io  

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Incorrect month order in zApply function

[Thiago V. dos Santos via R-sig-Geo]

Thanks everyone for their suggestions.

Reordering the output of zApply would be a WRONG solution, because the layers 
are already ordered from Jan to Dec.

It is just the NAMES that are misaligned, which can cause a bit of confusion.

I already contacted Robert about this issue and we should see a solution soon.
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Friday, July 29, 2016 2:09 PM, Kay Kilpatrick <kkilpatr...@rsmas.miami.edu> 
wrote:
Hi

The various Apply functions always return the output in simple 
alphabetical order when the factor label is a character, regardless of 
the order of the factor in the calling argument.

Ben had a correct solution --- simply reorder the output from Zapply.

Alternatively you could use a numeric representation for the month label 
using by = as.numeric(months(x). The output would then be in increasing 
numeric order from 1 -12.

I do not believe Vijay's suggestion is proper. It does not reorder the 
output --- instead it changes the label of a given monthly mean eg. The 
value for Jan. would then be labeled as Apr.

k

On 7/28/16 10:52 PM, Thiago V. dos Santos via R-sig-Geo wrote:
> Dear all,
>
> I am using the raster package to calculate monthly averages of climatic 
> variables.
>
> Essentially, this is the function I use:
>
> library(raster)
>
> # Create date sequence
> idx <- seq(as.Date("1996/1/1"), as.Date("2010/12/31"), by = "day")
>
> # Create raster stack and assign dates
> r <- raster(ncol=5, nrow=5)
> s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
> s <- setZ(s, idx)
>
> # Calculate monthly average
> x <- zApply(s, by=months, mean, name=month.abb[])
>
> names(x)
> [1] "April" "August" "December" "February" "January" "July" "June"
> [8] "March" "May" "November" "October" "September"
> getZ(x)
> [1] "April" "August" "December" "February" "January" "July" "June"
> [8] "March" "May" "November" "October" "September"
>
>
> The problem here is the output of both names(x) and getZ(x). It looks like a 
> random month order is returned (even though I provide the labels), which 
> makes me confused about the results.
>
>
> By performing the same calculation in a different software and comparing the 
> results, I came to realize that the order of months for the results by raster 
> should, in fact, be Jan-Feb-Mar-Apr-May-Jun-Jul-Aug-Sep-Oct-Nov-Dec
>  
> How can I control the way raster delivers the object names after using 
> zApply, in order to sort the months in the "natural" order?
>
> Greetings, -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Dsig-2Dgeo=CwICAg=y2w-uYmhgFWijp_IQN0DhA=VXzrzLaooo912e5bVVEO_z5cHQXu_22lipAXBkfCXjM=tzdb3kN3Nxjo7NbPkPVqtF269bJen1bf3AzSHEwd7CA=zuokuoKIDT6cOaBD2J6ePy1klmf1J3-f3BynP_G5I-g=


-- 
Katherine (Kay) Kilpatrick
University of Miami/RSMAS
Ocean Science Department
Satellite Remote Sensing Laboratory
MSC 231
4600 Rickenbacker Cswy
Miami, Fl
ph: 305-962-3069
kkilpatr...@rsmas.miami.edu


___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to calculate climatology in rasterbricks

[Thiago V. dos Santos via R-sig-Geo]

Thanks again Loïc and Vijay for helping to identify the issue.

Hopefully Robert or any of the contributors will see this question and kindly 
suggest a workaround.
 Cheers,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Saturday, June 4, 2016 10:30 AM, Vijay Lulla <vijaylu...@gmail.com> wrote:
Yes, it is because of na.rm argument.  It's stackApply that defaults
na.rm=TRUE not zApply.  You can check that from the following
interaction:

> args(mean.default)
function (x, trim = 0, na.rm = FALSE, ...)
NULL
> args(sum)
function (..., na.rm = FALSE)
NULL
> args(stackApply)
function (x, indices, fun, filename = "", na.rm = TRUE, ...)
NULL
> args(zApply)
function (x, by, fun = mean, name = "", ...)
NULL
>

It is surprising that stackApply defaults differently than most of the
R functions but it kinda makes sense too.  I suspect the defaults are
what they are because in remote sensing related image procesing work
NAs get introduced, especially when doing projection transformation
and this can trip up many users and the authors might've decided to be
helpful to them.  I don't know...I'm just speculating.

HTH,
Vijay.


On Sat, Jun 4, 2016 at 9:48 AM, Loïc Dutrieux <loic.dutri...@wur.nl> wrote:
> Hi Thiago,
>
> Yes, I can reproduce.
> It's probably because zApply() defaults to na.rm = TRUE (in fact the na.rm =
> argument is ignored if you pass it), and
>
> sum(c(NA, NA), na.rm = TRUE)
> [1] 0 (I was surprised by that).
>
> Perhaps the stackApply() call in the zApply() function definition is missing
> an ellipsis.
>
> Cheers,
> Loïc
>
>
> On 06/03/2016 10:25 PM, Thiago V. dos Santos wrote:
>>
>> Dear Vijay and Loïc,
>>
>>
>> I have identified some artifacts in the results when using zApply on my
>> data. To reproduce it, please download this shapefile
>> (https://www.dropbox.com/s/in2z10mlerr2sfu/southern.zip?dl=0) and run the
>> code below (see the comments after the plot commands):
>>
>> ---
>> library(raster)
>> library(zoo)
>>
>> # Create the date sequence
>> idx <- seq(as.Date("1961/1/1"), as.Date("1990/12/31"), by = "day")
>>
>> # Create raster stack and assign dates
>> r <- raster(ncol=24, nrow=21, xmn=-58, xmx=-47.5, ymn=-34, ymx=-22,
>> resolution=0.5)
>> s <- stack(lapply(1:length(idx), function(x) setValues(r,
>> runif(ncell(r)
>> s <- setZ(s, idx)
>>
>> # Load shapefile and crop data
>> shp <- shapefile('~/Downloads/southern.shp', warnPRJ=F)
>> s.c <- crop(s, extent(shp))
>> s.c <- mask(s.c, shp)
>> plot(s.c,1)
>> plot(shp,add=T) # Looks good
>>
>> # Loïc's approach using zApply
>> sYM <- zApply(s.c, by = as.yearmon, sum)
>> sM <- zApply(sYM, by = months, mean)
>> plot(sYM,1)  # Note that the "empty" space across the extent was filled
>> with 0's
>> plot(sM,1)  # The same "gray area" here
>>
>> # Vijay's approach using calc
>> idxYM <- as.integer(strftime(idx,"%Y%m"))
>> idxM <- unique(idxYM)%%100
>> meanYM <- calc(s.c, fun=function(x) { by(x, idxYM, sum) })
>> meanM <- calc(meanYM, fun=function(x) { by(x, idxM, mean) })
>> plot(meanYM,1)  # Note that no 0 was added across the extent
>> plot(meanM,1)  # This is OK too, no grey area
>> ---
>>
>> Basically, a portion of the raster extent was filled with 0's after using
>> zApply. It doesn't happen when using calc.
>>
>> Any ideas on what can be causing this?
>>
>>   Greetings,
>>   -- Thiago V. dos Santos
>>
>> PhD student
>> Land and Atmospheric Science
>> University of Minnesota
>>
>>
>>
>> On Friday, June 3, 2016 1:47 PM, Vijay Lulla <vijaylu...@gmail.com> wrote:
>> Looking at Loïc's answer and from reading ?zApply I learned of
>> stackApply, which is used internally by zApply.  While zApply is for
>> time series of layers, stackApply is more general and can be used for
>> applying functions to groups of layers in a stack/brick.  It is very
>> similar to base R's `tapply`.  And it is also fast!
>>
>>> system.time(meanYM1 <- stackApply(s, idxYM, mean))
>>
>> user  system elapsed
>> 0.450.030.48
>>>
>>> system.time(meanYM <- calc(s, fun=function(x) { by(x, idxYM, mean)}))
>>
>> user  system elapsed
>>17.500.01   17.61
>>>
>>> all(meanYM[] == meanYM1[])
>>
>> [1] TRUE
>>
>> Thanks Loïc for pointing out zApply.
>>
>>
>&

Re: [R-sig-Geo] How to calculate climatology in rasterbricks

[Thiago V. dos Santos via R-sig-Geo]

Dear Vijay and Loïc,


I have identified some artifacts in the results when using zApply on my data. 
To reproduce it, please download this shapefile 
(https://www.dropbox.com/s/in2z10mlerr2sfu/southern.zip?dl=0) and run the code 
below (see the comments after the plot commands):

---
library(raster)
library(zoo)

# Create the date sequence
idx <- seq(as.Date("1961/1/1"), as.Date("1990/12/31"), by = "day")

# Create raster stack and assign dates
r <- raster(ncol=24, nrow=21, xmn=-58, xmx=-47.5, ymn=-34, ymx=-22, 
resolution=0.5)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
s <- setZ(s, idx)

# Load shapefile and crop data
shp <- shapefile('~/Downloads/southern.shp', warnPRJ=F)
s.c <- crop(s, extent(shp))
s.c <- mask(s.c, shp)
plot(s.c,1)
plot(shp,add=T) # Looks good

# Loïc's approach using zApply
sYM <- zApply(s.c, by = as.yearmon, sum)
sM <- zApply(sYM, by = months, mean)
plot(sYM,1)  # Note that the "empty" space across the extent was filled with 0's
plot(sM,1)  # The same "gray area" here

# Vijay's approach using calc
idxYM <- as.integer(strftime(idx,"%Y%m"))
idxM <- unique(idxYM)%%100
meanYM <- calc(s.c, fun=function(x) { by(x, idxYM, sum) })
meanM <- calc(meanYM, fun=function(x) { by(x, idxM, mean) })
plot(meanYM,1)  # Note that no 0 was added across the extent
plot(meanM,1)  # This is OK too, no grey area
---

Basically, a portion of the raster extent was filled with 0's after using 
zApply. It doesn't happen when using calc.

Any ideas on what can be causing this?

 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Friday, June 3, 2016 1:47 PM, Vijay Lulla <vijaylu...@gmail.com> wrote:
Looking at Loïc's answer and from reading ?zApply I learned of
stackApply, which is used internally by zApply.  While zApply is for
time series of layers, stackApply is more general and can be used for
applying functions to groups of layers in a stack/brick.  It is very
similar to base R's `tapply`.  And it is also fast!

> system.time(meanYM1 <- stackApply(s, idxYM, mean))
   user  system elapsed
   0.450.030.48
> system.time(meanYM <- calc(s, fun=function(x) { by(x, idxYM, mean)}))
   user  system elapsed
  17.500.01   17.61
> all(meanYM[] == meanYM1[])
[1] TRUE

Thanks Loïc for pointing out zApply.


On Fri, Jun 3, 2016 at 11:30 AM, Thiago V. dos Santos via R-sig-Geo
<r-sig-geo@r-project.org> wrote:
> Cool Loïc, thanks for showing one more option.
>
> By the way, in this case zApply is surprisingly faster than calc:
>
>> system.time(meanYM <- calc(s,fun=function(x) { by(x, idxYM, sum) }))
> user  system elapsed
> 21.700   0.248  22.095
>
>
>> system.time(sYM <- zApply(s, by = as.yearmon, sum))
> user  system elapsed
> 0.811   0.047   0.866
>
>  Cheers,
>  -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
>
>
> On Friday, June 3, 2016 4:25 AM, Loïc Dutrieux <loic.dutri...@wur.nl> wrote:
> This can also be done with zApply:
>
> library(zoo)
>
> sYM <- zApply(s, by = as.yearmon, sum)
> sM <- zApply(sYM, by = months, mean)
>
> Cheers,
> Loïc
>
> On 06/03/2016 02:02 AM, Vijay Lulla wrote:
>> I think the following StackOverflow question has the answer:
>> http://stackoverflow.com/questions/16135877/applying-a-function-to-a-multidimensional-array-with-grouping-variable/16136775#16136775
>>
>> Following the instructions listed on that page for your case might go
>> something like below:
>>
>>> idxYM <- as.integer(strftime(idx,"%Y%m"))
>>> idxM <- unique(idxYM)%%100
>>> meanYM <- calc(s,fun=function(x) { by(x, idxYM, mean) })
>>> meanYM
>> class   : RasterBrick
>> dimensions  : 20, 20, 400, 360  (nrow, ncol, ncell, nlayers)
>> resolution  : 18, 9  (x, y)
>> extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
>> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
>> data source : in memory
>> names   : X196101, X196102, X196103, X196104, X196105, X196106,
>> X196107, X196108, X196109, X196110, X196111, X196112, X196201,
>> X196202, X196203, ...
>> min values  :  0.3728,  0.2725,  0.3421,  0.3652,  0.3342,  0.3185,
>> 0.3130,  0.3780,  0.3376,  0.3727,  0.3537,  0.3737,  0.3515,  0.3588,
>>   0.3334, ...
>> max values  :  0.6399,  0.6652,  0.6583,  0.6640,  0.6359,  0.6761,
>> 0.6442,  0.6800,  0.6397,  0.6769,  0.6489,  0.6388,  0.6471,  0.6661,
>>   0.6255, ...
>>
>>> meanM <- calc(meanYM, fun=function(x) { by(x, idxM, mean) })
>>> meanM
>> c

Re: [R-sig-Geo] How to calculate climatology in rasterbricks

[Thiago V. dos Santos via R-sig-Geo]

Cool Loïc, thanks for showing one more option.

By the way, in this case zApply is surprisingly faster than calc:

> system.time(meanYM <- calc(s,fun=function(x) { by(x, idxYM, sum) }))
user  system elapsed 
21.700   0.248  22.095 


> system.time(sYM <- zApply(s, by = as.yearmon, sum))
user  system elapsed 
0.811   0.047   0.866

 Cheers,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Friday, June 3, 2016 4:25 AM, Loïc Dutrieux <loic.dutri...@wur.nl> wrote:
This can also be done with zApply:

library(zoo)

sYM <- zApply(s, by = as.yearmon, sum)
sM <- zApply(sYM, by = months, mean)

Cheers,
Loïc

On 06/03/2016 02:02 AM, Vijay Lulla wrote:
> I think the following StackOverflow question has the answer:
> http://stackoverflow.com/questions/16135877/applying-a-function-to-a-multidimensional-array-with-grouping-variable/16136775#16136775
>
> Following the instructions listed on that page for your case might go
> something like below:
>
>> idxYM <- as.integer(strftime(idx,"%Y%m"))
>> idxM <- unique(idxYM)%%100
>> meanYM <- calc(s,fun=function(x) { by(x, idxYM, mean) })
>> meanYM
> class   : RasterBrick
> dimensions  : 20, 20, 400, 360  (nrow, ncol, ncell, nlayers)
> resolution  : 18, 9  (x, y)
> extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : in memory
> names   : X196101, X196102, X196103, X196104, X196105, X196106,
> X196107, X196108, X196109, X196110, X196111, X196112, X196201,
> X196202, X196203, ...
> min values  :  0.3728,  0.2725,  0.3421,  0.3652,  0.3342,  0.3185,
> 0.3130,  0.3780,  0.3376,  0.3727,  0.3537,  0.3737,  0.3515,  0.3588,
>   0.3334, ...
> max values  :  0.6399,  0.6652,  0.6583,  0.6640,  0.6359,  0.6761,
> 0.6442,  0.6800,  0.6397,  0.6769,  0.6489,  0.6388,  0.6471,  0.6661,
>   0.6255, ...
>
>> meanM <- calc(meanYM, fun=function(x) { by(x, idxM, mean) })
>> meanM
> class   : RasterBrick
> dimensions  : 20, 20, 400, 12  (nrow, ncol, ncell, nlayers)
> resolution  : 18, 9  (x, y)
> extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
> data source : in memory
> names   : X1, X2, X3, X4, X5, X6, X7,
> X8, X9,X10,X11,X12
> min values  : 0.4645, 0.4715, 0.4768, 0.4717, 0.4749, 0.4705, 0.4697,
> 0.4724, 0.4629, 0.4774, 0.4736, 0.4708
> max values  : 0.5274, 0.5275, 0.5293, 0.5259, 0.5285, 0.5276, 0.5269,
> 0.5260, 0.5256, 0.5281, 0.5279, 0.5286
>
>>
>
> I'm not sure how [in]efficient this is for actual (i.e. not toy
> example) data.  Maybe others more experienced, and knowledgeable,
> members can provide better answers.
>
> HTH,
> Vijay.
>
> On Thu, Jun 2, 2016 at 4:30 PM, Thiago V. dos Santos via R-sig-Geo
> <r-sig-geo@r-project.org> wrote:
>> Dear all,
>>
>> I am working with daily time series of meteorological variables. This is an 
>> example of the dataset:
>>
>> library(raster)
>>
>> # Create date sequence
>> idx <- seq(as.Date("1961/1/1"), as.Date("1990/12/31"), by = "day")
>>
>> # Create raster stack and assign dates
>> r <- raster(ncol=20, nrow=20)
>> s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
>> s <- setZ(s, idx)
>>
>>
>> Now, let's assume those values represent daily precipitation. What I need to 
>> do is to integrate daily to monthly values,
>> and then take a monthly climatology. Climatology in this case means 
>> multi-year average of selected months, e.g., an average of the 30 Octobers 
>> from 1961 to 1990, an average of the 30 Novembers from 1961 to 1990 and etc.
>>
>> On the other hand, let's assume the raster values represent daily 
>> temperature. Integrating daily to monthly temperature doesn't make sense. 
>> Hence, instead of integrating daily values, I need to take monthly means 
>> (e.g. mean value of all days in every month), and then calculate the 
>> climatology.
>>
>> What would be the best approach to achieve that using the raster package?
>>
>>   Greetings,
>>   -- Thiago V. dos Santos
>>
>> PhD student
>> Land and Atmospheric Science
>> University of Minnesota
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] How to calculate climatology in rasterbricks

[Thiago V. dos Santos via R-sig-Geo]

Dear Vijay,

Thank you so much for your answer - it was exactly what I was looking for.

In fact, I compared your code (applied to my actual data) to the code I wrote 
in cdo (a tool to deal with climate data in netcdf format) and the results 
match to the eleventh decimal place. With the obvious advantage of not breaking 
my workflow in R.
 Cheers,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Thursday, June 2, 2016 7:02 PM, Vijay Lulla <vijaylu...@gmail.com> wrote:
I think the following StackOverflow question has the answer:
http://stackoverflow.com/questions/16135877/applying-a-function-to-a-multidimensional-array-with-grouping-variable/16136775#16136775

Following the instructions listed on that page for your case might go
something like below:

> idxYM <- as.integer(strftime(idx,"%Y%m"))
> idxM <- unique(idxYM)%%100
> meanYM <- calc(s,fun=function(x) { by(x, idxYM, mean) })
> meanYM
class   : RasterBrick
dimensions  : 20, 20, 400, 360  (nrow, ncol, ncell, nlayers)
resolution  : 18, 9  (x, y)
extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names   : X196101, X196102, X196103, X196104, X196105, X196106,
X196107, X196108, X196109, X196110, X196111, X196112, X196201,
X196202, X196203, ...
min values  :  0.3728,  0.2725,  0.3421,  0.3652,  0.3342,  0.3185,
0.3130,  0.3780,  0.3376,  0.3727,  0.3537,  0.3737,  0.3515,  0.3588,
0.3334, ...
max values  :  0.6399,  0.6652,  0.6583,  0.6640,  0.6359,  0.6761,
0.6442,  0.6800,  0.6397,  0.6769,  0.6489,  0.6388,  0.6471,  0.6661,
0.6255, ...

> meanM <- calc(meanYM, fun=function(x) { by(x, idxM, mean) })
> meanM
class   : RasterBrick
dimensions  : 20, 20, 400, 12  (nrow, ncol, ncell, nlayers)
resolution  : 18, 9  (x, y)
extent  : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names   : X1, X2, X3, X4, X5, X6, X7,
   X8, X9,X10,X11,X12
min values  : 0.4645, 0.4715, 0.4768, 0.4717, 0.4749, 0.4705, 0.4697,
0.4724, 0.4629, 0.4774, 0.4736, 0.4708
max values  : 0.5274, 0.5275, 0.5293, 0.5259, 0.5285, 0.5276, 0.5269,
0.5260, 0.5256, 0.5281, 0.5279, 0.5286

>

I'm not sure how [in]efficient this is for actual (i.e. not toy
example) data.  Maybe others more experienced, and knowledgeable,
members can provide better answers.

HTH,
Vijay.


On Thu, Jun 2, 2016 at 4:30 PM, Thiago V. dos Santos via R-sig-Geo
<r-sig-geo@r-project.org> wrote:
> Dear all,
>
> I am working with daily time series of meteorological variables. This is an 
> example of the dataset:
>
> library(raster)
>
> # Create date sequence
> idx <- seq(as.Date("1961/1/1"), as.Date("1990/12/31"), by = "day")
>
> # Create raster stack and assign dates
> r <- raster(ncol=20, nrow=20)
> s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
> s <- setZ(s, idx)
>
>
> Now, let's assume those values represent daily precipitation. What I need to 
> do is to integrate daily to monthly values,
> and then take a monthly climatology. Climatology in this case means 
> multi-year average of selected months, e.g., an average of the 30 Octobers 
> from 1961 to 1990, an average of the 30 Novembers from 1961 to 1990 and etc.
>
> On the other hand, let's assume the raster values represent daily 
> temperature. Integrating daily to monthly temperature doesn't make sense. 
> Hence, instead of integrating daily values, I need to take monthly means 
> (e.g. mean value of all days in every month), and then calculate the 
> climatology.
>
> What would be the best approach to achieve that using the raster package?
>
>  Greetings,
>  -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] How to calculate climatology in rasterbricks

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I am working with daily time series of meteorological variables. This is an 
example of the dataset:

library(raster)

# Create date sequence
idx <- seq(as.Date("1961/1/1"), as.Date("1990/12/31"), by = "day")

# Create raster stack and assign dates
r <- raster(ncol=20, nrow=20)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
s <- setZ(s, idx)


Now, let's assume those values represent daily precipitation. What I need to do 
is to integrate daily to monthly values, 
and then take a monthly climatology. Climatology in this case means multi-year 
average of selected months, e.g., an average of the 30 Octobers from 1961 to 
1990, an average of the 30 Novembers from 1961 to 1990 and etc.

On the other hand, let's assume the raster values represent daily temperature. 
Integrating daily to monthly temperature doesn't make sense. Hence, instead of 
integrating daily values, I need to take monthly means (e.g. mean value of all 
days in every month), and then calculate the climatology.

What would be the best approach to achieve that using the raster package?

 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Run focal function on a multi-layer raster

[Thiago V. dos Santos via R-sig-Geo]

Hi Loïc,

Thanks for the hint - it worked like a charm!
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Saturday, May 21, 2016 3:54 PM, Loïc Dutrieux <loic.dutri...@wur.nl> wrote:
Hi Thiago,

calc is not aware of neighbouring pixels in the x,y dimensions, so 
functions like focal cannot work.

I found the function below in an old repos of mine. It still seems to 
work fine :)

#' Focal for RasterBrick or RasterStack
#'
#' @param x RasterBrick/Stack or character pointing to multilayer raster 
object written on disk
#' @param w See \code{\link{focal}}
#' @param ... Arguments to be passed to \code{\link{focal}}
#'
#' @import raster
#' @export
#'

multiFocal <- function(x, w=matrix(1, nr=3, nc=3), ...) {

   if(is.character(x)) {
 x <- brick(x)
   }
   # The function to be applied to each individual layer
   fun <- function(ind, x, w, ...){
 focal(x[[ind]], w=w, ...)
   }

   n <- seq(nlayers(x))
   list <- lapply(X=n, FUN=fun, x=x, w=w, ...)

   out <- stack(list)
   return(out)
}

On 05/21/2016 10:02 PM, Thiago V. dos Santos via R-sig-Geo wrote:
> I have a multi-layer raster, in this case a rasterbrick, on which I would 
> like to apply a focal function on each layer, and get another rasterbrick as 
> a result.
>
> I am trying to use calc, but apparently I am not assembling my function in 
> the proper way:
>
> require(raster)
> r <- raster(ncol=50, nrow=50)
> r[]=1:ncell(r)
> b <- brick(r,r,r,r,r,r)
> b <- b * 1:6
>

multiFocal(b, w=matrix(1, 5, 5), mean)

Cheers,
Loïc

> myfocalfun <- function(x){
> b.f <- focal(x, w=matrix(1, 5, 5), mean)
> return(b.f)
> }
>
> b1 <- calc(b, myfocalfun)
>
> Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
> cannot use this function
>
> What am I missing here?
>   Greetings,
>   -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Run focal function on a multi-layer raster

[Thiago V. dos Santos via R-sig-Geo]

I have a multi-layer raster, in this case a rasterbrick, on which I would like 
to apply a focal function on each layer, and get another rasterbrick as a 
result.

I am trying to use calc, but apparently I am not assembling my function in the 
proper way:

require(raster)
r <- raster(ncol=50, nrow=50)
r[]=1:ncell(r)
b <- brick(r,r,r,r,r,r)
b <- b * 1:6

myfocalfun <- function(x){
b.f <- focal(x, w=matrix(1, 5, 5), mean)
return(b.f)
}

b1 <- calc(b, myfocalfun)

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
cannot use this function

What am I missing here?
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Warning when stacking, masking or cropping raster objects

[Thiago V. dos Santos via R-sig-Geo]

I also noticed those warnings after the last raster update. Here is one simple 
way to reproduce them:


library(raster)
r <- raster(nrows=100, ncols=100, xmn=5, xmx=10, ymn=46, ymx=51)
r[] <- 1:ncell(r)

# crop raster with a SpatialPolygon* object
if (require(rgdal) & require(rgeos)) {
p <- shapefile(system.file("external/lux.shp", package="raster"))
}

# Crop and mask
rc <- crop(r, extent(p))
rc <- mask(rc, p)


I wonder if there are new instructions to crop/mask rasters that are not yet 
documented...
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Tuesday, May 10, 2016 1:45 AM, Adam H Sparks  wrote:
Hi all,
I'm working on a vignette from the R-OpenSci auunconf,
https://github.com/saundersk1/auunconf16/blob/master/Vignette%20BoM.Rmd,
and am suddenly getting warning messages where I've never had any before
with this process.

In the "Australian Water Availability Project (AWAP)" section, at the very
end of the file, some grid files are downloaded and saved and GADM level 2
data are fetched to do cropping and masking.

Two weeks ago when we put this together it all worked with no error
messages. Today as I'm working on cleaning up the vignette I'm getting
error messages I've never seen before when doing these things.

For example, the line:
cfiles <- stack(dir(pattern = "grid$"))

Now returns these warning messages:
Warning messages:
1: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
2: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
3: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
4: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
5: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
6: In `[<-`(`*tmp*`, j, value = ) :
  implicit list embedding of S4 objects is deprecated
7: In `[<-`(`*tmp*`, nl, value = r) :
  implicit list embedding of S4 objects is deprecated

When I try to mask I get similar messages about the implicit list embedding
of S4 objects being deprecated.

For now in the vignette I've suppressed the warning messages as the code
appears to work, it just gives these warnings, but I'd like to understand
why I'm getting them.



-- 
Adam Sparks
Associate Professor Field Crops Pathology
Centre for Crop Health
University of Southern Queensland
Toowoomba, QLD, AU

@adamhsparks 

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Mask a map using statistical significance

[Thiago V. dos Santos via R-sig-Geo]

Ruben, 

This is a tentative reproducible example using "real" spatial data, rather than 
the synthetic data I proposed in the original question.

(I should have thought of this example earlier, in order to spatially please 
Barry lol)

require(raster)
require(rasterVis)

# Scratch raster objects
data(volcano)
r1 <- raster(volcano)

over <- ifelse(volcano >=160 & volcano <=180, 1, NA) # This is the "mask" raster
r2 <- raster(over)

# And this is the key step:
# To convert the "mask" raster to spatial points
r.mask <- rasterToPoints(r2, spatial=TRUE)

# Plot
levelplot(r1, margin=F) +
layer(sp.points(r.mask, pch=20, cex=0.3, alpha=0.8))


It looks a little better now. To control the parameters of the points (such as 
color, size, type etc), documentation is available in ?sp.points . 
 Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota


On Monday, May 9, 2016 4:59 PM, Rubén Fernández-Casal <rubenfca...@gmail.com> 
wrote:



I will appreciate instructions to reproduce this map...
Thanks...
Ruben FC

El 05/05/2016 00:06, "Thiago V. dos Santos via R-sig-Geo" 
<r-sig-geo@r-project.org> escribió: Thanks Barry and Arnaud (off-list) for the 
valuable hints. I ended up managing to use rasterVis' levelplot to produce the 
map I was looking for. The key step was to convert the raster I wanted to use 
as a mask to SpatialPoints, and then pass it as an additional layer to 
levelplot. The final map, produced using my actual data, can be visualized 
here: https://dl. dropboxusercontent.com/u/27700634/rainfall_trend.png. 
Stippling indicates regions exceeding the 95% statistical significance.

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Mask a map using statistical significance

[Thiago V. dos Santos via R-sig-Geo]

Thanks Barry and Arnaud (off-list) for the valuable hints.
I ended up managing to use rasterVis' levelplot to produce the map I was 
looking for.
The key step was to convert the raster I wanted to use as a mask to 
SpatialPoints, and then pass it as an additional layer to levelplot.
The final map, produced using my actual data, can be visualized here: 
https://dl.dropboxusercontent.com/u/27700634/rainfall_trend.png. Stippling 
indicates regions exceeding the 95% statistical significance.
 
|   |
|   |  |   |   |   |   |   |
|  |
|  |
| View on dl.dropboxuserconten... | Preview by Yahoo |
|  |
|   |


Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota


On Wednesday, May 4, 2016 11:02 AM, Barry Rowlingson 
<b.rowling...@lancaster.ac.uk> wrote:



Your example looks a bit rubbish because you have mostly isolated
pixels. Let's make an example with an area so we can see the shading.

> r.pvalue = raster(outer(-24:25,-24:25, function(a,b){a^2+b^2})/12)
> range(r.pvalue[])
[1] 0. 0.01041667

now convert your raster to polygons at the threshold you are
interested in. Let's imagine the area less than 0.001 (the central
circle here):

p_poly = rasterToPolygons(r.pvalue<0.001, dissolve=TRUE)

and overlay hatched polygons:

plot(r.pvalue) # or whatever your underlying data is
plot(p_poly[p_poly$layer==1,],density=5,add=TRUE, border=NA)

see help(polygon) for parameters to control the shading angle and
density. The "border=NA" above hides the outline which looks like
those maps you linked to.

Barry




On Wed, May 4, 2016 at 12:11 AM, Thiago V. dos Santos via R-sig-Geo
<r-sig-geo@r-project.org> wrote:
> Dear all,
>
> In climate studies, it is a common practice to perform some statistical test 
> between two fields (maps), and then plot the resulting map using a 
> significance mask. This masking is usually done by adding some kind of 
> pattern (shading, crosshatching etc) on top of the actual color palette.
>
> Examples can be seen here in this image 
> https://www3.epa.gov/climatechange/images/science/GlobalPrecipMap-large.png 
> and in the left images of this panel:
> http://www.nature.com/nclimate/journal/vaop/ncurrent/images_article/nclimate2996-f3.jpg
> In my case, I ran a statistical test for detecting trend on a time-series 
> raster and I now have one raster with the trend for rainfall (in degree C per 
> year) and one with the p-values associated to the test.
>
> My data looks roughly like this:
>
> require(raster)
>
> ## scratch raster objects and fill them with some random values
> r.trend <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
> r.trend[] <- round(runif(50 * 50, min=0, max=3), digits=2)
>
> r.pvalue <- raster(nrows=50, ncols=50, xmn=-58, xmx=-48, ymn=-33, ymx=-22)
> r.pvalue[] <- round(runif(50 * 50, min=0.0001, max=0.1), digits=5)
>
>
> What I would like to do is to plot r.trend, and on top of it plot r.pvalue 
> (as a pattern) where r.pvalues < 0.01 (corresponding to a significance level 
> of 99%).
>
> Has anyone here ever tried to do a similar plot and can point me to any 
> direction?
>
> I usually use rasterVis and ggplot2 to plot maps, but I would be open to try 
> some other package and even other SIG software other than R.
>
> Many thanks in advance, -- Thiago V. dos Santos
>
> PhD student
> Land and Atmospheric Science
> University of Minnesota
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Faster way to extract raster values using multiple polygons?

[Thiago V. dos Santos via R-sig-Geo]

Hi Ben,


Thank you so much for the suggestion. 

It is pretty fast, and I guess it's a good assumption for temperature. I'll 
have to do the same with precipitation, which is way more spatially variable 
than temperature. In that case, I think I'll have to wait the 10 hours it takes 
to run "extract" using a mean value of all cells covering a polygon.

It is a one-time processing anyway, so I should not worry too much about the 
waiting time.

Greetings,
-- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota




On Wednesday, April 6, 2016 10:48 AM, Ben Tupper <btup...@bigelow.org> wrote:
Hi,

The resolution of your raster data (1 degree) is much more coarse than what 
your polygons represent.  Could you short-circuit the process by assuming that 
the temp at the centroid of each polygon would suitably represent the mean 
temperature across each polygon?  Unless you have some much bigger polygons, I 
can't imagine it will be very far off. If so, then you could pretty quickly 
extract the values for each layer in the raster at each centroid.  Perhaps like 
this?

cents <- coordinates(br_sub)
v <- extract(b, cents)

Is that close enough?

Cheers,
Ben


> On Apr 6, 2016, at 8:00 AM, Thiago V. dos Santos via R-sig-Geo 
> <r-sig-geo@r-project.org> wrote:
> 
> Dear all,
> 
> I am trying to extract a time series from a raster brick using multiple 
> polygons.
> 
> The raster brick is a global temperature netcdf file with 1200 layers 
> (months) and the polygons are each of the 5567 municipalities in Brazil. I 
> intend to extract the temperature time series for each municipality.
> 
> As a final result, I would like to have a data frame containing: -the name 
> and coordinates of the municipality, -the date tag from the raster, and -the 
> extracted temperature value.
> 
> My initial, VERY SLOW attempt was something like this:
> 
> 
> library(raster)
> 
> # Create some sample raster data
> idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month")
> r <- raster(ncol=360, nrow=180)
> b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
> b <- setZ(b, idx)
> 
> # Import polygon data
> br <- getData("GADM", country="Brazil", level=2) # about 7MB to download
> 
> # Subset the SMALLEST state - only 75 municipalities
> br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),]
> plot(br_sub)
> 
> # Now let's extract the data for each municipality
> beginCluster()
> e <- extract(b, br_sub, fun=mean, df=T)
> endCluster()
> 
> 
> Even using the smallest state possible, this example takes about 20 minutes 
> to run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, 
> there are states with over 850 municipalities. And remember, in total there 
> are over 5500 municipalities I need to extract the data for.
> 
> I have the feeling that my code is not very optimized.
> 
> Has anybody ever dealt with this amount of data in this kind of raster 
> operation? Is there any fastest way to achieve my expected result?
> Thanks in advance,
> -- Thiago V. dos Santos
> 
> PhD student
> Land and Atmospheric Science
> 
> University of Minnesota
> 
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

[R-sig-Geo] Faster way to extract raster values using multiple polygons?

[Thiago V. dos Santos via R-sig-Geo]

Dear all,

I am trying to extract a time series from a raster brick using multiple 
polygons.

The raster brick is a global temperature netcdf file with 1200 layers (months) 
and the polygons are each of the 5567 municipalities in Brazil. I intend to 
extract the temperature time series for each municipality.

As a final result, I would like to have a data frame containing: -the name and 
coordinates of the municipality, -the date tag from the raster, and -the 
extracted temperature value.

My initial, VERY SLOW attempt was something like this:


library(raster)

# Create some sample raster data
idx <- seq(as.Date("1960/1/1"), as.Date("1990/12/01"), by = "month")
r <- raster(ncol=360, nrow=180)
b <- brick(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r)
b <- setZ(b, idx)

# Import polygon data
br <- getData("GADM", country="Brazil", level=2) # about 7MB to download

# Subset the SMALLEST state - only 75 municipalities
br_sub <- br[br$NAME_1=="Sergipe" & !is.na(br$NAME_1),]
plot(br_sub)

# Now let's extract the data for each municipality
beginCluster()
e <- extract(b, br_sub, fun=mean, df=T)
endCluster()


Even using the smallest state possible, this example takes about 20 minutes to 
run on a dual-core 2.5GHz Macbook Pro using four threads. As a reference, there 
are states with over 850 municipalities. And remember, in total there are over 
5500 municipalities I need to extract the data for.

I have the feeling that my code is not very optimized.

Has anybody ever dealt with this amount of data in this kind of raster 
operation? Is there any fastest way to achieve my expected result?
Thanks in advance,
-- Thiago V. dos Santos

PhD student
Land and Atmospheric Science

University of Minnesota

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo