Re: [R-sig-Geo] spplot with legend only, and legend inside the plotting area
Sorry: I forgot to say that you need latticeExtrato use c.trellis.Then: library(latticeExtra) c(dummy, dummy, dummy, dummy). Oscar El miércoles 7 de marzo de 2012, Mauricio Zambrano-Bigiarini hzambran.newsgro...@gmail.com escribió: 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 [[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] counting the depressions(group of negative numbers) over time from a raster brick
Dear all I am trying to count the number of group of depressions (negative values) from a climate data set and have least idea on how to go about it. Let me explain the scenario. I have a raster brick with 468 layers and and each layer has 7458 cells. cntnegclass : RasterBrick dimensions : 66, 113, 7458, 468 (nrow, ncol, ncell, nlayers) resolution : 0.108, 0.108 (x, y) extent : 77.946, 90.15, 24.946, 32.074 (xmin, xmax, ymin, ymax) coord. ref. : NA values : in memory min values : -359.51 -341.21 -315.45 -148.10 -187.39 -52.87 -66.72 -52.17 -286.81 -306.74 ... max values : -7.589 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 ... Now for example lets take the 5000th pixel cntneg[5000] Which will give me 468 values of that pixel over time. [1] -90.795107 -89.990016 -94.8407540.00 -15.0855170.00 [7]0.000.000.00 -12.469657 -114.757702 -115.372023 [13] -107.194478 -92.916680 -115.105817 -113.205776 -115.003430 -62.175070 [19]0.000.000.00 -72.358073 -105.006508 -115.372023 [25] -48.836959 -102.314928 -113.271826 -115.372023 -79.5300550.00 [31]0.000.000.00 -15.048987 -115.208204 -115.372023 [37] -115.003430 -108.757617 -113.122594 -115.372023 -111.699048 -17.618498 [43]0.000.00 Now here I need to do two tasks 1) count number of times the rainfall went below the average - those with negative values. And zeros are with positive RF values (which i converted to zero using reclass) for the ease of calculation. In the above example I want to pick the group of negative numbers and count. ie, (-90.795107 -89.990016 -94.840754), (-15.085517), ( -12.469657 -114.757702 -115.372023, -107.194478 -92.916680 -115.105817 -113.205776 -115.003430 -62.175070), (-72.358073 -105.006508 -115.372023, -48.836959 -102.314928 -113.271826 -115.372023 -79.530055) etc. The resulted layer pixel value should be the count of these groups, which in this case is 5. Like wise need to do for all the pixels along time dimension. 2) For each group I want to pick the minimum values and resulted pixel will have the sum of those minimum values. If a group has one value, keep the same. I am stuck to start with this process. I am assuming I need to convert the brick into dataframe and do this. Can any one help me in giving a lead on how to go about it? Really appreciate any help. Sorry if the explanation is confusing. Regards Sajid -- Sajid Pareeth IWMI, Colombo [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] 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
[R-sig-Geo] Dismo Maxent - Error with evaluate(), missing layers (or wrong names)
Hi, I am running multiple Maxent models for a list of 40 species in a loop. Within the loop I am splitting the presence points into a train dataset and a test dataset. All the model fitting and the prediction works fine, there is just a problem with some species and evaluate(). What I am doing and the error: #witholding a 20%(1/5) sample for testing fold - kfold(subset.occ, k=5) occtest - subset.occ[fold == 1, ] occtrain - subset.occ[fold != 1, ] occtest coordinates species 4025 (3525400, 6055500) N.Gesamt.Meerforelle occtrain coordinates species 4021 (3527250, 6060480) N.Gesamt.Meerforelle 4023 (3526100, 6058550) N.Gesamt.Meerforelle 4056 (3522320, 6048620) N.Gesamt.Meerforelle 4066 (3526840, 6060240) N.Gesamt.Meerforelle 4090 (3527920, 6063450) N.Gesamt.Meerforelle evaluate(me, p=occtest, a=bg, x=predictors) Fehler in .local(object, ...) : missing layers (or wrong names) I guess that is related to the number of occtest-samples. As there is only one sample I might get this error. I just wanted to get a proof if I am right with this guess, as the error message is difficult to interpret to me. thank you! cheers, /johannes -- Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] counting the depressions(group of negative numbers) over time from a raster brick
Hi Sajid,
you don't need to transform anything to a data.frame.
Instead you may use the function 'calc' with a function that calculate
the two values you are looking for.
(see the help page on 'calc').
I think the following function will do what you expect to :
theFun - function(x)
{
start.grp - which (x 0 c(TRUE, x[-length(x)]) = 0)
end.grp - which (x 0 c(x[-1], TRUE) = 0)
nb.grps - length(start.grp)
if (nb.grps 0) {
grps - rep (NA, length(x))
grps[x0] - unlist(mapply (rep, 1:nb.grps,
end.grp-start.grp+1, SIMPLIFY=FALSE))
sum.val - sum (tapply(x[x0], grps[x0], min))
return (c(nb.grps, sum.val))
} else {
return(c(0, 0))
}
}
calc (cntneg, theFun)
Best,
Vlad
---
Vladislav Navel
SCAN
www.scan-datamining.com
Le 08/03/2012 11:46, sajid pareeth a écrit :
Dear all
I am trying to count the number of group of depressions (negative
values) from a climate data set and have least idea on how to go about it.
Let me explain the scenario.
I have a raster brick with 468 layers and and each layer has 7458 cells.
cntnegclass : RasterBrick
dimensions : 66, 113, 7458, 468 (nrow, ncol, ncell, nlayers)
resolution : 0.108, 0.108 (x, y)
extent : 77.946, 90.15, 24.946, 32.074 (xmin, xmax, ymin, ymax)
coord. ref. : NA
values : in memory
min values : -359.51 -341.21 -315.45 -148.10 -187.39 -52.87 -66.72
-52.17 -286.81 -306.74 ...
max values : -7.589 0.000 0.000 0.000 0.000 0.000 0.000
0.000 0.000 0.000 ...
Now for example lets take the 5000th pixel
cntneg[5000]
Which will give me 468 values of that pixel over time.
[1] -90.795107 -89.990016 -94.8407540.00 -15.0855170.00
[7]0.000.000.00 -12.469657 -114.757702 -115.372023
[13] -107.194478 -92.916680 -115.105817 -113.205776 -115.003430 -62.175070
[19]0.000.000.00 -72.358073 -105.006508 -115.372023
[25] -48.836959 -102.314928 -113.271826 -115.372023 -79.5300550.00
[31]0.000.000.00 -15.048987 -115.208204 -115.372023
[37] -115.003430 -108.757617 -113.122594 -115.372023 -111.699048 -17.618498
[43]0.000.00
Now here I need to do two tasks
1) count number of times the rainfall went below the average - those with
negative values. And zeros are with positive RF values (which i converted
to zero using reclass) for the ease of calculation.
In the above example I want to pick the group of negative numbers and
count. ie, (-90.795107 -89.990016 -94.840754), (-15.085517), (
-12.469657 -114.757702 -115.372023, -107.194478 -92.916680 -115.105817
-113.205776 -115.003430 -62.175070), (-72.358073 -105.006508 -115.372023,
-48.836959 -102.314928 -113.271826 -115.372023 -79.530055) etc. The
resulted layer pixel value should be the count of these groups, which in
this case is 5. Like wise need to do for all the pixels along time
dimension.
2) For each group I want to pick the minimum values and resulted pixel will
have the sum of those minimum values. If a group has one value, keep the
same.
I am stuck to start with this process. I am assuming I need to convert the
brick into dataframe and do this.
Can any one help me in giving a lead on how to go about it?
Really appreciate any help.
Sorry if the explanation is confusing.
Regards
Sajid
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] problems in shortestPath()?
Hello Melanie, It's really hard to tell what your problem is with no error message or reproducible example. However here some info that may help and some trick to help you figure out the problem First, gdistance build a transition matrix from a CONDUCTANCE matrix which is not a COST matrix like ArcGIS. Therefore you initial raster should be 0 for land (no conductance) and 1 for water (conductance) and not 1 for water and 999 for land. Second, I recently had trouble with shortestPath() when the shortest path involve to much of a detour. So you should try first with too simple points where is not too hard to find the shortest path. Plot your transition matrix transi with plot(raster(transi)) and add your point to it with points(x,y), it will help you see if your trying to do something that make sense. Don't forget the geoCorrection(), it was necessary for me. Good luck, Bastien Ferland-Raymond, M.Sc. Stat., M.Sc. Biol. Division des orientations et projets spéciaux Direction des inventaires forestiers Ministère des Ressources naturelles et de la Faune [[alternative HTML version deleted]] ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] redundancy of coordinates of pixels
Hi everybody, I try to use the advice from Roger Bivand to get the redundant coordinates of my pixels. In the r help I found this example code: nn2(data, query, k=min(10,nrow(data)),treetype=c(kd,bd), searchtype=c(standard,priority,radius),radius=0.0,eps=0.0) I tried to adapt it according my aim: redon-read.table(C:\\Users\\Documents\\RepPixel.txt,sep=,dec=,,header=TRUE) nn2(redon,query,k=min(10,nrow(redon)),treetype=c(kd,bd), searchtype=c(priority),eps0.0) The result is not good. Can you help me solve this problem Thank you in advance KOMINE -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/redundancy-of-coordinates-of-pixels-tp7169408p7355090.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] redundancy of coordinates of pixels
Am 08.03.2012 15:48, schrieb Komine: Hi everybody, I try to use the advice from Roger Bivand to get the redundant coordinates of my pixels. In the r help I found this example code: nn2(data, query, k=min(10,nrow(data)),treetype=c(kd,bd), searchtype=c(standard,priority,radius),radius=0.0,eps=0.0) I tried to adapt it according my aim: redon-read.table(C:\\Users\\Documents\\RepPixel.txt,sep=,dec=,,header=TRUE) nn2(redon,query,k=min(10,nrow(redon)),treetype=c(kd,bd), searchtype=c(priority),eps0.0) The result is not good. Can you help me solve this problem Which problem? Maybe you can help us help you by first reading the posting guide: http://www.r-project.org/posting-guide.html HTH Tom Thank you in advance KOMINE -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/redundancy-of-coordinates-of-pixels-tp7169408p7355090.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Technische Universität München Department für Pflanzenwissenschaften Lehrstuhl für Grünlandlehre Alte Akademie 12 85350 Freising / Germany Phone: ++49 (0)8161 715324 Fax: ++49 (0)8161 713243 email: tom.gottfr...@wzw.tum.de http://www.wzw.tum.de/gruenland ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] redundancy of coordinates of pixels
Thank Tom for your observation.You are right, my question was not clear. When I run my code. I have this error message that I translated from french to english: Error in nrow (query): object'query' not found -- View this message in context: http://r-sig-geo.2731867.n2.nabble.com/redundancy-of-coordinates-of-pixels-tp7169408p7355260.html Sent from the R-sig-geo mailing list archive at Nabble.com. ___ 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
Mauricio Zambrano-Bigiarini mauricio.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,
Oscar.
--
Oscar Perpiñán Lamigueiro
Dpto. Ingeniería Eléctrica
EUITI-UPM
http://procomun.wordpress.com
___
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
[R-sig-Geo] GWR: set initial bandwidth value for gwr.sel()
Is there a way to set the initial bandwidth for gwr.sel() when doing geographically weighted regression? It looks like gwr.sel() calls optimize() which accepts a min/max range for parameters but not an initial value. I'd rather not constrain the bandwidth. However, I often send different subsets of the same dataset to gwr.sel and would like to speed up the optimization process by passing a value from the last call to the function. Solutions? thanks, Burton Shank ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] GWR: set initial bandwidth value for gwr.sel()
On Thu, 8 Mar 2012, Burton Shank wrote: Is there a way to set the initial bandwidth for gwr.sel() when doing geographically weighted regression? It looks like gwr.sel() calls optimize() which accepts a min/max range for parameters but not an initial value. I'd rather not constrain the bandwidth. However, I often send different subsets of the same dataset to gwr.sel and would like to speed up the optimization process by passing a value from the last call to the function. Solutions? Just set the bandwidth manually. gwr.sel() does use optimize(), which does not take a start value (optim() does, but does not work well for line search). Hope this clarifies, Roger thanks, Burton Shank ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Roger Bivand Department of Economics, NHH Norwegian School of Economics, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: roger.biv...@nhh.no ___ R-sig-Geo mailing list R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Error in function gbif package dismo
Dear list members,
I am trying to get the records for 1 species in Costa Rica using the gbif
function and I created an extent for the region (see extent object),
however, there is an error message when I run the gbif function:
e = extent(c(-9598065, -9199982, 901214, 1267766))
a=gbif(Centurio, senex, geo = T, ext = e)
Error in if (n == 0) { : argumento tiene longitud cero
What the error mean?
Best,
Manuel
--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.ac.cr
mspinol...@gmail.com
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río https://sites.google.com/site/lobitoderio/
Institutional website: ICOMVIS http://www.icomvis.una.ac.cr/
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[R-sig-Geo] Sum value of points that intersect a subset of a network
Hello list,
I would like to calculate the Getis-Ord Gi* clustering
statistic for points (mileposts) on a network (roads) and I donât think there
is a function in one of the spatial packages to do so. Hence I decided to code
it myself, but Iâm having some difficulties. Here is the general framework I
was hoping to implement:
1)     For each focal point i =1,2,â¦, N, grow a subset of the network
within distance k
2)Â Â Â Â Â Determine which points intersect the resultant
network subset and sum the values of all these points plus the focal point i
3)Â Â Â Â Â Calculate Gi* by dividing âneighborhoodâ sum by
study extent sum and redefine as a normal standard deviate as per Getis and
Ordâs
procedures (Getis and Ord 1992).
I was envisioning writing this code as a loop incrementing i  from 1 to N,
calculating the sum of values within each network âneighborhoodâ and
attributing that to each point. Then the rest are just dataframe calculations.
However, I got stuck in my trial and error and I have 2 questions please in
regards to implementing the above procedures that Iâm hoping other list
members
can help with:
First, I started with bullet number 1 creating the network âneighborhoodâ
using the lineardisc function in spatstat. For example the following code shows
the network subsets within distance k = 0.2, for focal points 2 and 5:
library(spatstat)
data(simplenet)
set.seed(100)
X - rpoislpp(5, simplenet)
xp - ppp(X$data$x,X$data$y)
plot(simplenet)
points(xp,col=c(blue,red,rep(blue,2),red,rep(blue,7)),pch=16,cex=1)
text(X$data$x,X$data$y,adj=0.01,cex=1.5)
lineardisc(simplenet,xp[2],0.2,cols=c(transparent,
red,transparent))
lineardisc(simplenet,xp[5],0.2,cols=c(transparent,
red,transparent))
Â
IÂ donât know how to proceed from here in terms of determining
which points intersect each network âneighborhoodâ so I can sum them. In
this
example the values of pts 2, 3, 10, and 11 will be summed and attributed to
focal
pt 2 and values of pts 5, 6, and 8 will be summed and attributed to focal pt
5. I tried to search for functions that could 1)
convert the output of lineardisc into a spatial lines object and 2) intersect
spatial
points with lines, but with no success⦠I would welcome any suggestions on
completing such tasks.
Second, my road network is an ESRI shp file which I
successfully converted into a psp object using maptools functions. However, I
donât know how to extract an edge list so I can convert it to a linnet object
using the linnet{spatstat} function. I would welcome any suggestions here as
well.
Thanks in advance for your help and apologies if Iâm missing
something obviousâ¦
Sharon
Sharon Baruch-Mordo
Ph.D. Candidate in Ecology
Dept. of Fish, Wildlife, and Conservation Biology
Colorado State University
Fort Collins, CO 80523-1474
Phone: 970-491-7020Â
Fax: 970-491-5091
E-mail: sharo...@warnercnr.colostate.edu
http://www.warnercnr.colostate.edu/research/bear
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] Error in function gbif package dismo
Manuel,
If you use an extent, you need to use longitude/latitude coordinates. This
is not a reasonable extent for this function:
e = extent(c(-9598065, -9199982, 901214, 1267766))
Robert
On Thu, Mar 8, 2012 at 3:06 PM, Manuel Spínola mspinol...@gmail.com wrote:
Dear list members,
I am trying to get the records for 1 species in Costa Rica using the gbif
function and I created an extent for the region (see extent object),
however, there is an error message when I run the gbif function:
e = extent(c(-9598065, -9199982, 901214, 1267766))
a=gbif(Centurio, senex, geo = T, ext = e)
Error in if (n == 0) { : argumento tiene longitud cero
What the error mean?
Best,
Manuel
--
*Manuel Spínola, Ph.D.*
Instituto Internacional en Conservación y Manejo de Vida Silvestre
Universidad Nacional
Apartado 1350-3000
Heredia
COSTA RICA
mspin...@una.ac.cr
mspinol...@gmail.com
Teléfono: (506) 2277-3598
Fax: (506) 2237-7036
Personal website: Lobito de río
https://sites.google.com/site/lobitoderio/
Institutional website: ICOMVIS http://www.icomvis.una.ac.cr/
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Re: [R-sig-Geo] writeRaster and coord ref problem
Johannes, Yes, that smells like something is wrong with your GDAL/PROJ4. Can you try this: library(rgdal) crs - +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0 CRS(crs) For me, it returns: CRS arguments: +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0 Robert On Wed, Mar 7, 2012 at 2:29 AM, Johannes Radinger jradin...@gmx.at wrote: Hi again, I also tried your example and I can reproduce the problem in contrast to you. After opening a fresh session, here is what I get: library(raster) Lade notiges Paket: sp raster version 1.9-70 (27-February-2012) crs - +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0 r - raster(nc=10, nr=10, crs=crs) r[] - 1:ncell(r) b - brick(r,r,r) x - writeRaster(b, 'datum.tif', overwrite=TRUE) Geospatial Data Abstraction Library extensions to R successfully loaded Loaded GDAL runtime: GDAL 1.9.0, released 2011/12/29 Path to GDAL shared files: /Users/Johannes Radinger/Library/R/2.14/library/rgdal/gdal Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470] Path to PROJ.4 shared files: /Users/Johannes Radinger/Library/R/2.14/library/rgdal/proj x class : RasterBrick dimensions : 10, 10, 100, 3 (nrow, ncol, ncell, nlayers) resolution : 36, 18 (x, y) extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +ellps=bessel +towgs84=606,23,413,0,0,0,0 +units=m +no_defs values : /Users/Johannes Radinger/Documents/Eclipse_workspace/datum.tif min values : 1 1 1 max values : 100 100 100 sessionInfo() R version 2.14.2 (2012-02-29) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rgdal_0.7-8 raster_1.9-70 sp_0.9-95 rj_1.0.0-3 loaded via a namespace (and not attached): [1] grid_2.14.2lattice_0.20-0 rj.gd_1.0.0-1 tools_2.14.2 What might cause the problem? Is there any bug I should report? Maybe the problem is platform specific (I use a Mac) or has something to do with my version of GDAL/PROJ4? Any suggestions? In the meanwhile I'll try your approach with resetting the CRS: projection(x) - crs /Johannes Original-Nachricht Datum: Tue, 6 Mar 2012 09:20:09 -0800 Von: Robert J. Hijmans r.hijm...@gmail.com An: Johannes Radinger jradin...@gmx.at CC: r-sig-geo@r-project.org Betreff: Re: [R-sig-Geo] writeRaster and coord ref problem Johannes, it works for me: library(raster) crs - +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0 r - raster(nc=10, nr=10, crs=crs) r[] - 1:ncell(r) b - brick(r,r,r) x - writeRaster(b, 'datum.tif', overwrite=TRUE) x class : RasterBrick dimensions : 10, 10, 100, 3 (nrow, ncol, ncell, nlayers) resolution : 36, 18 (x, y) extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) coord. ref. : +proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=350 +y_0=0 +datum=potsdam +units=m +no_defs +ellps=bessel +towgs84=606.0,23.0,413.0 values : datum.tif min values : 1 1 1 max values : 100 100 100 If for some reason part of the proj.4 description is lost, you can set it again with something like projection(x) - crs sessionInfo() R version 2.14.1 (2011-12-22) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United States.1252LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C LC_TIME=English_United States.1252 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rgdal_0.7-5 raster_1.9-71 sp_0.9-93 loaded via a namespace (and not attached): [1] grid_2.14.1lattice_0.20-0 Best, Robert On Tue, Mar 6, 2012 at 6:36 AM, Johannes Radinger jradin...@gmx.at wrote: Hi, I need to load some data from the Worldclim Database and reproject them to some other, already existing rasters. Because I don't want to use the reprojected rasters from the memory I want to save the brick as a GeoTiff (like the other rasters) using writeRaster. But somehow writeRaster() changes the coord. ref. line and drops +datum=Potsdam. Here what I am doing: WC.bio - projectRaster(getData(worldclim, var=bio, res=0.5, lon=9.3, lat=54.5),second.map) WC.bio class : RasterBrick dimensions : 912, 840, 766080, 19 (nrow, ncol, ncell, nlayers) resolution : 50, 50 (x, y) extent : 3505300, 3547300, 6025800, 6071400 (xmin, xmax, ymin,
