Re: [R-sig-Geo] rgdal 1.5-10 published on CRAN
Thank you very much for all the great work you and all the developers constantly do. Mauricio Mauricio Zambrano-Bigiarini, PhD = Department of Civil Engineering Faculty of Engineering and Sciences Universidad de La Frontera, Temuco, Chile http://hzambran.github.io/ https://orcid.org/-0002-9536-643X = mailto : mauricio.zambr...@ufrontera.cl work-phone : +56 45 259 2812 = "What makes the desert beautiful is that somewhere it hides a well" (Antoine de Saint-Exupery = 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
[R-sig-Geo] How to correctly compare the CRS of 'sf' and 'raster' objects?
Dear mailing list, I'm wondering what is the correct way of comparing the CRS of a vectorial (simple feature) and a raster object. Please look a the following reproducible example: --- CODE START library(sf) library(raster) # 1) Getting the character corresponding to lat/lon CRS string latlon.crs <- st_crs(4326)$proj4string # 2) Creating a simple feature point geometry and assigning lat/lon CRS vec <- st_sfc(st_point(c(0,0)), st_point(c(1,1))) st_crs(vec) <- 4326 # 3) Creating a simple feature point geometry and assigning lat/lon CRS rast <- raster(nrows=108, ncols=21, xmn=0, xmx=10) crs(rast) <- latlon.crs # 4) Getting the CRS of the point and raster objects vect.crs <- sf::st_crs(vec)$proj4string rast.crs <- sf::st_crs(rast)$proj4string # 5) Making the CRS comparison are.equal <- vect.crs == rast.crs --- CODE END The previous comparison works fine in R 4.0.0 (raster_3.1-5 sp_1.4-1, sf_0.9-2) with the following set of libraries: GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2 but the same line of code fails with the newer set of libraries: GEOS 3.8.1, GDAL 3.0.4, PROJ 7.0.0 Is there any "correct" way of comparing the CRS of the previous vectorial and raster object in 'sf', something similar to sp::identicalCRS? If so, what is the minimum required version of the GEOS, GDAL, and PROJ libraries? Thanks in advance for any help, Mauricio Zambrano-Bigiarini, PhD = Department of Civil Engineering Faculty of Engineering and Sciences Universidad de La Frontera, Temuco, Chile http://hzambran.github.io/ = mailto : mauricio.zambr...@ufrontera.cl work-phone : +56 45 259 2812 = "The world is changed by your example, not your opinion". (Paulo Coelho) = 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
Re: [R-sig-Geo] Fwd: Spatiotemporal analysis r command
https://www.r-project.org/posting-guide.html IHTH Mauricio Zambrano-Bigiarini, PhD = Department of Civil Engineering Faculty of Engineering and Sciences Universidad de La Frontera, Temuco, Chile http://hzambran.github.io/ = mailto : mauricio.zambr...@ufrontera.cl work-phone : +56 45 259 2812 = "Focus is about saying No" (Steve Jobs) = Linux user #454569 -- Linux Mint user On 1 July 2018 at 06:05, Yalemzewod Gelaw wrote: > Hello everyone, > I am planning to do a spatiotemporal distribution of hiv/tb coinfection and > the effects of socio-demographic variables. I have a shapefile for 128 > districts. For this ecological area analysis I am planning to use > Bayesian Poisson regression analysis. > I here with looking for your help tutorial/notes/ r commands for the > spatiotemporal analysis and Bayesian Poisson regression. > > Thank you for any help. > > *Regards, * > > Yalem > -- > Best, Yalemzewod Assefa Gelaw (Yalem) > > [[alternative HTML version deleted]] > > ___ > R-sig-Geo mailing list > R-sig-Geo@r-project.org > https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- 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
Re: [R-sig-Geo] Fwd: Re: Adding text to rasterVis levelplot
Thank you very much you both Thiago and Florian for your help.
Just for future searches on Google, the minimal example that allows
creating the plot i was looking for was:
-- START -
library(rasterVis)
f <- system.file("external/test.grd", package="raster")
r <- raster(f)
s <- stack(r, r+500, r-500, r+200)
col.titles = c('col1','col2','','')
row.titles = c('row1','row2')
levelplot(s, layout=c(2,2),
names.attr=col.titles,
ylab=row.titles,
scales = list(y = list(rot = 90) )
)
-- END -
However, the example provided by Florian opened up new possibilities
with the grid package for more customisation.
Kind regards,
Mauricio
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
On 21 July 2017 at 04:21, Florian Detsch
<florian.det...@staff.uni-marburg.de> wrote:
> Ah, I forgot to include the definition of text vectors which must be
> inserted before the for() loop - sorry for that.
>
> ---
> ## labels
> cls <- c("col1", "col2")
> rws <- c("row1", "row2")
> ---
>
>
> Forwarded Message
> Subject:Re: [R-sig-Geo] Adding text to rasterVis levelplot
> Date: Fri, 21 Jul 2017 10:05:08 +0200
> From: Florian Detsch <florian.det...@staff.uni-marburg.de>
> To: r-sig-geo@r-project.org
>
>
>
> Hi Mauricio,
>
> Complementing Thiago's answer using 'names.att', here is a more manual
> approach which might come in handy when a finer control over text
> annotations is required. It is built upon lattice::trellis.focus() which
> lets you select relevant panels of your trellis graph in an iterative
> manner.
>
> ---
> p <- levelplot(s, layout = c(2, 2), names.att = rep("", 4),
> scales = list(y = list(rot = 90)))
>
> library(grid)
>
> png("~/rasterVis.png", width = 14, height = 16, units = "cm", res = 300L)
> grid.newpage()
> print(p, newpage = FALSE)
>
> ## loop over panels to be labelled (ie 1:3)
> panels <- trellis.currentLayout()
> for (i in 1:3) {
>
># focus on current panel of interest and disable clipping
>ids <- which(panels == i, arr.ind = TRUE)
>trellis.focus("panel", ids[2], ids[1], clip.off = TRUE)
>
># add labels
>if (i %in% c(1, 3)) {
> if (i == 1) {
>grid.text(cls[1], x = .5, y = 1.1) # add 'col1'
>grid.text(rws[1], x = -.35, y = .5, rot = 90) # add 'row1'
> } else {
>grid.text(rws[2], x = -.35, y = .5, rot = 90) # add 'row2'
> }
>} else {
> grid.text(cls[2], x = .5, y = 1.1) # add 'col2'
>}
>
>trellis.unfocus()
> }
>
> dev.off()
> ---
>
> A very similar use case, from which the above approach definitely took
> some inspiration, can be found here:
> https://stackoverflow.com/questions/7649737/is-it-possible-to-update-a-lattice-panel-in-r.
> Have a look at ?grid.text() and ?gpar() to learn more about
> customization options.
>
> Cheers,
> Florian
>
> --
> Dr. Florian Detsch
> Environmental Informatics
> Department of Geography
> Philipps-Universität Marburg
> Deutschhausstraße 12
> 35032 (parcel post: 35037) Marburg, Germany
>
> Phone: +49 (0) 6421 28-25323
> Web:
> http://www.uni-marburg.de/fb19/fachgebiete/umweltinformatik/detschf/index.html
>
>
> On 20.07.2017 05:07, Thiago V. dos Santos via R-sig-Geo wrote:
>> 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
>> <mauricio.zambr...@ufrontera.cl> 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
[R-sig-Geo] Adding text to rasterVis levelplot
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
[R-sig-Geo] raster: crop error ('nrows > 0 is not TRUE')
Dear list,
I'm applying the 'crop' command of the raster pacakge to subset a
RasterStack ('r') with a SpatialPolygonsDataFrame ('boundary'):
raster.crop <- crop(r, boundary, snap="out")
and what I get is:
Error: nrows > 0 is not TRUE
The summary of the two previous objects are:
> r
class : RasterStack
dimensions : 666, 211, 140526, 8 (nrow, ncol, ncell, nlayers)
resolution : 0.05, 0.05 (x, y)
extent : -76.3, -65.75, -50, -16.7 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
names : Year.2003, Year.2004, Year.2005, Year.2006, Year.2007,
Year.2008, Year.2009, Year.2010
min values : 0, 0, 0, 0, 0,
0, 0, 0
max values :844.79,706.53,616.73,888.75,701.53,
853.47,642.07,986.31
> boundary
class : SpatialPolygonsDataFrame
features: 16
extent : -75.693, -66.419, -55.914, -17.498 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0
variables : 1
names : NOM_REG
min values : 1
max values : 9
the previous 'crop' command worked perfectly with a different
RasterStack object, with a similar spatial extent, and I don't now
what is wrong here.
The output of my sessionInfo() is:
sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.2 LTS
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_GB.UTF-8
LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] colorspace_1.2-6reshape2_1.4.1 rasterVis_0.40
latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-34
rgdal_1.1-10hydroGOF_0.3-8-7hydroTSM_0.4-8 xts_0.9-7
zoo_1.7-13
[12] raster_2.5-8sp_1.2-3
loaded via a namespace (and not attached):
[1] Rcpp_0.12.6 magrittr_1.5 maptools_0.8-39
viridisLite_0.1.3 FNN_1.1 stringr_1.0.0 plyr_1.8.4
tools_3.3.1 parallel_3.3.1grid_3.3.1gstat_1.1-3
e1071_1.6-7
[13] class_7.3-14 intervals_0.15.1 automap_1.0-14ncdf4_1.15
stringi_1.1.1 spacetime_1.1-5 reshape_0.8.5
foreign_0.8-66hexbin_1.27.1
Any idea about how to solve this issue is highly appreciated.
Thanks in advance,
Mauricio Zambrano-Bigiarini, PhD
=
Department of Civil Engineering
Faculty of Engineering and Sciences
Universidad de La Frontera
PO Box 54-D, Temuco, Chile
http://hzambran.github.io/
=
mailto : mauricio.zambr...@ufrontera.cl
work-phone : +56 45 259 2812
=
"To accomplish great things, we must not only act,
but also dream; not only plan, but also believe"
(Anatole France)
=
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
Re: [R-sig-Geo] raster pkg: overwriting the dimensions detected by netcdf4 (PERSIANN-CDR)
On 3 February 2016 at 11:34, Michael Sumner <mdsum...@gmail.com> wrote:
>
>
> On Thu, 4 Feb 2016 at 01:26 MAURICIO ZAMBRANO BIGIARINI
> <mauricio.zambr...@ufrontera.cl> wrote:
>>
>> Dear spatial community,
>>
>> While reading some netCDFfiles with the raster package , I got a
>> "rotated" file, where the columns where read as rows and vice versa:
>>
>> - START
>> x <- raster("PERSIANN-CDR_v01r01_20090101_c20140523.nc")
>> Loading required namespace: ncdf4
>> plot(x)
>> x
>> class : RasterLayer
>> dimensions : 1440, 480, 691200 (nrow, ncol, ncell)
>> resolution : 0.25, 0.25 (x, y)
>> extent : -60, 60, 0, 360 (xmin, xmax, ymin, ymax)
>> coord. ref. : NA
>> data source : /home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
>> names : NOAA.Climate.Data.Record.of.PERSIANN.CDR.daily.precipitation
>> z-value : 2009-01-01
>> zvar: precipitation
>> - END
>>
>
> Transpose might be sufficient?
>
> tx <- t(x)
>
> plot(tx)
> tx
Thanks Michale,
In fact, that is what I'm doing for reading the files, and the spatial
extent is correct and the precipitation values also seems to be right.
However, my question was more "conceptual", in the sense that the
provider says that the original file is NOT rotated, while the only
way I have to get the correct files is rotating the original file
with the 't' command you mentioned.
In addition, while reading the file with the latest version of GRASS
GIS (7.0.3), and I got the following error message:
"
(Tue Feb 2 19:41:29 2016)
r.in.gdal input=/home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
output=PERSIANN_CDR_v01r01_20090101_c20140523 -e
Warning 1: dimension #2 (lat) is not a Longitude/X
dimension.
Warning 1: dimension #1 (lon) is not a Latitude/Y dimension.
ERROR: Input raster map is flipped or rotated - cannot import. You may
use 'gdalwarp' to transform the map to North-up.
"
So, I'm quite sure that the original file IS rotated, but I didn't
have more technical arguments to give to the data provider to
demonstrate that the file is rotated, because they insist that the
rotation is because R is not reading the file in the proper way
mzb
=
"In the end, it's not the years in your life that count.
It's the life in your years". (Abraham Lincoln)
=
Linux user #454569 -- Linux Mint user
>> 60W-60E and 0-360N, while the correct extent should be: 60S, 60N, 0E,
>> 360E.
>>
>> The file can be downloaded from:
>>
>> ftp://data.ncdc.noaa.gov/cdr/persiann/files/2009/PERSIANN-CDR_v01r01_20090101_c20140523.nc
>>
>>
>> After contacting the providers of the file, they mentioned that when
>> reading the file I have to provide the dimension information, i.e.,
>> 480 rows and 1440 columns, while the raster package detects the other
>> way around.
>>
>> I would highly appreciate if you could tell me:
>>
>> 1) is the previous situation a bug of my netCDF driver (netCDF 4.1.3,
>> GDAL 1.11.2, released 2015/02/10) or a problem of the dimensions
>> actually recorded in the netCDF file ?
>>
>> 2) is there any way to overwrite the dimensions detected by netcdf4
>> when reading the file with the 'raster' command ?
>>
>>
>> My sessionInfo():
>> R version 3.2.3 (2015-12-10)
>> Platform: x86_64-pc-linux-gnu (64-bit)
>> Running under: Ubuntu 14.04.2 LTS
>>
>> locale:
>> [1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
>> [3] LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
>> [5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_GB.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
>>
>> other attached packages:
>> [1] raster_2.5-2 sp_1.2-1
>>
>> loaded via a namespace (and not attached):
>> [1] tools_3.2.3 Rcpp_0.12.3 ncdf4_1.15 grid_3.2.3
>> [5] lattice_0.20-33
>>
>>
>> Thanks in advance and kind regards
>>
>> Mauricio Zambrano-Bigiarini, PhD
>>
>> =
>> Dept. of Civil Engineering
>> Faculty of Engineering and Sciences
>> Universidad de La Frontera
>> PO Box 54-D, Temuco, Chile
>> http://ingenieriacivil.ufro.cl/
>> =
>> mailto : mauricio.zambr...@ufrontera.cl
&
[R-sig-Geo] raster pkg: overwriting the dimensions detected by netcdf4 (PERSIANN-CDR)
Dear spatial community,
While reading some netCDFfiles with the raster package , I got a
"rotated" file, where the columns where read as rows and vice versa:
- START
x <- raster("PERSIANN-CDR_v01r01_20090101_c20140523.nc")
Loading required namespace: ncdf4
plot(x)
x
class : RasterLayer
dimensions : 1440, 480, 691200 (nrow, ncol, ncell)
resolution : 0.25, 0.25 (x, y)
extent : -60, 60, 0, 360 (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
names : NOAA.Climate.Data.Record.of.PERSIANN.CDR.daily.precipitation
z-value : 2009-01-01
zvar: precipitation
- END
The previous summary shows that the spatial extent of the file is
60W-60E and 0-360N, while the correct extent should be: 60S, 60N, 0E,
360E.
The file can be downloaded from:
ftp://data.ncdc.noaa.gov/cdr/persiann/files/2009/PERSIANN-CDR_v01r01_20090101_c20140523.nc
After contacting the providers of the file, they mentioned that when
reading the file I have to provide the dimension information, i.e.,
480 rows and 1440 columns, while the raster package detects the other
way around.
I would highly appreciate if you could tell me:
1) is the previous situation a bug of my netCDF driver (netCDF 4.1.3,
GDAL 1.11.2, released 2015/02/10) or a problem of the dimensions
actually recorded in the netCDF file ?
2) is there any way to overwrite the dimensions detected by netcdf4
when reading the file with the 'raster' command ?
My sessionInfo():
R version 3.2.3 (2015-12-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.2 LTS
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8LC_COLLATE=en_GB.UTF-8
[5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_GB.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
other attached packages:
[1] raster_2.5-2 sp_1.2-1
loaded via a namespace (and not attached):
[1] tools_3.2.3 Rcpp_0.12.3 ncdf4_1.15 grid_3.2.3
[5] lattice_0.20-33
Thanks in advance and kind regards
Mauricio Zambrano-Bigiarini, PhD
=
Dept. of Civil Engineering
Faculty of Engineering and Sciences
Universidad de La Frontera
PO Box 54-D, Temuco, Chile
http://ingenieriacivil.ufro.cl/
=
mailto : mauricio.zambr...@ufrontera.cl
work-phone : +56 45 259 2812
=
"In the end, it's not the years in your life that count.
It's the life in your years". (Abraham Lincoln)
=
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
Re: [R-sig-Geo] raster pkg: overwriting the dimensions detected by netcdf4 (PERSIANN-CDR)
On 3 February 2016 at 12:56, Michael Sumner <mdsum...@gmail.com> wrote:
>
>
> On Thu, 4 Feb 2016 at 02:20 MAURICIO ZAMBRANO BIGIARINI
> <mauricio.zambr...@ufrontera.cl> wrote:
>>
>> On 3 February 2016 at 11:34, Michael Sumner <mdsum...@gmail.com> wrote:
>> >
>> >
>> > On Thu, 4 Feb 2016 at 01:26 MAURICIO ZAMBRANO BIGIARINI
>> > <mauricio.zambr...@ufrontera.cl> wrote:
>> >>
>> >> Dear spatial community,
>> >>
>> >> While reading some netCDFfiles with the raster package , I got a
>> >> "rotated" file, where the columns where read as rows and vice versa:
>> >>
>> >> - START
>> >> x <- raster("PERSIANN-CDR_v01r01_20090101_c20140523.nc")
>> >> Loading required namespace: ncdf4
>> >> plot(x)
>> >> x
>> >> class : RasterLayer
>> >> dimensions : 1440, 480, 691200 (nrow, ncol, ncell)
>> >> resolution : 0.25, 0.25 (x, y)
>> >> extent : -60, 60, 0, 360 (xmin, xmax, ymin, ymax)
>> >> coord. ref. : NA
>> >> data source : /home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
>> >> names :
>> >> NOAA.Climate.Data.Record.of.PERSIANN.CDR.daily.precipitation
>> >> z-value : 2009-01-01
>> >> zvar: precipitation
>> >> - END
>> >>
>> >
>> > Transpose might be sufficient?
>> >
>> > tx <- t(x)
>> >
>> > plot(tx)
>> > tx
>>
>>
>> Thanks Michale,
>>
>> In fact, that is what I'm doing for reading the files, and the spatial
>> extent is correct and the precipitation values also seems to be right.
>>
>> However, my question was more "conceptual", in the sense that the
>> provider says that the original file is NOT rotated, while the only
>> way I have to get the correct files is rotating the original file
>> with the 't' command you mentioned.
>>
>
>> In addition, while reading the file with the latest version of GRASS
>> GIS (7.0.3), and I got the following error message:
>>
>> "
>> (Tue Feb 2 19:41:29 2016)
>> r.in.gdal input=/home/hzambran/PERSIANN-CDR_v01r01_20090101_c20140523.nc
>> output=PERSIANN_CDR_v01r01_20090101_c20140523 -e
>> Warning 1: dimension #2 (lat) is not a Longitude/X
>> dimension.
>> Warning 1: dimension #1 (lon) is not a Latitude/Y dimension.
>> ERROR: Input raster map is flipped or rotated - cannot import. You may
>> use 'gdalwarp' to transform the map to North-up.
>> "
>>
>>
>> So, I'm quite sure that the original file IS rotated, but I didn't
>> have more technical arguments to give to the data provider to
>> demonstrate that the file is rotated, because they insist that the
>> rotation is because R is not reading the file in the proper way
>>
>
>
> I wouldn't bother with the conversation if I were you. Just get your toolkit
> in order and make sure you know what it's doing and why. If you explore
> their tools you might quickly understand why they're so adamant, but it's
> not worth the time IMO.
>
> library(fortunes)
> fortune("illogical")
>
> If there's a predictable pattern in the file you might explore a patch for
> raster to auto-detect and handle it? I'd be happy to do that - if you
> provide the data and as much information about it and its kind as you can
> find. If you want to know more I would start with ncdf4 directly, so you
> see what raster has to do down there.
Thanks your very much Michael and Chris for your comments.
In fact, reading the file with gdal_translate I got a warning message
about the dimensions:
Checking gdal_installation...
Scanning for GDAL installations...
Checking Sys.which...
GDAL version 1.11.2
GDAL command being used: "/usr/bin/gdal_translate" -of "GTiff" -b "1"
"NETCDF:PERSIANN-CDR_v01r01_20090101_c20140523.nc:precipitation"
"PERSIANN-CDR_v01r01_20090101_c20140523.nc.tif"
Warning 1: dimension #2 (lat) is not a Longitude/X dimension.
Warning 1: dimension #1 (lon) is not a Latitude/Y dimension.
also, 'ncview' indicates that the Y variable corresponds to longitude,
going from 0.125 to 359.875, while the X variable corresponds to
latitude, with values going from 59.75 to -59.875.
However, 'ncdump -h ' shows that the attributes are correctly stored
within the file (but with no extent):
netcdf PERSIANN-CDR_v01r01_20090101_c20140523 {
dimensions:
time = 1 ;
lat = 480 ;
lon = 1440 ;
nv = 2 ;
variables:
int time(tim
Re: [R-sig-Geo] mapView: basic interactive viewing of spatial data in R
On 24 July 2015 at 09:26, Edzer Pebesma edzer.pebe...@uni-muenster.de wrote: Hi Tim, thanks for starting the discussion. Thanks Tim for this nice implementation and to all of you for this very interesting discussion. I started a similar discussion (off-list, with package maintainers) for the case where google map or openstreetmap maps are used as background in map plots created in R, support for which is now scattered over ggmap::get_map, RgoogleMaps::GetMap, dismo::gmap and https://github.com/hadley/rastermap. (on r-forge, sp::plot and sp::spplot can now deal with these, see sp::demo(webmap)). by the way Edzer, the sp::demo(webmap) din't work for me with sp_1.1-1: library(sp) demo(meuse)# working demo(webmap) # not working could you tell me what I'm doing wrong ? Thanks in advance, Mauricio Zambrano-Bigiarini, PhD = Dept. of Civil Engineering Faculty of Engineering and Sciences Universidad de La Frontera PO Box 54-D, Temuco, Chile = mailto : mauricio.zambr...@ufrontera.cl work-phone : +56 45 259 2812 http://ingenieriacivil.ufro.cl/ = When the pupil is ready, the master arrives. (Zen proverb) = Linux user #454569 -- Linux Mint use ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] raster/rgdal- problem: Too many open files (Linux)
On 08/05/2013 04:37 PM, Jon Olav Skoien wrote:
Dear list,
We have a problem which appears to be a bug in either rgdal or raster,
although it could also be a bug in base R or in our understanding of how
to deal with connections.
We have a process which is writing a rather large (~10-20.000) number of
geoTiffs via writeRaster. However, the process has frequently stopped
with an error of the type:
Error in .local(.Object, ...) :
TIFFOpen:/local0/skoiejo/hri/test.tif: Too many open files
The issue seems to be the creation of temp-files in the temp directory
which is given by tempdir(), not by raster:::.tmpdir(). These temp-files
seem to be created by the call
transient - new(GDALTransientDataset, driver=driver, rows=r@nrows,
cols=r@ncols, bands=nbands, type=dataformat, fname=filename,
options=options, handle=NULL)
from raster:::.getGDALtransient
The temp-files are deleted after writing the geoTiff, but are not
removed from the list of open files in Linux, which on our system was
limited to 1024 files (ulimit -n) per process. Below is a script which
can replicate the issue (takes a few minutes to reach 1024) and
sessionInfo().
Currently we are trying to solve the issue by increasing the limit of
file connections, but we would prefer a solution where the connections
are properly deleted, either before writeRaster finishes, or a command
which we can include in our script, either R-code or a call to System().
The connections are not visible via showConnections(), and
closeAllConnections() does not help.
Thanks,
Jon
I stumbled across the same problem (with exactly the same configuration
reported by Jon with 'sessionInfo()'), while trying to change the values
of some pixels in more than 6000 maps.
Thank you very much Jon for the detailed report about the problem,
which helped me to find a workaround to this problem (so far, just to
split the 6000 maps in smaller groups).
r - raster(system.file(external/test.grd, package=raster))
for (ifile in 1:2000) {
writeRaster(r, test.tif, format = GTiff, overwrite = TRUE)
print(ifile)
}
After trying the previous reproducible code, I don't understand why I
got the error when ifile=1019 and not 1024:
[1] 1018
[1] 1019
Error in .local(.Object, ...) :
TIFFOpen:/home/hzambran/test.tif: Too many open files
Thanks again Jon for sharing your findings about this.
All the best,
Mauricio Zambrano-Bigiarini, Ph.D
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission
=
Sometimes life's going to hit you in the head
with a brick. Don't lose faith (Steve Jobs)
After the script stops, I checked the open files from the process, and
got the following:
lsof -aPn -p 596 | more
COMMAND PIDUSER FD TYPE DEVICE SIZE/OFF NODE NAME
R 596 skoiejo cwdDIR 8,17 139264 4213416
/local0/skoiejo/Raster_temp
R 596 skoiejo rtdDIR 253,0 40962 /
R 596 skoiejo txtREG 253,014192 2500779
/usr/lib64/R/bin/exec/R
R 596 skoiejo memREG 253,0 2679280 2501177
/usr/lib64/R/lib/libR.so
..
R 596 skoiejo 1015u REG 253,0 238 4983213
/tmp/RtmpNyOKCk/toptest.tif (deleted)
R 596 skoiejo 1016u REG 253,0 238 4983214
/tmp/RtmpNyOKCk/qxdtest.tif (deleted)
R 596 skoiejo 1017u REG 253,0 238 4983215
/tmp/RtmpNyOKCk/zwotest.tif (deleted)
R 596 skoiejo 1018u REG 253,0 238 4983216
/tmp/RtmpNyOKCk/cnqtest.tif (deleted)
R 596 skoiejo 1019u REG 253,0 238 4983217
/tmp/RtmpNyOKCk/lottest.tif (deleted)
R 596 skoiejo 1020u REG 253,0 238 4983218
/tmp/RtmpNyOKCk/fartest.tif (deleted)
R 596 skoiejo 1021u REG 253,0 238 4983219
/tmp/RtmpNyOKCk/vsqtest.tif (deleted)
R 596 skoiejo 1022u REG 253,0 238 4983220
/tmp/RtmpNyOKCk/czptest.tif (deleted)
Even if tested by someone with a limit higher than 2000, it should still
be possible to see the long list of open connections, as above.
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-redhat-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=C 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
other attached packages:
[1] rgdal_0.8-10 raster_2.1-49 sp_1.0-11
loaded
Re: [R-sig-Geo] working with large ascii files in R
On 07/18/2013 09:26 AM, Roman Luštrik wrote:
Have you considered updating your R and associated packages?
As Roman mentioned, I also suggest you, as first step, to update R (and
all the packages) up to version 3 or higher, because one of the
significant user changes introduced in R 3.0.0 was related to memory
(from the NEWS file):
It is now possible for 64-bit builds to allocate amounts of
memory limited only by the OS. It may be wise to use OS
facilities (e.g. ulimit in a bash shell, limit in csh), to set
limits on overall memory consumption of an R process,
particularly in a multi-user environment. A number of packages
need a limit of at least 4GB of virtual memory to load.
64-bit Windows builds of R are by default limited in memory usage
to the amount of RAM installed: this limit can be changed by
command-line option --max-mem-size or setting environment
variable R_MAX_MEM_SIZE.
Kind regards,
Mauricio Zambrano-Bigiarini, Ph.D
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
Work Phone : +39 0332 789588
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission
=
The journey is the reward (Steve Jobs)
There's also a plethora of other GUI based GIS systems. Just the other day,
I found Quantum GIS to be very user friendly.
That being said, I suggest you use R. :)
Cheers,
Roman
On Thu, Jul 18, 2013 at 9:05 AM, Rocio Ponce r.ponce.re...@gmail.comwrote:
Hi there,
I'm a newie in using R to make maps... I used to work with ArcGis but I ran
out of a valid license...
I'm trying to crop the future Worldclim datasets layers which are huge
ascii files (of the whole world at 1km resolution, so they are between
3-4GB each) to the extent of Mexico using the library (raster).
I used the following code:
setwd('G:/future_A2a/2020/bccr_bcm2_0_sres_a2_2020s/')
maask- read.asc(C:/xalapa/worldclim_pres/bio01.asc,gz=FALSE)
mask-raster.from.asc(maask,projs=NA)
files - list.files(pattern='\\.asc$')
s - stack(files)
crop(s, mask)
Error: cannot allocate vector of size 253.2 Mb
I'm using R.2.14.1 on a computer running on Windows XP
Any help that you could provide me would be greatly appreciated,
Cheers,
Rocio
--
Dr. Rocio Ponce-Reyes
http://rocioponcereyes.wordpress.com/
(disculpe la falta de acentos en este correo electronico)
[[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] writeRaster to ascii file (.asc)
On 10/07/13 18:38, Manuel Spínola wrote: Dear list members, I am trying to export a raster file to a ascii file (.asc) but the min and max values are not the same. Any idea why? Hi Manuel, Could you provide your sessionInfo() ? bio1res class : RasterLayer dimensions : 382, 407, 155474 (nrow, ncol, ncell) resolution : 0.00833, 0.00833 (x, y) extent : -85.95, -82.55833, 8.041667, 11.225 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : bio1_23 values : 47, 272 (min, max) writeRaster(bio1res, filename=temp.asc, format=ascii,overwrite=TRUE) class : RasterLayer dimensions : 382, 407, 155474 (nrow, ncol, ncell) resolution : 0.00833, 0.00833 (x, y) extent : -85.95, -82.55833, 8.041667, 11.225 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : /Users/manuelspinola/Dropbox/DistribucioÌ�n de especies/World_Clim_Costa_Rica/temp.asc names : temp values : -2147483648, 2147483647 (min, max) Unfortunately, your example is not reproducible, so I tried with a different example and the values in the original and the written ascii file were the same (rgdal_0.8-10, raster_2.1-48, sp_1.0-11): -- START -- library(raster) Loading required package: sp # Creation of a raster file for testing r - raster(system.file(external/test.grd, package=raster)) r #class : RasterLayer #dimensions : 115, 80, 9200 (nrow, ncol, ncell) #resolution : 40, 40 (x, y) #extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax) #coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.156160 +lon_0=5.387639 +k=0.079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs #data source : /usr/lib64/R/library/raster/external/test.grd #names : test #values : 128.434, 1805.78 (min, max) writeRaster(r, filename=temp.asc, format=ascii,overwrite=TRUE) #rgdal: version: 0.8-10, (SVN revision 478) #Geospatial Data Abstraction Library extensions to R successfully loaded #Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08 #Path to GDAL shared files: /usr/share/gdal #Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480] #Path to PROJ.4 shared files: (autodetected) #class : RasterLayer #dimensions : 115, 80, 9200 (nrow, ncol, ncell) #resolution : 40, 40 (x, y) #extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax) #coord. ref. : NA #data source : /home/zambrhe/temp.asc #names : temp # reading the .asc file form disk r2 - raster(temp.asc) r2 #class : RasterLayer #dimensions : 115, 80, 9200 (nrow, ncol, ncell) #resolution : 40, 40 (x, y) #extent : 178400, 181600, 329400, 334000 (xmin, xmax, ymin, ymax) #coord. ref. : NA #data source : /home/zambrhe/temp.asc #names : temp # Summary of the ascii file in the hard disk #summary(r2) # temp #Min. 128.4340 #1st Qu. 293.2325 #Median 371.4120 #3rd Qu. 499.8195 #Max.1805.7800 #NA's6097. # Summary of the original raster file summary(r) # test #Min. 128.4340 #1st Qu. 293.2325 #Median 371.4120 #3rd Qu. 499.8195 #Max.1805.7800 #NA's6097. -- END -- Kind regards, Mauricio Zambrano-Bigiarini, Ph.D -- = Water Resources Unit Institute for Environment and Sustainability (IES) Joint Research Centre (JRC), European Commission TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT Work Phone : +39 0332 789588 webinfo: http://floods.jrc.ec.europa.eu/ = DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as sta- ting an official position of the European Commission = The journey is the reward (Steve Jobs) ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] R Package for Climatology particular for netCDF file
On 26/03/13 18:13, ping yang wrote:
Dear Dr. Pierce and R specialists,
I am much appreciated your work on ncdf(ncdf4 and RnetCDF) packages,
however, I am wondering if those package can go further.
Is there a package for calculating climate statistics ,e.g. monthly mean,
seasonal mean, yearly mean,spatial mean( like gridboxmean in CDO), and
climatology for multiyear analysis in R or CRAN ? and Is there a plotting
library for the climate variables particularly for NetCDF.
I found a little trivial when I plan to do statistical based on daily data
for long history (I do it using some tools but not efficient, everytime I
need to generate a file and then use ncdf funcitons to read and then plot.
Thanks,
Ping
Hi Ping,
Regarding the climate statistics, you may have a look to the hydroTSM
package, in particular to the functions:
?monthlyfunction
?seasonalfunction
?dm2seasonal
and the vignette in:
http://cran.r-project.org/web/packages/hydroTSM/vignettes/hydroTSM_Vignette.pdf
IHTH,
Kinds,
Mauricio Zambrano-Bigiarini, Ph.D
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:9}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] raster directions to vectorial lines (Stream to Feature)
On 22/03/13 05:31, Gavan McGrath wrote:
Hi Mauricio,
Below is a simple pair of functions for plotting a flow network over an image
which should achieve what you were after. You may want to play around with how
your flow direction numbers 1 -9 match up with those in the lookupfidr
function. Its a bit cluncky but gets the job done.
The function fdirPlot takes as input a matrix imData the background image, and
a matrix fdirs of flow direction integers.
Good luck.
Thank you very much Gavan for you valuable help, and apologise for my
late reply.
With slight modifications to the code you gave me, I managed to do what
I want on the screen device:
- START reproducible example
lookupfdir -function(fdir) {
switch(paste(fdir,,sep=),
'1' = c(-1,-1),
'4' = c(-1,0),
'7' = c(-1,1),
'8' = c(0,1),
'9' = c(1,1),
'6' = c(1,0),
'3' = c(1,-1),
'2' = c(0,-1),
'5' = c(0,0),
c(-2,-2)
)
} # 'lookupfdir' end
fdirPlot - function(imData,fdirs) {
dimDat- dim(imData)
Ly- dimDat[1] - 1
Lx- dimDat[2] - 1
centres.x - seq(0,1, length.out=dimDat[2])
centres.y - seq(0,1, length.out=dimDat[1])
image(t(imData),axes=FALSE)
for (i in centres.x) { # x direction
for (j in centres.y) { # y direction
dxdy - lookupfdir(fdirs[i+1,j+1])
dx - dxdy[1]/Lx
dy - dxdy[2]/Ly
if (dxdy[1] != -2)
lines( rbind( c(i, j), c(i + dx, j + dy) ) )
} # IF end
} # FOR j end
} # 'fdirPlot' end
# Testing
imData - matrix(rnorm(30), ncol=3, nrow=2)
fdirs - matrix(9, ncol=3, nrow=2)
fdirPlot(imData, fdirs)
- END reproducible example
However, I would like to do the same as before but over a raster file:
# creating the two rasters from the scratch
library(raster)
imData - raster(nrows=2, ncols=3, xmn=0, xmx=10, ymn=0, ymx=5)
values(imData) - rnorm(6)
fdirs - raster(nrows=2, ncols=3, xmn=0, xmx=10, ymn=0, ymx=5)
values(r) - 9
However, I do not know how to plot the lines over the raster and then
save the resulting lines as a shapefile.
Do you -or somebody else- have any advice about how to do it ?
Thank in advance for your help,
Mauricio
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission
=
Linux user #454569 -- Ubuntu user #17469
=
It is not enough to have knowledge;
one must also apply it (Goethe)
lookupfdir -function(fdir) {
switch(paste(fdir,,sep=),
'1' = c(-1,-1),
'2' = c(-1,0),
'3' = c(-1,1),
'4' = c(0,1),
'5' = c(1,1),
'6' = c(1,0),
'7' = c(1,-1),
'8' = c(0,-1),
'9' = c(0,0),
c(-2,-2)
)
}
fdirPlot - function(imData,fdirs) {
dimDat- dim(imData)
image(imData,axes=FALSE)
for (i in 1:dimDat[1]) {
for (j in 1:dimDat[2]) {
dxdy - lookupfdir(fdirs[i,j])
if (dxdy[1] != -2) {
lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-dxdy[1]-0.5,j-dxdy[2]-0.5)/dimDat[1]))
} else {
if (i==1) {
lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-1.5,j-0.5)/dimDat[1]))
} else {
if (i==dimDat[1]) {
lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i+1-0.5,j-0.5)/dimDat[1]))
} else {
if (j==1) {
lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-0.5,j-1-0.5)/dimDat[1]))
} else {
if (j==dimDat[2]) {
lines(rbind(c(i-0.5,j-0.5)/dimDat[1],c(i-0.5,j+1-0.5)/dimDat[1]))
}
}
}
}
}
}
}
}
Gavan McGrath
___
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] raster directions to vectorial lines (Stream to Feature)
Dear List,
I have a raster file representing flow directions, and I would like to
know if is there any function for transforming this raster file into a
vectorial representation (with lines)
The original raster only has 9 possible values:
1: 225 deg (South-West)
2: 270 deg (South)
3: 315 deg (South-East)
4: 180 deg (West)
5: NA (no flow)
6: 0 deg (East)
7: 135 deg (North-West)
8: 90 deg (North)
9: 45 deg (North-East)
I can use the function 'reclassify' provided by the raster package to
re-class the 9 values into the corresponding decimal degrees.
However, I don't know how to create a line from each cell centre to the
next cell, following the direction stored in each cell.
My colleagues do this with the 'Stream to Feature' tool of Spatial
Analyst (ArcGis 10.0,
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z005800.htm),
but I would like to do it in R.
I would highly appreciate any ideas about how to carry out this task in R.
Thanks in advance,
Mauricio Zambrano-Bigiarini
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:10}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] make_EPSG chile
On 27/01/13 22:09, Orietta Nicolis wrote:
Dear R users,
I'm trying to run the following example with data of Chile, but I didn't
find a EPSG code for this country which could be used in R.
Hi Orietta,
AFAIK, in Chile there are two main local projections:
PSAD56 / UTM zone 18S, and
PSAD56 / UTM zone 19S
which are used depending on the latitude you are working on.
For PSAD56_Z19S, you may try:
library(sp)
psad56.p4s - CRS(+proj=utm +zone=19 +south +ellps=intl +units=m +no_defs)
which is equivalent to use EPSG:24879
if you need higher precision, you may need to add the '+towgs84'
argument to 'psad56.p4s'
Saludos,
Mauricio
PS,
I'm happy to know there are people in Chile using R !
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as sta-
ting an official position of the European Commission
=
Linux user #454569 -- Ubuntu user #17469
=
True courage is a result of reasoning.
A brave mind is always impregnable
(Jeremy Collier)
datahttp://inside-r.org/r-doc/utils/data(meuse)
coordinates(meuse)-~x+y proj4string(meuse)- CRS('+init=epsg:28992')
# Line data
datahttp://inside-r.org/r-doc/utils/data(meuse.grid)
coordinates(meuse.grid)-chttp://inside-r.org/r-doc/base/c('x','y')
meuse.grid-ashttp://inside-r.org/r-doc/methods/as(meuse.grid,
'SpatialPixelsDataFrame')
im-as.image.SpatialGridDataFrame(meuse.grid['dist'])
cl-ContourLines2SLDF(contourLineshttp://inside-r.org/r-doc/grDevices/contourLines
(im))
proj4string(cl)- CRS('+init=epsg:28992')
I run the script:
make_EPSG(Chile)
but it gives the following error:
In file(file, open = r) :
cannot open file 'Chile': No such file or directory
Please, can anyone suggest me how to solve this problem?
Thank you very much,
Orietta
[[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] make_EPSG chile
On 28/01/13 11:48, Roger Bivand wrote: On Mon, 28 Jan 2013, Mauricio Zambrano-Bigiarini wrote: On 27/01/13 22:09, Orietta Nicolis wrote: Dear R users, I'm trying to run the following example with data of Chile, but I didn't find a EPSG code for this country which could be used in R. Hi Orietta, AFAIK, in Chile there are two main local projections: PSAD56 / UTM zone 18S, and PSAD56 / UTM zone 19S But: EPSG - make_EPSG() EPSG[grep(Chile, EPSG$note), 1:2] code note 450 5360 # SIRGAS-Chile 2697 5361 # SIRGAS-Chile / UTM zone 19S 2698 5362 # SIRGAS-Chile / UTM zone 18S I do not get those values: library(rgdal) Loading required package: sp rgdal: version: 0.8-4, (SVN revision 431) Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 1.8.1, released 2011/07/09 Path to GDAL shared files: /usr/share/gdal GDAL does not use iconv for recoding strings. Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470] Path to PROJ.4 shared files: (autodetected) EPSG - make_EPSG() EPSG[grep(Chile, EPSG$note), 1:2] [1] code note 0 rows (or 0-length row.names) which is very likely due to the fact I'm using PROJ4 4.7.1 :( all of which are WGS84 - do we know whether the +ellps in Orietta's data is GRS80 or your earlier intl? She didn't mention it. I was working in Chile up to 2006, and at that time most of the cartography I saw from public institutions (related to water) used intl as ellipsoid. However, last December I realized the cartography (related to water) used WGS84. With PROJ.4 4.8.0, we get a later EPSG too, with: CRS(+init=epsg:24879) CRS arguments: +init=epsg:24879 +proj=utm +zone=19 +south +ellps=intl +towgs84=-288,175,-376,0,0,0,0 +units=m +no_defs Thank you very much for this info. I'll try to update my PROJ4. which can be transformed to WGS84. I recently used that projection. Then, when I transformed some features to WGS84 for visualizing in Google Earth, all of them were in close agreement with GE. Mauricio Roger which are used depending on the latitude you are working on. For PSAD56_Z19S, you may try: library(sp) psad56.p4s - CRS(+proj=utm +zone=19 +south +ellps=intl +units=m +no_defs) which is equivalent to use EPSG:24879 if you need higher precision, you may need to add the '+towgs84' argument to 'psad56.p4s' Saludos, Mauricio PS, I'm happy to know there are people in Chile using R ! ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] [plain text] Error: cannot allocate vector of size 274.7 Mb
On 13/11/12 11:20, Raffaele Morelli wrote:
Hi,
sorry for double posting but I saw my messages were scrubbed from html
garbage and don't know if nobody replies for that, but just in case...
apologize if not.
I am having trouble with this error when calculating a density on a
point pattern dataset with 34480 points.
My script exit when calculating
z-density(mypattern, 3000, dimyx=3000)
The machine the script it's running on it's a debian/amd64 OS with
16Gb RAM, I have set to unlimited the memory limits (the output of
ulimit -a follows)
Any ideas?
Try:
?memory.limit
?gc
and read 'Circle 2' of the R_inferno
(www.burns-stat.com/pages/Tutor/R_inferno.pdf)
IHTH,
Mauricio
--
=
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:11}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Saving a variogram in gstat
The following reproducible example works for me (please run it in a new session): START -- library(gstat) data(meuse) # no trend: coordinates(meuse) = ~x+y myvgm - variogram(log(zinc)~1, meuse) # saving the variogram to your home directory setwd(~) save(myvgm, file=Myvgm.RData) ## remove (almost) everything in the working environment. ## IMPORTANT:: You will get no warning, so ## don't do this unless you are really sure. rm(list = ls()) ## Loading the variogram load(Myvgm.RData) summary(myvgm) END -- IHTH, Mauricio -- = Water Resources Unit Institute for Environment and Sustainability (IES) Joint Research Centre (JRC), European Commission TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), IT webinfo: http://floods.jrc.ec.europa.eu/ = DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as sta- ting an official position of the European Commission = Linux user #454569 -- Ubuntu user #17469 = It is a mistake to look too far ahead. Only one link in the chain of destiny can be handled at a time (Sir Winston Churchill) On 05/11/12 15:04, carolina monmany wrote: Thanks, Sarah. I've tried that and it looks like it is loaded but when I try to do anything such as plot(variog1) or fit it it does not work. When I ask class(variog1) it says character 2012/11/5 Sarah Gosleesarah.gos...@gmail.com If you use save() to save your variograms, you need to use load() to load them. load(variog1.Rdata) Sarah On Monday, November 5, 2012, carolina monmany wrote: Dear all, I am using R 2.15.1 and the package gstat in Linux to calculate variograms from large datasets (200x200m high resolution satellite images). I need to save all the variograms for later modelling and others and I realized late that I was using the wrong command: save(variog1, file=variog1.Rdata) the variograms were being indeed saved but I cannot open them again because they not being saved as variograms but as characters. I tried to save them using the suggested command: data(variog1): 'variog1'; with no success. My main question is: is there a way I can recover the already saved variograms? And the other question : what is the correct syntaxis for saving variograms in gstat? Thanks a lot. -- --- A. CAROLINA MONMANY Universidad de Puerto Rico Departamento de Biologia - CN 235 POBOX 23360 San Juan, Puerto Rico 00931-3360 Tel: +1 787 764 x2847 Fax: +1 787 764 2610 [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Sarah Goslee http://www.stringpage.com http://www.sarahgoslee.com http://www.functionaldiversity.org ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Saving a variogram in gstat
Hi Carolina,
You can save all the variables of your current session by using 'save.image':
save.image(myVGMsession.RData)
and then load it as before (from Windows or GNU/Linux):
load(myVGMsession.RData)
Kinds,
Mauricio
2012/11/5 carolina monmany acmonm...@gmail.com:
I see, Sarah. Thank you!
Now I need to transfer those variograms to my PC running windows and that
was the reason I was doing
save(variog11.2, file=tuc11.2.200.var.Rdata)
and transfering those files from Linux to Windows. Should I transfer the
.Rdata in the directory instead? Where are the objects themselves?
2012/11/5 Sarah Goslee sarah.gos...@gmail.com
On Mon, Nov 5, 2012 at 10:26 AM, carolina monmany acmonm...@gmail.com
wrote:
Thanks Sarah and Mauricio,
here's the example
library(gstat)
setwd(/home/acmonmany/matricesDN)
tuc11.2.frm = read.table(plot11_2.txt, header=T)
coordinates(tuc11.2.frm) - c(X, Y)
proj4string(tuc11.2.frm) -CRS(+proj=utm +zone=20 +south +ellps=intl
+towgs84=-355,21,72,0,0,0,0 +units=m +no_defs)
variog11.2 - variogram(reflectance~1, data=tuc11.2.frm)
save(variog11.2, file=tuc11.2.200.var.Rdata)
NOW when I try to load the saved file (I do not change directory and I
can
see the files there, in the directory stated above):
load(tuc11.2.200.var.Rdata)
class(tuc11.2.200.var.Rdata)
[1] character
Of course it is:
tuc11.2.200.var.Rdata is a character string, as is anything in quotes.
summary(tuc11.2.200.var)
Error in summary(tuc11.2.200.var) :
error in evaluating the argument 'object' in selecting a method for
function 'summary': Error: object 'tuc11.2.200.var' not found
That's not anything, just the name of your Rdata file. The R object
itself keeps the same name as it had when you saved it: variog11.2
summary(variog11.2)
will show you its properties. If you ls() you'll see there is no R object
called
tuc11.2.200.var
If you want your saved object to have that name, you need to rename it
before saving.
Sarah
--
Sarah Goslee
http://www.functionaldiversity.org
--
---
A. CAROLINA MONMANY
Universidad de Puerto Rico
Departamento de Biologia - CN 235
POBOX 23360
San Juan, Puerto Rico 00931-3360
Tel: +1 787 764 x2847
Fax: +1 787 764 2610
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
=
Water Resources Unit
Institute for Environment and Sustainability
Joint Research Centre, European Commission
webinfo: http://floods.jrc.ec.europa.eu/
=
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:13}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] testing for a valid raster format before reading
Dear list,
Do you know if there is any way for testing if a file can be read by the
'raster' command of the 'raster' package ?
something like:
is.raster(myfile.asc)
I'm asking because I have to read many files within several directories.
The filenames of the maps go from .001 up to XX99.999, with
some additional files that should be skipped of the reading process
(which extension varies from directory to directory, e.g., .csv, .txt,
.zip, etc).
The error I'm getting now, which appears when trying to get information
about the file with 'GDALinfo', is the following:
Error in .local(.Object, ...) :
`lai.zip' not recognised as a supported file format.
1 traceback()
15: .Call(RGDAL_OpenDataset, as.character(filename), TRUE, silent,
PACKAGE = rgdal)
14: .local(.Object, ...)
13: initialize(value, ...)
12: initialize(value, ...)
11: new(GDALReadOnlyDataset, filename, silent = silent)
10: GDAL.open(fname, silent = silent)
9: GDALinfo(maps.list[i])
Thanks in advance for any help,
Kinds,
Mauricio Zambrano-Bigiarini
--
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo: http://floods.jrc.ec.europa.eu/
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:10}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] testing for a valid raster format before reading
On 27/07/12 11:35, Barry Rowlingson wrote:
On Fri, Jul 27, 2012 at 10:08 AM, Mauricio Zambrano-Bigiarini
mauricio.zambr...@jrc.ec.europa.eu wrote:
Dear list,
Do you know if there is any way for testing if a file can be read by the
'raster' command of the 'raster' package ?
something like:
is.raster(myfile.asc)
I'm asking because I have to read many files within several directories.
The filenames of the maps go from .001 up to XX99.999, with
some additional files that should be skipped of the reading process
(which extension varies from directory to directory, e.g., .csv, .txt,
.zip, etc).
The error I'm getting now, which appears when trying to get information
about the file with 'GDALinfo', is the following:
Error in .local(.Object, ...) :
`lai.zip' not recognised as a supported file format.
1 traceback()
15: .Call(RGDAL_OpenDataset, as.character(filename), TRUE, silent,
PACKAGE = rgdal)
14: .local(.Object, ...)
13: initialize(value, ...)
12: initialize(value, ...)
11: new(GDALReadOnlyDataset, filename, silent = silent)
10: GDAL.open(fname, silent = silent)
9: GDALinfo(maps.list[i])
Thanks in advance for any help,
You can use 'try' to run code and catch errors. See help(try) for more:
for(f in files){
r = try(raster(f))
if(inherits(r, try-error)){
warning(couldnt read ,f)
}else{
print(summary(r))
}
}
Thanks Barry for your feedback.
Adding the 'silent' argument to 'try' produced the behaviour I wanted:
for(f in files){
r = try(raster(f), silent=TRUE)
if(inherits(r, try-error)){
warning(couldnt read ,f)
}else{
print(summary(r))
}
}
Length ClassMode
66548 RasterLayer S4
Warning messages:
1: In eval(expr, envir, enclos) : couldnt read LHOAT-AGNPS2012.txt
2: In eval(expr, envir, enclos) : couldnt read SPSO.zip
All the best,
Mauricio
--
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo: http://floods.jrc.ec.europa.eu/
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:10}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Plot spatial time series
On 12/06/12 08:40, Tomislav Hengl wrote:
Connected with this topic, there is now an (experimental) plotKML method
for spatial series of raster bricks that allows visualization of
time-series rasters in Google Earth:
http://plotkml.r-forge.r-project.org/RasterBrickTimeSeries-class.html
Here is an example:
http://plotkml.r-forge.r-project.org/Fig_RasterBrickTimeSeries.jpg
http://plotkml.r-forge.r-project.org/LST.ts.kml
plotKML should be soon available via CRAN.
Thank you very much Tom for this info.
This is just what I was looking for.
Cheers,
Mauricio Zambrano-Bigiarini
--
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo: http://floods.jrc.ec.europa.eu/
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:12}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Crop a raster using a shapefile
On 01/06/12 20:49, Thiago Veloso wrote: Dear all, When cropping a raster using a shapefile, through 'crop' function from raster package, only the extent is taken into account: r-raster(lai2011361.Lai_1km.tif) extent(r) class : Extent xmin: -104.4326 xmax: -29.99736 ymin: -40.00064 ymax: 10 br-readShapePoly(/mnt/disco3/MODIS/shapes/brazil.shp) extent(br) class : Extent xmin: -73.83943 xmax: -34.8581 ymin: -33.77086 ymax: 5.38289 c-crop(r,br) extent(c) class : Extent xmin: -73.83772 xmax: -34.85554 ymin: -33.76852 ymax: 5.38428 However, the cropped raster is rectangular. Is there any way to keep also the shape of a shapefile (and not only its extent) when performing this kind of crop operation? You may use a combination of crop and mask: START --- library(raster) # Reading the shapefile myshp - readOGR(myshapefile.shp, layer=basename(myshapefile) ) # Getting the spatial extent of the shapefile e - extent(myshp) # Reading the raster you want to crop myraster - raster(myraster.filename) # Cropping the raster to the shapefile spatial extent myraster.crop - crop(raster, e, snap=out) # Dummy raster with a spatial extension equal to the cropped raster, # but full of NA values crop - setValues(myraster.crop, NA) # Rasterize the catchment boundaries, with NA outside the catchment boundaries myshp.r - rasterize(myshp, crop) # Putting NA values in all the raster cells outside the shapefile boundaries myraster.masked - mask(x=myraster.crop, mask=myshp.r) END --- IHTH, Mauricio -- Water Resources Unit Institute for Environment and Sustainability (IES) Joint Research Centre (JRC), European Commission webinfo: http://floods.jrc.ec.europa.eu/ DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as sta- ting an official position of the European Commission Linux user #454569 -- Ubuntu user #17469 The only things in life you regret, Are the risks that you didn't take (Anonymous) Thanks in advance, Thiago. ___ 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] writeOGR: ID field for exporting a KML file
Dear list,
I'm trying to export a point shapefile into a KML file by using the
'writeOGR' function of the 'rgdal' package.
Everything works fine, except that I'm not able to select the field to
be used for labelling the KML points from the table of the original
shapefile.
I'm able to select the ID field to be used in the KML file by using the
'ogr2ogr' program of GDAL, but I would like to be able to do the same
directly within R with the 'writeOGR' function.
Below you can find a reproducible example:
START -
require(rgdal)
# creating the initial shapefile for the example
cities - readOGR(system.file(vectors, package = rgdal)[1], cities)
writeOGR(cities, cities.shp, cities, driver=ESRI Shapefile)
OGRstring - paste(ogr2ogr -f KML , cities.kml, , cities.shp,
-dsco NameField=COUNTRY, sep = )
system(OGRstring)
END -
I would like to achieve the same final KML file, but using:
writeOGR(cities, dsn=cities2.kml, layer=cities, driver=KML )
but I don't know how to pass the -dsco NameField=COUNTRY argument to
the 'writeOGR' function.
Any help will be very much appreciated.
sessionInfo()
R version 2.15.0 (2012-03-30)
Platform: x86_64-redhat-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_GB.utf8 LC_NUMERIC=C
[3] LC_TIME=en_GB.utf8LC_COLLATE=en_GB.utf8
[5] LC_MONETARY=en_GB.utf8LC_MESSAGES=en_GB.utf8
[7] LC_PAPER=CLC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.utf8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] rgdal_0.7-11 sp_0.9-99
loaded via a namespace (and not attached):
[1] grid_2.15.0lattice_0.20-6 tools_2.15.0
Thanks in advance,
Mauricio Zambrano-Bigiarini
--
Water Resources Unit
Institute for Environment and Sustainability (IES)
Joint Research Centre (JRC), European Commission
webinfo: http://floods.jrc.ec.europa.eu/
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:11}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot question in sp_0.9-98
2012/4/3 Edzer Pebesma edzer.pebe...@uni-muenster.de: Thanks for notifying this. Removing the cex=.7 solves this, as does replicating cex, as in spplot(meuse, c(cadmium, copper, lead, zinc), do.log = TRUE, key.space = right, as.table = TRUE, sp.layout=list(rv, scale, text1, text2, arrow), main = Heavy metals (top soil), ppm, cex=rep(.7,4*length(meuse)), cuts = cuts) This, and returning the color vector, are obviously bugs, and will be addressed in the next version of sp. Thanks for your answer Edzer. Yesterday I could confirm that this issue is not present in sp_0.9-97 All the best, Mauricio Zambrano-Bigiarini -- === Water Resources Unit Institute for Environment and Sustainability Joint Research Centre, European Commission webinfo: http://floods.jrc.ec.europa.eu/ === DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as stating an official position of the European Commission === Linux user #454569 -- Ubuntu user #17469 There is only one pretty child in the world, and every mother has it. (Chinese Proverb) http://c2.com/cgi/wiki?HowToAskQuestionsTheSmartWay On 04/03/2012 08:37 PM, Mauricio Zambrano-Bigiarini wrote: Dear List, I just notice an strange behaviour of the 'spplot' function, and I would like your help to figure out if it is something related to my system or is something related to the 'spplot' function. Below goes a reproducible example taken from http://r-spatial.sourceforge.net/gallery/ : # fig06.R multi-panel plot, scales + north arrow only in last plot. --START--- library(sp) library(lattice) # required for trellis.par.set(): trellis.par.set(sp.theme()) # sets color ramp to bpy.colors() data(meuse) coordinates(meuse)=~x+y data(meuse.riv) meuse.sr = SpatialPolygons(list(Polygons(list(Polygon(meuse.riv)),meuse.riv))) rv = list(sp.polygons, meuse.sr, fill = lightblue) ## multi-panel plot, scales + north arrow only in last plot: ## using the which argument in a layout component ## (if which=4 was set as list component of sp.layout, the river ## would as well be drawn only in that (last) panel) scale = list(SpatialPolygonsRescale, layout.scale.bar(), offset = c(180500,329800), scale = 500, fill=c(transparent,black), which = 4) text1 = list(sp.text, c(180500,329900), 0, cex = .5, which = 4) text2 = list(sp.text, c(181000,329900), 500 m, cex = .5, which = 4) arrow = list(SpatialPolygonsRescale, layout.north.arrow(), offset = c(181300,329800), scale = 400, which = 4) cuts = c(.2,.5,1,2,5,10,20,50,100,200,500,1000,2000) spplot(meuse, c(cadmium, copper, lead, zinc), do.log = TRUE, key.space = right, as.table = TRUE, sp.layout=list(rv, scale, text1, text2, arrow), main = Heavy metals (top soil), ppm, cex = .7, cuts = cuts) --END--- After the last spplot command, I only get one point in the upper left figure, which is very different from the Figure 6 shown in http://r-spatial.sourceforge.net/gallery/ In addition, I get the following output after the spplot command: [1] #A714EBFF #6500 #6500 #2400 #2400 #2400 [7] #2400 #2400 #2400 #DAFF #DAFF #DAFF [13] #A714EBFF #2400 #DAFF #6500 #6500 #6500 [19] #6500 #A714EBFF #6500 #2400 #2400 #DAFF [25] #DAFF #DAFF #DAFF #DAFF #DAFF #DAFF [31] #DAFF #DAFF #2400 #DAFF #DAFF #2400 [37] #6500 #6500 #6500 #A714EBFF #2400 #DAFF [43] #DAFF #DAFF #2400 #2400 #2400 #DAFF [49] #DAFF #2400 #DAFF #6500 #A714EBFF #A714EBFF [55] #6500 #6500 #2400 #2400 #A714EBFF #6500 [61] #6500 #6500 #6500 #6500 #6500 #6500 [67] #6500 #86FF #2400 #2400 #2400 #2400 [73] #2400 #2400 #2400 #2400 #2400 #2400 [79] #6500 #6500 #A714EBFF #A714EBFF #6500 #2400 [85] #DAFF #2400 #2400 #2400 #2400 #DAFF [91] #DAFF #2400 #2400 #86FF #86FF #33FF [97] #86FF #33FF #33FF #33FF #86FF #33FF [103] #33FF #33FF #33FF #33FF #33FF #33FF [109] #33FF #33FF #33FF #33FF #33FF #33FF [115] #33FF #33FF #33FF #2400 #33FF #33FF [121] #33FF #33FF #DAFF #2400 #33FF #33FF [127] #33FF #33FF #33FF #DAFF #86FF #DAFF [133] #33FF #33FF #86FF #86FF #86FF #DAFF [139] #DAFF #DAFF #DAFF #86FF #86FF #2400 [145] #2400
Re: [R-sig-Geo] Generating a precise map from coordinates.
On 19/03/12 04:50, Agus Camacho wrote: Hello all, I am still new to plot maps in R and trying to plot the object (wrld_simpl) from maptools. However found it difficult to plot a map with precise coordinates. The following series of plots shows what I am talking about: The margins of the plot are those specified only in some cases but others not. Anybody could give me a hint to generate maps from coordinates without having to be fishing a set of coordinates that looks good? Did you look at : http://r-spatial.sourceforge.net/gallery/ http://cran.r-project.org/web/views/Spatial.html ? IHTH, Mauricio Zambrano-Bigiarini -- === FLOODS Action Water Resources Unit (H01) Institute for Environment and Sustainability (IES) European Commission, Joint Research Centre (JRC) TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), Italy webinfo: http://floods.jrc.ec.europa.eu/ work-phone : (+39)-(0332)-789588 work-fax : (+39)-(0332)-786653 === DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as stating an official position of the European Commission. === Linux user #454569 -- Ubuntu user #17469 === If Columbus had turned back, no one would have blamed him. Of course, no one would have remembered him either. (Source Unknown) plot(wrld_simpl, xlim=c(100,160), ylim=c(-45,10), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(100,160), ylim=c(-40,-10), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(100,160), ylim=c(-40,-15), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(100,160), ylim=c(-40,-20), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(100,160), ylim=c(-40,-25), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(115,160), ylim=c(-40,-25), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(115,155), ylim=c(-40,-25), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(115,155), ylim=c(-40,-20), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(115,155), ylim=c(-40,-15), axes=TRUE,col='light yellow') plot(wrld_simpl, xlim=c(115,155), ylim=c(-40,-14), axes=TRUE,col='light yellow') ___ 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] spplot with legend only, and legend inside the plotting area
On 08/03/12 09:21, Oscar Perpiñan wrote: Sorry: I forgot to say that you need latticeExtrato use c.trellis.Then: library(latticeExtra) c(dummy, dummy, dummy, dummy). Thanks Oscar !. I didn't realise about 'latticeExtra' instead of lattice... One last question. By using c(), is it possible to put different axis labels for each plot ? -- START library(sp) xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100)) coordinates(xyz)=~x+y # Different spatial extent for the 2nd dummy set xyz2 = data.frame(expand.grid(x=100:200,y=100:200),rnorm(10201)) coordinates(xyz2)=~x+y # only to get the default values from 'auto.key' dummy - spplot(xyz, cex=2, alpha=0.7, xlab=X1, ylab=Y1, scales=list(draw=TRUE)) dummy2 - spplot(xyz2, cex=2, alpha=0.7, xlab=X2, ylab=Y2) # checking that 'dummy' is a trellis class(dummy) library(lattice) library(latticeExtra) cObj - c(dummy, dummy2, dummy, dummy2, layout=c(2,4), merge.legends=FALSE) update(cObj) -- END in the previous figure, xlab and ylab are set the same for all the plots, even if for two of them they are different ... Thanks in advance, Mauricio -- === FLOODS Action Water Resources Unit (H01) Institute for Environment and Sustainability (IES) European Commission, Joint Research Centre (JRC) webinfo: http://floods.jrc.ec.europa.eu/ === DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as stating an official position of the European Commission. === Linux user #454569 -- Ubuntu user #17469 === If Columbus had turned back, no one would have blamed him. Of course, no one would have remembered him either. (Source Unknown) ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot with legend only, and legend inside the plotting area
On 07/03/12 17:38, Mauricio Zambrano-Bigiarini wrote:
I'm almost there...
The real problem appears when I try to produce several spplots within
the same figure (each one of them for a different spatial area), with a
unique legend for all the figures (that is the reason why I can not use
'legend' command).
For example, while trying to produce 4 different plots within the same
figure:
library(sp)
xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100))
coordinates(xyz)=~x+y
# only to get the default values from 'auto.key'
dummy - spplot(xyz, cex=10, alpha=0.7)
# This has auto.key = bottom by default.
cols = dummy$legend$bottom$args$key$points$col
txt = dummy$legend$bottom$args$key$text[[1]]
# number of elements in the legend
n - length(cols)
p1 - spplot(xyz, cex=10, alpha=0.7, auto.key=FALSE)
p2 - p1
p3 - p1
# only the legend
p4 - spplot(xyz, cex=0, auto.key=list(x=.5, y=.5, corner = c(0.5, 0.5),
title=my legend) )
print(p1, split = c(1, 1, 2, 2), more = TRUE)
print(p2, split = c(2, 1, 2, 2), more = TRUE)
print(p3, split = c(1, 2, 2, 2), more = TRUE)
print(p4, split = c(2, 2, 2, 2), more = FALSE)
The previous solution has two problems (for me) that I couldn't solve:
1) the colors for the symbols are not correct
2) The symbols are plotted at the right side of the legend, while I
would like to have them on the left side
I just solved the two previous problems, and I share the solution now
just in case it be useful for somebody else in the future:
library(sp)
mydata - rnorm(100)
xyz = data.frame(expand.grid(x=1:10,y=1:10), mydata)
coordinates(xyz)=~x+y
# only to get the default values from 'auto.key'
dummy - spplot(xyz, cex=10, alpha=0.7)
# This has auto.key = bottom by default.
cols = dummy$legend$bottom$args$key$points$col
txt = dummy$legend$bottom$args$key$text[[1]]
# number of elements in the legend
n - length(cols)
p1 - spplot(xyz, cex=10, alpha=0.7, auto.key=FALSE)
p2 - p1
p3 - p1
# legend levels
cuts - unique( quantile( as.numeric(mydata), probs=c(0.25, 0.5, 0.75,
1), na.rm=TRUE) )
gof.levels - cut(mydata, cuts)
# only the legend, in an empty plot
library(lattice)
p4 - xyplot(1~1, type=n, xlab=, ylab=,
groups=gof.levels,
scales=list(draw=FALSE),
# automatic legend
key = list(x = .5, y = .5, corner = c(0.5, 0.5),
title=legend,
points = list(pch=16, col=cols, cex=1.5),
text = list(as.character(txt))
),
# removing outter box.
#From:
https://stat.ethz.ch/pipermail/r-help/2007-September/140098.html
par.settings = list(axis.line = list(col = transparent)),
axis = function(side, ...) {
axis.default(side = side, ...)
},
)
# Creating the final figure
print(p1, split = c(1, 1, 2, 2), more = TRUE)
print(p2, split = c(2, 1, 2, 2), more = TRUE)
print(p3, split = c(1, 2, 2, 2), more = TRUE)
print(p4, split = c(2, 2, 2, 2), more = FALSE)
I tried to solve the 2nd problem by using 'key' instead of 'auto.key':
p4 - spplot(xyz, cex=0,
key = list(x = 0.5, y = 0.5, corner = c(0.5, 0.5),
title=my legend, text = list(txt),
points = list(pch = rep(16, length = n),
col = cols,
cex = rep(1.5, length = n)
)
)
)
print(p1, split = c(1, 1, 2, 2), more = TRUE)
print(p2, split = c(2, 1, 2, 2), more = TRUE)
print(p3, split = c(1, 2, 2, 2), more = TRUE)
print(p4, split = c(2, 2, 2, 2), more = FALSE)
and now the colors are correct, but the legend appears at the bottom of
the plotting area, not in the centre as before.
Does it mean that spplot ignores the x, y, and corner arguments in 'key' ?
sessionInfo() # part of it
R version 2.14.1 (2011-12-22)
Platform: x86_64-redhat-linux-gnu (64-bit)
.
.
.
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sp_0.9-95 lattice_0.20-0
loaded via a namespace (and not attached):
[1] cluster_1.14.2 grid_2.14.1 Hmisc_3.9-2
Thanks in advance,
Mauricio Zambrano-Bigiarini
--
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
TP 261, Via Enrico Fermi 2749, 21027 Ispra (VA), Italy
webinfo: http://floods.jrc.ec.europa.eu/
work-phone : (+39)-(0332)-789588
work-fax : (+39)-(0332)-786653
===
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:11}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot with legend only, and legend inside the plotting area
On 08/03/12 18:22, Oscar Perpiñán Lamigueiro wrote:
Mauricio Zambrano-Bigiarinimauricio.zambr...@jrc.ec.europa.eu writes:
Thanks Oscar !. I didn't realise about 'latticeExtra' instead of lattice...
One last question. By using c(), is it possible to put different axis
labels for each plot ?
-- START
library(sp)
xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100))
coordinates(xyz)=~x+y
# Different spatial extent for the 2nd dummy set
xyz2 = data.frame(expand.grid(x=100:200,y=100:200),rnorm(10201))
coordinates(xyz2)=~x+y
# only to get the default values from 'auto.key'
dummy- spplot(xyz, cex=2, alpha=0.7, xlab=X1, ylab=Y1,
scales=list(draw=TRUE))
dummy2- spplot(xyz2, cex=2, alpha=0.7, xlab=X2, ylab=Y2)
# checking that 'dummy' is a trellis
class(dummy)
library(lattice)
library(latticeExtra)
cObj- c(dummy, dummy2, dummy, dummy2, layout=c(2,4), merge.legends=FALSE)
update(cObj)
-- END
Yes, it is possible. You have to modify the parameters with update:
cObj- c(dummy, dummy2, dummy, dummy2)
update(cObj, xlab=c('X1', 'X2'), ylab.right='Y2')
Best,
Thanks Oscar !.
Now I know two solutions for the initial problem :)
All the best,
Mauricio
--
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as stating
an official position of the European Commission.
===
Linux user #454569 -- Ubuntu user #17469
===
If Columbus had turned back, no one would have blamed him.
Of course, no one would have remembered him either.
(Source Unknown)
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot with legend only, and legend inside the plotting area
On 07/03/12 13:11, Oscar Perpiñán Lamigueiro wrote: Mauricio Zambrano-Bigiarinimauricio.zambr...@jrc.ec.europa.eu writes: I need to produce two different spplot figures of spatial points, and I would like to ask your help for the creation of the second one: 1) the first figure has to show the values but without any legend: library(sp) xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100)) coordinates(xyz)=~x+y spplot(xyz, cex=10, alpha=0.7, auto.key=FALSE) (I don't have any problem with this) 2) The second figure has to show ONLY the legend for the values of the previous figure (without plotting the points), inside the plotting area. I tried with: spplot(xyz, cex=0, alpha=0.7, auto.key=TRUE, key.space=center) Hello, Try this: spplot(xyz, auto.key=list(x=.5, y=.5), cex=0) Change the values of x,y for your convenience. You can also find useful corner instead of x,y. For more information read the auto.key, key and legend parts of the help page of xyplot. Cheers, Oscar. Thank you very much you both Oscar and Jon for your help. I'm almost there... The real problem appears when I try to produce several spplots within the same figure (each one of them for a different spatial area), with a unique legend for all the figures (that is the reason why I can not use 'legend' command). For example, while trying to produce 4 different plots within the same figure: library(sp) xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100)) coordinates(xyz)=~x+y # only to get the default values from 'auto.key' dummy - spplot(xyz, cex=10, alpha=0.7) # This has auto.key = bottom by default. cols = dummy$legend$bottom$args$key$points$col txt = dummy$legend$bottom$args$key$text[[1]] # number of elements in the legend n - length(cols) p1 - spplot(xyz, cex=10, alpha=0.7, auto.key=FALSE) p2 - p1 p3 - p1 # only the legend p4 - spplot(xyz, cex=0, auto.key=list(x=.5, y=.5, corner = c(0.5, 0.5), title=my legend) ) print(p1, split = c(1, 1, 2, 2), more = TRUE) print(p2, split = c(2, 1, 2, 2), more = TRUE) print(p3, split = c(1, 2, 2, 2), more = TRUE) print(p4, split = c(2, 2, 2, 2), more = FALSE) The previous solution has two problems (for me) that I couldn't solve: 1) the colors for the symbols are not correct 2) The symbols are plotted at the right side of the legend, while I would like to have them on the left side I tried to solve the 2nd problem by using 'key' instead of 'auto.key': p4 - spplot(xyz, cex=0, key = list(x = 0.5, y = 0.5, corner = c(0.5, 0.5), title=my legend, text = list(txt), points = list(pch = rep(16, length = n), col = cols, cex = rep(1.5, length = n) ) ) ) print(p1, split = c(1, 1, 2, 2), more = TRUE) print(p2, split = c(2, 1, 2, 2), more = TRUE) print(p3, split = c(1, 2, 2, 2), more = TRUE) print(p4, split = c(2, 2, 2, 2), more = FALSE) and now the colors are correct, but the legend appears at the bottom of the plotting area, not in the centre as before. Does it mean that spplot ignores the x, y, and corner arguments in 'key' ? sessionInfo() # part of it R version 2.14.1 (2011-12-22) Platform: x86_64-redhat-linux-gnu (64-bit) . . . attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] sp_0.9-95 lattice_0.20-0 loaded via a namespace (and not attached): [1] cluster_1.14.2 grid_2.14.1Hmisc_3.9-2 Thanks in advance, Mauricio Zambrano-Bigiarini -- === FLOODS Action Water Resources Unit (H01) Institute for Environment and Sustainability (IES) European Commission, Joint Research Centre (JRC) webinfo: http://floods.jrc.ec.europa.eu/ === DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as stating an official position of the European Commission. === Linux user #454569 -- Ubuntu user #17469 === If Columbus had turned back, no one would have blamed him. Of course, no one would have remembered him either. (Source Unknown) ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot with legend only, and legend inside the plotting area
2012/3/7 Oscar Perpiñán Lamigueiro oscar.perpi...@upm.es: Mauricio Zambrano-Bigiarini mauricio.zambr...@jrc.ec.europa.eu writes: On 07/03/12 13:11, Oscar Perpiñán Lamigueiro wrote: Mauricio Zambrano-Bigiarinimauricio.zambr...@jrc.ec.europa.eu writes: The real problem appears when I try to produce several spplots within the same figure (each one of them for a different spatial area), with a unique legend for all the figures (that is the reason why I can not use 'legend' command). Then you could use c.trellis: With your example: c(dummy, dummy, dummy, dummy). Inside c() you can define the 'layout' and decide if you want to merge legends with 'merge.legends'. Oscar. Thanks Oscar. i followed the examples in c.trellis, but it doesn't work: library(sp) xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100)) coordinates(xyz)=~x+y # only to get the default values from 'auto.key' dummy - spplot(xyz, cex=10, alpha=0.7) # checking that 'dummy' is a trellis class(dummy) library(lattice) cObj - c(dummy, dummy) update(cObj) Error en eval(expr, envir, enclos) : '...' usado en un contexto incorrecto # EN: '...' used in an invalid context Also, in the documentation of c.trellis, it is mentioned that Many properties of the display, such as titles, axis settings and aspect ratio will be taken from the first object only., so if the spatial extent of the figures are different, 'c.trellis' shouldn't do the work, right ? Thanks again, Mauricio -- === FLOODS Action Water Resources Unit (H01) Institute for Environment and Sustainability (IES) European Commission, Joint Research Centre (JRC) webinfo : http://floods.jrc.ec.europa.eu/ === DISCLAIMER: The views expressed are purely those of the writer and may not in any circumstances be regarded as stating an official position of the European Commission. === Linux user #454569 -- Ubuntu user #17469 === There is only one pretty child in the world, and every mother has it. (Chinese Proverb) === http://c2.com/cgi/wiki?HowToAskQuestionsTheSmartWay ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot handling of overlapping points?
On 05/03/12 14:40, Edzer Pebesma wrote:
This has now been changed/repared in sp on r-forge (svn), and will
appear from sp 0.9-97 on. A test run would be:
library(sp)
xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100))
coordinates(xyz)=~x+y
spplot(xyz, cex=10)
which used to plot points in value order, and does now in data record order.
I put up resulting graphs, for comparison, at
http://ifgi.uni-muenster.de/~epebe_01/xyz.html
As this change may make some of your spplot's for points different, if
there are objections against this change (which I consider an
improvement, if not a bug fix), it is time now to let me / us know.
Would it be possible to set a degree of transparency for the
overlapping points in spplot ?
Thanks in advance,
Mauricio Zambrano-Bigiarini
-
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as stating
an official position of the European Commission.
===
Linux user #454569 -- Ubuntu user #17469
===
If Columbus had turned back, no one would have blamed him.
Of course, no one would have remembered him either.
(Source Unknown)
On 03/02/2012 09:03 PM, MacQueen, Don wrote:
I have a SpatialPointsDataFrame object in which many points are very close
together, such that the markers (plotting characters) tend to overlap. (Of
course, this depends on marker size and the scale at which I plot; if I
zoom in there is less overlap.)
It appears that when markers overlap, spplot() places a marker associated
with the larger value after, and therefore on top of, a marker associated
with a smaller value.
Am I correct? Or more generally, what is the algorithm that determines the
order in which markers are added?
I have searched ?spplot and related help pages and haven't found an
explanation (at least, not yet).
I can add that I don't believe markers are placed in the order in which
they appear in the SpatialPointsDataFrame. I say this because when I
attempt to reproduce an spplot using base graphics plot() I have to sort
from smallest to largest value to succeed.
I can probably provide a small reproducible example if necessary, but I'm
hoping it's not necessary.
Here are my actual commands.
tmps is the SpatialPointsDataFrame.
tmp is coordinates(tmps)
(they have the same number of rows in the same order)
Note that I'm using the cuts argument to break a continuous variable
('cpm2') into bins. I have carefully matched the colors in tmps$col2 with
the cbin.cols object passed to spplot(), and the break points for the bin
boundaries, so I believe that everything else that could affect the final
appearance, other than the order in which the markers are placed, is
controlled.
#1 using spplot
spplot(tmps,c('cpm2'),
key.space='right',
legendEntries=cbin.lbls,
cuts=cbin.brks,
col.regions=cbin.cols,
cex=0.4)
#2 using base graphics
plot(tmp[,1],tmp[,2], asp=1, cex=0.6, pch=16, col=tmps$col2)
Thanks
-Don
-
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] spplot handling of overlapping points?
On 06/03/12 16:55, Edzer Pebesma wrote:
Yes, look for argument alpha e.g. in ?rgb or ?bpy.colors
Thank you very much Edzer. It works perfectly:
library(sp)
xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100))
coordinates(xyz)=~x+y
spplot(xyz, cex=10, alpha=0.7)
Cheers,
Mauricio
--
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as stating
an official position of the European Commission.
===
Linux user #454569 -- Ubuntu user #17469
===
If Columbus had turned back, no one would have blamed him.
Of course, no one would have remembered him either.
(Source Unknown)
On 03/06/2012 04:49 PM, Mauricio Zambrano-Bigiarini wrote:
On 05/03/12 14:40, Edzer Pebesma wrote:
This has now been changed/repared in sp on r-forge (svn), and will
appear from sp 0.9-97 on. A test run would be:
library(sp)
xyz = data.frame(expand.grid(x=1:10,y=1:10),rnorm(100))
coordinates(xyz)=~x+y
spplot(xyz, cex=10)
which used to plot points in value order, and does now in data record
order.
I put up resulting graphs, for comparison, at
http://ifgi.uni-muenster.de/~epebe_01/xyz.html
As this change may make some of your spplot's for points different, if
there are objections against this change (which I consider an
improvement, if not a bug fix), it is time now to let me / us know.
Would it be possible to set a degree of transparency for the
overlapping points in spplot ?
Thanks in advance,
Mauricio Zambrano-Bigiarini
-
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:
The views expressed are purely those of the writer
and may not in any circumstances be regarded as stating
an official position of the European Commission.
===
Linux user #454569 -- Ubuntu user #17469
===
If Columbus had turned back, no one would have blamed him.
Of course, no one would have remembered him either.
(Source Unknown)
On 03/02/2012 09:03 PM, MacQueen, Don wrote:
I have a SpatialPointsDataFrame object in which many points are very
close
together, such that the markers (plotting characters) tend to
overlap. (Of
course, this depends on marker size and the scale at which I plot; if I
zoom in there is less overlap.)
It appears that when markers overlap, spplot() places a marker
associated
with the larger value after, and therefore on top of, a marker
associated
with a smaller value.
Am I correct? Or more generally, what is the algorithm that
determines the
order in which markers are added?
I have searched ?spplot and related help pages and haven't found an
explanation (at least, not yet).
I can add that I don't believe markers are placed in the order in which
they appear in the SpatialPointsDataFrame. I say this because when I
attempt to reproduce an spplot using base graphics plot() I have to sort
from smallest to largest value to succeed.
I can probably provide a small reproducible example if necessary, but
I'm
hoping it's not necessary.
Here are my actual commands.
tmps is the SpatialPointsDataFrame.
tmp is coordinates(tmps)
(they have the same number of rows in the same order)
Note that I'm using the cuts argument to break a continuous variable
('cpm2') into bins. I have carefully matched the colors in tmps$col2
with
the cbin.cols object passed to spplot(), and the break points for the
bin
boundaries, so I believe that everything else that could affect the
final
appearance, other than the order in which the markers are placed, is
controlled.
#1 using spplot
spplot(tmps,c('cpm2'),
key.space='right',
legendEntries=cbin.lbls,
cuts=cbin.brks,
col.regions=cbin.cols,
cex=0.4)
#2 using base graphics
plot(tmp[,1],tmp[,2], asp=1, cex=0.6, pch=16, col=tmps$col2)
Thanks
-Don
-
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] R / ArcGIS
On 17/02/12 08:15, Edzer Pebesma wrote:
Hi,
I might succeed meeting with some ESRI / ArcGIS folks at AAG next week;
I'm in touch with them in particular about linking R to ArcGIS.
I would like to find out if any of you has been successful (or
unsuccessful) in reading or writing ArcGIS geodatabase files directly
from R, by using the rgdal driver, or by some other means.
Some days ago Barry Rowlingson share a link regarding this topic:
http://www.structuralknowledge.com/2012/02/03/why-esri-as-is-cant-be-part-of-the-open-government-movement/
In my opinion it would be great for the spatial community if ESRI helps
to make it possible the open access to their new geodatabase format.
My 2 cents.
All the best,
Mauricio Zambrano-Bigiarini
-
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:20}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Building Windows binary of rgdal package to include PGeo driver for reading ESRI personal geodatabases.
On 07/02/12 16:57, Barry Rowlingson wrote:
On Tue, Feb 7, 2012 at 1:24 PM, Michael Denslow
michael.dens...@gmail.com wrote:
Sorry to join this discuss late, but I have also been curious about
reading ESRI geodatabases since I teach with ArcGIS and work with
rgdal for my own research purposes. I have two questions that I hope
are not too off topic.
1. I have the most recent rgdal version (rgdal_0.7-8) build from
source and Pgeo driver list says FALSE. I am using MAC OS 10.6. Can
you point me to documentation to use this on Mac OS?
2. Alternatively, I noticed that GDAL 1.9 has a new driver for .mdb
(http://trac.osgeo.org/gdal/wiki/Release/1.9.0-News). However this
makes reference to Pgeo.
Will this driver change things significantly?
I am still using GDAL 1.8 downloaded from
http://www.kyngchaos.com/software/frameworks.
Thanks for any thoughts and/or advice,
Michael
I just read your email just after reading this blog entry:
http://www.structuralknowledge.com/2012/02/03/why-esri-as-is-cant-be-part-of-the-open-government-movement/
Thank you very much Barry for sharing such interesting post about the
new proprietary format of ESRI !.
All the best,
Mauricio Zambrano-Bigiarini
--
===
FLOODS Action
Water Resources Unit (H01)
Institute for Environment and Sustainability (IES)
European Commission, Joint Research Centre (JRC)
webinfo: http://floods.jrc.ec.europa.eu/
===
DISCLAIMER:\ The views expressed are purely those of th...{{dropped:22}}
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
