Re: [R] German Map in package maps

2007-05-22 Thread Ray Brownrigg
On Wed, 23 May 2007, Maja Schröter wrote:
> Hello everybody,
>
> I'm an absolute newbe in R, so please be gentle to me.
>
> I am looking for a german map in the package maps by Becker and Wilks.
>
> After that I want to draw the German cities from world.cities of this
> package in the map and want to draw further a line from Berlin to Munich.
>
> Has anybody an idea?
>
There is no germany-only map in the maps databases, but the following might be 
close to what you want:

> library(mapdata)  # the world database in the maps package
# probably isn't high enough resolution
> map("worldHires", "Germany")
> map.cities(country="Germany", minpop=1e6, capital=1)
> lines(world.cities[world.cities$name %in% c("Berlin", "Munich") & 
world.cities$country=="Germany", c("long", "lat")])
>  
HTH
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] User defined grid on map

2007-04-11 Thread Ray Brownrigg
On Thu, 12 Apr 2007, Adriana Bejarano wrote:
> This may be a trivial question for some of you: Is there a way to add
> defined grid sizes on specific areas of a map?. I want to add 10kmx10km
> grids to all California coastal counties:
>
> "california,alameda", "california,contra costa", "california,del norte",
> "california,humboldt", "california,los angeles", "california,marin",
> "california,mendocino" ,"california,monterey", "california,napa",
> "california,orange", "california,riverside", "california,sacramento",
> "california,san benito", "california,san bernardino", "california,san
> diego" , "california,san francisco", "california,san joaquin",
> "california,san luis obispo", "california,san mateo", "california,santa
> barbara", "california,santa clara", "california,santa cruz",
> "california,siskiyou", "california,solano", "california,sonoma",
> "california,sutter", "california,trinity", "california,ventura",
> "california,yolo"
>
The mapproj package has a map.grid() function, but that only allows you to 
specify latitude- and longitude-based grids.  I don't think it would be too 
difficult to roll your own using standard trigonometry if you weren't too 
worried about the non-sphericity of the earth.

Otherwise I guess there would be something 'out there' to provide you with 
some conversion between lat/long and distance from a fixed point.

HTH,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] How can i do the same thing in the China map?

2007-03-16 Thread Ray Brownrigg
usstata wrote:
> The "maps" package has a function called "match.map", which is for map 
> coloring .
> Its example is followed:
> 
> # filled map showing Republican vote in 1900
> # (figure 6 in the reference)
> data(state, package = "datasets")
> data(votes.repub)
> state.to.map <- match.map("state", state.name)
> x <- votes.repub[state.to.map, "1900"]
> gray.colors <- function(n) gray(rev(0:(n - 1))/n)
> color <- gray.colors(100)[floor(x)]
> map("state", fill = TRUE, col = color); map("state", add = TRUE)
> 
> I want to do the same thing in the China map, but I can't find the Provinces 
> name of China.
> Who can help me ?
> 
nobody
> 
> 
> a rookie
> 
Tell us who you are, and you may get a more substantial reply.

Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] circle fill problem

2007-02-08 Thread Ray Brownrigg
I suspect the number is:
trunc((x + 3*r)/(2*r)) * trunc((y + 3*r)/(2*r)) + trunc((x +
2*r)/(2*r)) * trunc((y + 2*r)/(2*r))
where x and y are the dimensions of the rectangle and r is the radius of the 
circle.

I have some code which graphs the overlapping circles if you want visually to 
check the algorithm (but although it is only a page, I won't post it here).

HTH
Ray Brownrigg

On Thursday 08 February 2007 22:39, MINI GHOSH wrote:
> Dear Ingmar and Robin,
>
> Thanks for you suggestions. I will see to it.
>
> Regards,
> Mini
>
> --- Ingmar Visser <[EMAIL PROTECTED]> wrote:
> > Robin & Mini,
> > For those interested, googling for the 'orange
> > packing problem' as it
> > is known, or more officially the sphere packing
> > problems gives you
> > quite a few hits on these and similar problems.
> > So at least the 3-d case the problem has been solved
> > (I imagine the
> > problem is easier in 2-d ...)
> > hth, Ingmar
> >
> > On 8 Feb 2007, at 09:52, Robin Hankin wrote:
> > > Mini
> > >
> > > This is a hard problem in general.
> > >
> > > Recreational mathematics has wrestled with
> > > this and similar problems over the years; the
> > > general field is the "set cover problem" but
> > > in your case the sets are uncountably infinite
> > > (and there are uncountably many of them).
> > >
> > > I would be surprised if your problem were not NP
> >
> > complete.
> >
> > > HTH
> > >
> > >
> > > Robin
> > >
> > > On 8 Feb 2007, at 05:15, MINI GHOSH wrote:
> > >> Dear R user,
> > >>
> > >> I want to know is there a way to find the minimum
> > >> number of circles (of given radius) required to
> >
> > fill a
> >
> > >> given area (say rectangular) where overlapping of
> > >> circles is allowed.
> > >>
> > >> Thanks,
> > >> Regards,
> > >> Mini Ghosh
> > >>
> > >> __
> > >> R-help@stat.math.ethz.ch mailing list
> > >> https://stat.ethz.ch/mailman/listinfo/r-help
> > >> PLEASE do read the posting guide
> >
> > http://www.R-project.org/posting-
> >
> > >> guide.html
> > >> and provide commented, minimal, self-contained,
> >
> > reproducible code.
> >
> > > --
> > > Robin Hankin
> > > Uncertainty Analyst
> > > National Oceanography Centre, Southampton
> > > European Way, Southampton SO14 3ZH, UK
> > >   tel  023-8059-7743
> > >
> > > __
> > > R-help@stat.math.ethz.ch mailing list
> > > https://stat.ethz.ch/mailman/listinfo/r-help
> > > PLEASE do read the posting guide
> >
> > http://www.R-project.org/posting-
> >
> > > guide.html
> > > and provide commented, minimal, self-contained,
> >
> > reproducible code.
>
> __
> R-help@stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html and provide commented, minimal,
> self-contained, reproducible code.

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Cannot scale my map

2006-12-11 Thread Ray Brownrigg
On Tuesday 12 December 2006 09:07, Mihalis Tsoukalos wrote:
> Dear list,
>
> I have the following problem:
>
> I want to plot my data and display a map at the same graph. The
> problem is that my map is not properly sized (is very small). It
> appears at the bottom-left side of the output.
> I am using the following R commands:
>
> {
> plot(Time, col="lightgrey")
> title("Connections per Minute")
> map("world", fill=TRUE, lwd=0.25, col="green", add=TRUE)
> box()
> }
>
> Please note that Time contains data in the following format:
> 00:00
> 00:01
> ...
> 22:04, etc.
>
> Can you please help me?
>
Well, the only way a map(..., add=TRUE) will work is if the original scale is 
in the order of latitude and longitude in degrees.  It is not clear exactly 
how the data you are plotting relates to the map you want.

If you just want the map as a 'background', then you will need to rescale one 
or other of the plots.

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] indexing question

2006-11-13 Thread Ray Brownrigg
On Tuesday 14 November 2006 09:28, Leeds, Mark (IED) wrote:
> thanks beilton but that won't work. A diff will also include 61-41 etc
> and I don't want to include those.
>
> I'm working on using lapply or sapply with a seq along 31, 61, etc.
> I'll let you know if it works.
>

Try looking at:
dim(idx) <- c(3, length(idx)/3)
then:
tmp[idx[2, ]] - temp[idx[1, ]]
temp[idx[3, ]] - temp[idx[2, ]]

HTH
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] RE : avoiding a loop: "cumsum-like"

2006-11-07 Thread Ray Brownrigg
Well, I do have a solution which works for the data set you provide, but 
possibly not in a more general case.

Firstly, tidying up your code, but using essentially the same looping 
algorithm, can provide a speed improvement of approximately 3:1.
Here is a first attempt:
mycode1 <- function(tab) {
  len <- diff(range(tab$Date)) + 1
  res <- numeric(len)
  val <- 0
  for (i in 1:len)
  {
if (is.na(tab$posit.lat[i]))
{
  val <- val + tab$x.jour[i]
}
else
{
  if (res[tab$posit.lat[i]] < 30)
  {
val <- val + tab$x.jour[i]
  }
  else
  {
val <- val + tab$x.jour[i] + 0.8*res[tab$posit.lat[i]]
  }
}
  res[i] <- val
  }
  return(res)
}

Then using a cumsum()-based algorithm can provide an overall 10:1 speed 
improvement:

mycode2 <- function(tab) {
  res0 <- cumsum(tab$x.jour)
  res1 <- cumsum(ifelse(is.na(tab$posit.lat), 0, 0.8*
(res0[tab$posit.lat] >= 30) * res0[tab$posit.lat]))
  res2 <- cumsum(ifelse(is.na(tab$posit.lat), 0, 0.8*res1[tab$posit.lat]))
  return(res0 + res1 + res2)
}

The condition is that:
res[tab$posit.lat[tab$posit.lat[tab$posit.lat[length(tab$posit.lat) < 30
where tab is the data and res is the result.  [There is also an implicit 
assumption that the result is monotonic.]

HTH
Ray Brownrigg

GOUACHE David wrote:
> Thanks Petr for taking a stab at it.
> I have yet to figure out a way to do it, but if I do I'll post it.
> Cheers
>
> David
>
> -Message d'origine-
> De : Petr Pikal [mailto:[EMAIL PROTECTED] 
> Envoyé : vendredi 3 novembre 2006 09:05
> À : GOUACHE David; r-help@stat.math.ethz.ch
> Objet : Re: [R] avoiding a loop: "cumsum-like"
>
> Hi
>
> I have not seen any answer yet so I wil try (partly).
>
> I believe that the loop can be vectorised but I am a little bit lost 
> in your fors and ifs. I found that first part of res is same as 
> cumsum(tab$x.jour) until about 81st value. However I did not decipher 
> how to compute the remaining part. I tried to add 
> cumsum(tab$posit.lat) (after changing NA to 0) what is not correct.
>
> Probably some combination of logical operation and summing can do 
> what you want. I thought that something like
> ((cumsum(tab$posit.lat)*0.8)*(cumsum(tab$x.jour)>30)+cumsum(tab$x.jour
> ))
>
> can do it but the result is defferent from your computation.
> Not much of help, but maybe you can do better with above suggestion.
>
> Petr
>
>
>
> On 2 Nov 2006 at 11:15, GOUACHE David wrote:
>
> Date sent:Thu, 2 Nov 2006 11:15:49 +0100
> From: "GOUACHE David" <[EMAIL PROTECTED]>
> To:   
> Subject:  [R] avoiding a loop: "cumsum-like"
>
>   
>> Hello Rhelpers,
>>
>> I need to run the following loop over a large number of data-sets, and
>> was wondering if it could somehow be vectorized. It's more or less a
>> cumulative sum, but slightly more complex. Here's the code, and an
>> example dataset (called tab in my code) follows. Thanks in advance for
>> any suggestions!
>>
>> res<-0
>> for (i in min(tab$Date):max(tab$Date))
>> {
>>  if (is.na(tab$posit.lat[tab$Date==i])==T)
>>  {
>>   res<-c(res,res[length(res)]+tab$x.jour[tab$Date==i])
>>  }
>>  else
>>  {
>>   if (res[tab$posit.lat[tab$Date==i]+1]<30)
>>   {
>>res<-c(res,res[length(res)]+tab$x.jour[tab$Date==i])
>>   }
>>   else
>>   {
>>res<-c(res,res[length(res)]+tab$x.jour[tab$Date==i]+0.8*res[tab$pos
>>it.lat[tab$Date==i]+1])
>>   }
>>  }
>> }
>> res[-1]
>>
>>
>> Date x.jour  posit.lat
>> 358040   NA
>> 358050   NA
>> 358060   NA
>> 358070   NA
>> 358080   NA
>> 358092.97338883  NA
>> 358102.796389915 NA
>> 358110   NA
>> 358120   NA
>> 358131.000711886 NA
>> 358140.894422571 NA
>> 358150   NA
>> 358160   NA
>> 358170   NA
>> 358180   NA
>> 358190   NA
>> 358200   NA
>> 358210   NA
>> 358220   NA
>> 358230   NA
>> 358240   NA
>> 358250   NA
>> 358260   NA
>> 358270   NA
>> 358280   NA
>> 358290   NA
>> 358300   NA
>> 358310   NA
>> 358320   NA
>> 358330   NA
>> 358340   NA
>> 35

Re: [R] avoiding a loop

2006-10-26 Thread Ray Brownrigg
On Wednesday 25 October 2006 07:36, Leeds, Mark (IED) wrote:
> I think I asked a similar question 3 years ago to the  Splus list and I
> think the answer was no or noone answered so noone should spend more
> than 5 minutes on this
> because it could definitely be a waste of time.
>
> My question is whether the function below can be rewritten without a for
> loop. apply is fine if it can be done that way but i doubt it. I call it
> a lot and would
> prefer to not loop.
>
> #---
> --
>
> constructLt<-function(invector) {
>
> outvector<-invector
>
>  for ( i in 2:length(invector) ) {
>  if ( invector[i] < 1 ) {
>   outvector[i]<-invector[i]*outvector[i-1]
>  }
> }
>
> return(outvector)
>
> }
>
Depending on the nature of your data, there is a faster way.  It still 
involves looping, but not over the entire vector.

Try the following:
constructLt <- function(invector) {
  outvector <- invector
  cs <- cumsum(rle(invector < 1)$lengths)
  if (invector[1] < 1)
cs <- c(1, cs)
  for (i in 0:(length(cs)%/%2 - 1)){
starti <- cs[2*i + 1]
stopi <- cs[starti + 1]
outvector[starti:stopi] <- cumprod(invector[starti:stopi])
  }
  return(outvector)
}

It is in the order of 3 times as fast for random vectors of considerable 
length (> 1000).  For random vectors of length 50 it is about the same speed 
as the full looping algorithm.  However if the data is such that there are 
longer runs (than N(1, 1)), then you might expect a better speedup.

HTH
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Possible bug?

2006-10-04 Thread Ray Brownrigg
On Thursday 05 October 2006 00:14, ronggui wrote:
> Yes, the installer missed some file(s) in RHOME/bin.
> When I choose "Message translations", the /bin directory contains the
> following files.
>
 -- snip --

I meant to mention that when I used the installer I didn't have to 
choose "Message translations", it was ticked by default.

HTH,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Possible bug?

2006-10-03 Thread Ray Brownrigg
On Wednesday 04 October 2006 14:17, Duncan Murdoch wrote:
> On 10/3/2006 9:00 PM, ronggui wrote:
> > This morning I downloaded R-2.4.0 and install in under Windows. I
> > customized the installation and choose "Message translations",but I could
> > not launch Rgui.exe successfully( Rterm.exe worked fine). If I did't
> > choose "Message translations", Rgui.exe worked fine.
> >
> > "Message translations" is Simplified Chinese.
>
> Did you test any of the betas or release candidates?  I run an English
> language version of Windows, and I don't even get offered the chance to
> install in Simplified Chinese, so I certainly didn't test this.
>
> Duncan Murdoch
>

> 7*runif(1)
[1] 2.897160

Well, I have a (self compiled) R-2.4.0 beta (2006-09-24 r39497) which doesn't 
exhibit the behaviour in the same environment (Windows XP Professional 
English version but with Simplified Chinese set as the locale) that the 
released R-2.4.0-win32.exe does.  [Perhaps my version didn't set up the 
equivalent of  "Message translations" so this may be a red herring, but it 
does output the startup message in Chinese, and error messages are also 
output in Chinese.]

What happens with the release version is an error message:
"R for Windows GUI front-end has encountered a problem and needs to close.  We 
are sorry for the inconvenience."  Then it offers to "Please tell Microsoft 
about the problem", and if you "click here"  the error signature has:
AppName: rgui.exe   AppVer: 2.40.39566.0ModName: r.dll
ModVer: 2.40.39566.0Offset: 000f22b3

There is a lot more information as well, but I don't know how much is 
relevant.

By "Simplified Chinese set as the locale" I mean 'Control Panel/Regional and 
Language Options' with Chinese (PRC) set in 'Regional Options/Standards and 
Formats' and 'Advanced/Language for non-Unicode Programs'.

[I don't rely on the Chinese environment, so this doesn't materially affect 
me, but I may be able to help with diagnosis.]

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] Speeding indexing and sub-sectioning of 3d array

2006-08-13 Thread Ray Brownrigg
> Date: Thu, 10 Aug 2006 14:34:27 -0400
> From: "Swidan, Firas" <[EMAIL PROTECTED]>
> 
> Hi Patrick,
> 
> Thanks for the help. The function I listed is just an example. I isolated
> and kept only the problematic part in my code for clarity sake. I ended up
> implementing the functionality in C and now it takes 22 seconds to calculate
> the objective.
> 
> Best regards,
> Firas.
> 
Interestingly, I was able to develop an algorithm in R that achieves the
same order-of-magnitude speedup as your C code, but at the expense of
greater memory requirements.  However it only works if the function you
are using is really is mean() [your code labels use Median].  It does
this by making use of cumsum() and logical indexing, working with sums
of values rather than calculationg the means and then dividing by the
numbers of values in the hypercube at the end.

If you want to try coding this algorithm in C for even greater
performance improvement (or for interest only), let me know.  I suspect
it will be difficult to code in C because of the vectorisation it takes
advantage of.

In the output below, cK3d() is your algorithm (slightly adjusted to
cover the whole matrix and to return something), and cK3dme is my
equivalent, running on a Pentium IV 3.2GHz, NetBSD system with 1GB
memory.

Regards,
Ray Brownrigg

> x <- rnorm(245*175*150)
> dim(x) <- c(245, 175, 150)
> unix.time(yme <- cK3dme(x, 3))
[1] 13.870  1.690 15.813  0.000  0.000
> unix.time(y <- cK3d(x, 3))
[1] 500.206   0.035 505.738   0.000   0.000
> all.equal(y, yme)
[1] TRUE
> 


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.


Re: [R] add city and point in the map

2006-04-27 Thread Ray Brownrigg
> From: YIHSU CHEN <[EMAIL PROTECTED]>
> Date: Thu, 27 Apr 2006 10:50:49 -0400
> 
> I'm trying to use packages "maps" and "mapdata" (see blow) to display
> the research resutls on map (Mid-Atlantic region).  In particular, I
> need to mark a number of points in the map by giving their latitude and
> longitude information.  For instance, I would like to mark a point on
> (long, lat) =(75.56027, 39.09271).  Also, I need to mark several cities
> that I suspect might not in the "world.cities" database such as
> Baltimore, etc.  Presumably, I can do the same thing by finding out the
> lat and long infor of Baltimroe and add text next to the point I mark.
> 
> Thanks a billion
> 
> Yihsu
> 
> 
> library(maps)
> library(mapdata)
> state.list <- c("maryland","new jersey","west virginia", 
> "virginia","ohio","delaware","pennsylvania","kentucky","indiana","north 
> carlina")
> map("state",state.list)
> map.cities(x=world.cities, Country="US",minpop=10)
> 
Several things:
1) the US is west of longitude 0, so you need longitude -75.56027
2) map.cities uses x=world.cities by default, so you don't need to
specify it in your call, but if you do specify it, you need to precede
it with a call to:
data(world.cities)
[I guess this is a bug, perhaps just at the documentation level]
3) map.cities takes country= as an argument, not Country=
4) The US is identified in the maps package as USA, so you would use
country="USA" (but you don't need country= anyway, beacause your map is
only of the USA, so there are no cities to exclude).
5) Baltimore is in world.cities, try map.cities(minpop=50, label=T)
6) If you really want North Carolina, you need to spell it correctly

So what you need is:
library(maps)
# library(mapdata)  # not needed for this
state.list <- c("maryland", "new jersey", "west virginia", "virginia",
  "ohio", "delaware", "pennsylvania", "kentucky", "indiana", "north carolina")
map("state", state.list)
map.cities(minpop=10)

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] matching a given sd range

2006-04-01 Thread Ray Brownrigg
> Date: Fri, 31 Mar 2006 14:44:05 -0800 (PST)
> From: "Fred J." <[EMAIL PROTECTED]>
>  
>  given a numeric array (a sequence of reals), I am interested in
>  finding the subsets of sequences (each with start and end index) which match 
> a given sd range.
>  
>  I read the docs on match and which and the "see also" but could not come up 
> with a way. I could loop with a stepping window over the sequence but that 
> would be limited to a fixed size window, I guess I as well loop through all 
> window sizes from 1:length(data) but that is not elegant.  any hints towards 
> a solution?
>  
Presumably your "array" is a vector, in which case you can use
vectorisation to eliminate one of the loops.  The following function
produces a matrix in upper triangular form for which res[i, j] is the
variance of the sequence dat[i:j] (where dat is your "array").  It is
(approximately) two orders of magnitude faster than the double loop,
but doesn't produce exactly the same numbers because of numerical
rounding because it uses the formula:  variance = (sum(xi^2) -
n*xbar^2)/(n - 1)

Now all you have to do is test this matrix against the range of
variances you want to isolate, and use the matching indices to provide
you with the endpoints of your subsets.

Hope this helps,
Ray Brownrigg

allvar <- function(dat) {
  len <- length(dat)
  res <- numeric((len - 1)*len/2)
  mym <- dat[1]
  myv <- 0
  for (j in 2:len) {
oldm <- mym
mym <- (((j - 1):1)*mym + dat[j])/(j:2)
myv <- (((j - 2):0)*myv + ((j - 1):1)*oldm^2 + dat[j]^2 - (j:2)*mym^2)/
  ((j - 1):1)
res[((j - 2)*(j - 1)/2 + 1):((j - 1)*j/2)] <- myv
mym <-c(mym, dat[j])
myv <- c(myv, 0)
  }
return(res)
}
# Now test it
> n <- 500; dat <- runif(n)
> unix.time({res2 <- matrix(0, n, n); res2[upper.tri(res2)] <- allvar(dat)})
[1] 0.64 0.00 0.64 0.00 0.00
> unix.time({res1 <- matrix(0, n, n); for (j in 2:n) for (i in 1:(j - 1)) 
> res1[i, j] <- var(dat[i:j])})
[1] 72.64 13.74 86.38  0.00  0.00
> max(abs(res1 - res2))
[1] 1.249001e-15
> 


__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Plotting shapefiles on existing maps

2006-03-29 Thread Ray Brownrigg
> From: [EMAIL PROTECTED]
> 
> I have plotted a map of the Barents Sea and surrounding coastline using:
> 
> map('worldHires',ylim=c(50,85),xlim=c(5,65),fill=T,resolution=0)
> map.axes()
> map.scale(x=30,metric=T)
> 
> Next, I imported a shapefile with depth contours for the sea:
> 
> contours<-read.shape("D://My Documents/BarentsSea.shp",dbf.data=T)
> 
> (This is in mercator projection).
> 
> Despite extensive searches of the help files and R site, I cannot find a
> way to plot the contours onto the map.  Does anyone have any suggestions?
> 
It is not clear to me exactly what you are having difficulty with, but
if you are able to produce a contour plot from that data, then note that
both contour() and map() have an "add" argument, so as long as the axis
limits are conformant (i.e. both longitude/latitude) you can add one
plot to the other.

Perhaps the trick is in converting your contours data into
longitude/latitude form for the x and y components.

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] bug in map('world') ?

2006-03-28 Thread Ray Brownrigg
> From: Joerg van den Hoff <[EMAIL PROTECTED]>
> 
> hi,
> 
> did'nt see anything in the archive:
> 
> map('world',pro='rectangular',para=0)
> 
> yields a strange artifact (horizontal bar) extending over the whole map 
> at a certain latitude range (approx 65 deg. north), whereas
> 
> map('world',pro='rectangular',para=180)
> 
> (which should be the same) does not show the artifact.
> 
> the artifact shows up in other projections as well, e.g.
> 
> map('world',pro='bonne',para=45)
> 
> 
> which seems to show that the problem might be in the data base of the 
> map (i.e. polygon definition)??
> 
> 
> any ideas?
> 
The short answer is RTFM, and now that I have again :-), I note that
there is a wrap= option.  Setting that to TRUE seems to have the
desired effect.

Ray

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] bug in map('world') ?

2006-03-08 Thread Ray Brownrigg
I seem to recall I came across something similar recently (well,
relatively recently, perhaps a couple of years ago).  The problem is
related to how the code handles polygons that are split across the map
boundaries.  As I recall, the fix was to modify the polygons that
straddle the date line.  I'll have to delve into this again, but don't
hold your breath, sorry.

Ray Brownrigg

> From: Joerg van den Hoff <[EMAIL PROTECTED]>
> 
> hi,
> 
> did'nt see anything in the archive:
> 
> map('world',pro='rectangular',para=0)
> 
> yields a strange artifact (horizontal bar) extending over the whole map 
> at a certain latitude range (approx 65 deg. north), whereas
> 
> map('world',pro='rectangular',para=180)
> 
> (which should be the same) does not show the artifact.
> 
> the artifact shows up in other projections as well, e.g.
> 
> map('world',pro='bonne',para=45)
> 
> 
> which seems to show that the problem might be in the data base of the 
> map (i.e. polygon definition)??
> 
> 
> any ideas?
> 
> regards,
> 
> joerg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Installing packages without clicking

2006-02-14 Thread Ray Brownrigg
> I need to install several (too many) packages from local *.zip
> files.
> 
Are you aware of being able to select multiple .zip files at one time
using SHIFT-click or CTRL-click?

> is there any form to do it  without clicking?
> 
> I'm asking for a R code that allow me perform these task.
> 
You could just do what the clicking does, i.e. call the function
install.packages() with appropriate parameters.  See
utils:::menuInstallLocal

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] expand.grid without expanding

2006-02-08 Thread Ray Brownrigg
> From: =?iso-8859-1?q?Lu=EDs_Torgo?= <[EMAIL PROTECTED]>
> Date: Wed, 8 Feb 2006 18:08:40 +
> 
> Dear list,
> I've recently came across a problem that I think I've solved and that I 
> wanted 
> to share with you for two reasons:
> - Maybe others come across the same problem.
> - Maybe someone has a much simpler solution that wants to share with me ;-)
> 
> The problem is as follows: expand.grid() allows you to generate a data.frame 
> with all combinations of a set of values, e.g.:
> > expand.grid(par1=-1:1,par2=c('a','b'))
>   par1 par2
> 1   -1a
> 20a
> 31a
> 4   -1b
> 50b
> 61b
> 
> There is nothing wrong with this nice function except when you have too many 
> combinations to fit in your computer memory, and that was my problem: I 
> wanted to do something for each combination of a set of variants, but this 
> set was to large for storing in memory in a data.frame generated by 
> expand.grid. A possible solution would be to have a set of nested for() 
> cycles but I preferred a solution that involved a single for() cycle going 
> from 1 to the number of combinations and then at each iteration having some 
> form of generating the combination "i". And this was the "real problem": how 
> to generate a function that picks the same style of arguments as 
> expand.grid() and provides me with the values corresponding to line "i" of 
> the data frame that would have been created bu expand.grid(). For instance, 
> if I wanted the line 4 of the above call to expand.grid() I should get the 
> same as doing:
> > expand.grid(par1=-1:1,par2=c('a','b'))[4,]
>   par1 par2
> 4   -1b
> 
> but obviously without having to use expand.grid() as that involves generating 
> a data frame that in my case wouldn't fit in the memory of my computer.
> 
> Now, the function I've created was the following:
> 
> getVariant <- function(id,vars) {
>   if (!is.list(vars)) stop('vars needs to be a list!')
>   nv <- length(vars)
>   lims <- sapply(vars,length)
>   if (id > prod(lims)) stop('id above the number of combinations!')
>   res <- vector("list",nv)
>   for(i in nv:2) {
> f <- prod(lims[1:(i-1)])
> res[[i]] <- vars[[i]][ceiling(id / f)]
> id <- id - (ceiling(id/f)-1)*f
>   }
>   res[[1]] <- vars[[1]][id]
>   names(res) <- names(vars)
>   res
> }
> --
> > expand.grid(par1=-1:1,par2=c('a','b'))[4,]
>   par1 par2
> 4   -1b
> > getVariant(4,list(par1=-1:1,par2=c('a','b')))
> $par1
> [1] -1
> 
> $par2
> [1] "b"
> 
> I would be glad to know if somebody came across the same problem and has a 
> better suggestion on how to solve this.
> 
A few minor improvements:
1) let id be a vector of indices
2) use %% and %/% instead of ceiling (perhaps debateable)
3) return a data frame as does expand.grid

So your function now looks like:

getVariant <- function(id, vars) {
  if (!is.list(vars)) stop('vars needs to be a list!')
  nv <- length(vars)
  lims <- sapply(vars, length)
  if (any(id > prod(lims))) stop('id above the number of combinations!')
  res <- vector("list", nv)
  for(i in nv:2) {
f <- prod(lims[1:(i-1)])
res[[i]] <- vars[[i]][(id - 1)%/%f + 1]
id <- (id - 1)%%f + 1
  }
  res[[1]] <- vars[[1]][id]
  names(res) <- names(vars)
  return(as.data.frame(res))
}

Now, for example, you get:

> expand.grid(par1=-1:1,par2=c('a','b'),par3=c('w','x','y','z'))[12:15,]
   par1 par2 par3
121bx
13   -1ay
140ay
151ay
> getVariant(12:15,list(par1=-1:1,par2=c('a','b'), par3=c('w','x','y','z')))
  par1 par2 par3
11bx
2   -1ay
30ay
41ay
>  

Note that you will run into trouble when the product of the lengths is
greater than the largest representable integer on your system.

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Number of replications of a term

2006-01-24 Thread Ray Brownrigg
There's an even faster one, which nobody seems to have mentioned yet:

rep(l <- rle(ids)$lengths, l)

Timing on my 2.8GHz NetBSD system shows:

> length(ids)
[1] 45150
> # Gabor:
> system.time(for (i in 1:100) ave(as.numeric(factor(ids)), ids, FUN =
length))
[1] 3.45 0.06 3.54 0.00 0.00
> # Barry (and others I think):
> system.time(for (i in 1:100) table(ids)[ids])
[1] 2.13 0.05 2.20 0.00 0.00
> Me:
> system.time(for (i in 1:100) rep(l <- rle(ids)$lengths, l))
[1] 1.60 0.00 1.62 0.00 0.00

Of course the difference between 21 milliseconds and 16 milliseconds is
not great, unless you are doing this a lot.

Ray Brownrigg

> From: Gabor Grothendieck <[EMAIL PROTECTED]>
> 
> Nice.  I timed it and its much faster than mine too.
> 
> On 1/24/06, Barry Rowlingson <[EMAIL PROTECTED]> wrote:
> > Laetitia Marisa wrote:
> > > Hello,
> > >
> > > Is there a simple and fast function that returns a vector of the number
> > > of replications for each object of a vector ?
> > > For example :
> > > I have a vector of IDs :
> > > ids <- c( "ID1", "ID2", "ID2", "ID3", "ID3","ID3", "ID5")
> > >
> > >  I want the function returns the following vector where each term is the
> > > number of replicates for the given id :
> > > c( 1, 2, 2, 3,3,3,1 )
> >
> > One-liner:
> >
> >  > table(ids)[ids]
> > ids
> > ID1 ID2 ID2 ID3 ID3 ID3 ID5
> >   1   2   2   3   3   3   1
> >
> >  'table(ids)' computes the counts, then the subscripting [ids] looks it
> > all up.
> >
> >  Now try it on your 40,000-long vector!
> >
> > Barry

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] elements in a matrix to a vector

2005-11-09 Thread Ray Brownrigg
> my apologies to everyone taking the time to answer this question.  i
> didn't explain myself very well.  
> 
> i have the indices of the elements of the matrix and the matrix won't
> necessarily have zeros everywhere else.  is there a way to use the
> indices to strip the elements i need?
> 
RTFM ?"[".

> x <- matrix(c(0, 2, 0, 0, 0, 4, 3, 0, 0), nrow=3)
> x
 [,1] [,2] [,3]
[1,]003
[2,]200
[3,]040
> ?"["
> xx <- matrix(c(1, 2, 3, 3, 1, 2), ncol=2)
> xx
 [,1] [,2]
[1,]13
[2,]21
[3,]32
> x[xx]
[1] 3 2 4
>   

HTH,
Ray

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Error in integrate

2005-10-06 Thread Ray Brownrigg
Murray Jorgensen <[EMAIL PROTECTED]> wrote:

> I'm using R 2.0.1 under Windows XP. I get the following message when I 
> run the code listed below. I don't seem to have any problems using the 
> function "slice" outside "integrate".
> 
> Error in integrate(slice, 0, Inf, , , , , , , a, b) * delta :
>  non-numeric argument to binary operator
> 
RTFM!  Integrate returns:
Value:

 A list of class '"integrate"' with components
 :
 :

You need:
A[j] <- integrate(slice,0,Inf,,,a,b)$value*delta

Ray

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] generating a vector from clusters of logicals

2005-09-08 Thread Ray Brownrigg
> From: Troels Ring <[EMAIL PROTECTED]>
> 
> I have a vector of clusters of TRUE and FALSE like 
> c(TRUE,TRUE,TRUE...,FALSE,FALSE,FALSE,TRUE,TRUE...) and want to make 
> that into a vector
> of c(1,1,1,1...2,2,2,2,.3,3,3,3) increasing the number assigned to each 
> cluster as they change.
> How would I do that ?
> 
How about:
> TF <- c(TRUE,TRUE,TRUE,FALSE,FALSE,FALSE,TRUE,TRUE,FALSE)
> rep(1:length(rlel <- rle(TF)$lengths), rlel)
[1] 1 1 1 2 2 2 3 3 4

HTH,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


[R] problem with Hershey fonts

2005-07-26 Thread Ray Brownrigg
This was reported to me by a colleague in China, so I may not be
reproducing exactly what they are seeing (which I suspect is rw2011), but
this is what I see:
> version
 _
platform i386--netbsdelf
arch i386
os   netbsdelf
system   i386, netbsdelf
status
major2
minor1.1
year 2005
month06
day  20
language R
> help(Hershey)
:
:
 If the 'vfont' argument to one of the text-drawing functions
 ('text', 'mtext', 'title', 'axis', and 'contour') is a character
 vector of length 2, Hershey vector fonts are used to render the
 text.
:
  The other useful escape sequences all begin with '\\'.  These
  are described below. Remember that backslashes have to be
  doubled in R character strings, so they need to be entered
  with _four_ backslashes.
:
:
> plot(runif(100))  # to get something on the screen
> text(0, 1, '\\Re', vfont=c('serif', 'plain')) # works
> title(main = '\\Re', vfont=c('serif','plain'))# doesn't work, and...
Warning message:
parameter "vfont" could not be set in high-level plot() function
> mtext('\\Re', 2, vfont=c('serif','plain'))# doesn't work, but...
Warning message:
Hershey fonts not yet implemented for mtext() in: mtext(text, side, line, 
outer, at, adj, padj, cex, col, font,
> axis(3, at=50, tick=F, labels='\\Re', vfont=c('serif','plain'))   # 
> doesn't work 

Now at least mtext() 'redeems' itself with its Warning, but title()
contradicts the help(Hershey), and axis() just gets it wrong by printing
the characters "\Re" (or "\\Re" when the string is specified as
'Re' as alluded to in the Hershey documentation - but somewhat
contradicted by the success of the text() call above).

Now perhaps the problem is just with the Hershey documentation (or my
reading of it), but I don't know what is supposed to work.  Or is there
some combination that I haven't tried?

Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Maps, Eastern Europe

2005-05-22 Thread Ray Brownrigg
> From: "Rau, Roland" <[EMAIL PROTECTED]>
> 
> I would like to employ a European map in a presentation. My idea was to
> use:
> library(mapdata)
> map("worldHires", c("Austria", "Switzerland", "Germany"))
> 
> where I include all countries from my analysis as a vector of character
> strings like in the example above. Unfortunately, I was unable to
> specify the Baltic States (Latvia, Lithuania, Estonia) or the Czech
> Republic. 
> Initially, I thought I misspelled the country names. To avoid this
> obvious error, I checked within 'mapdata_2.0-14.tar.gz' the file
> 'worldHires.names'. To my surprise, the spelling was not the source of
> the error - it was rather the case that those countries were not
> included.
> 
> Did some people experience similar problems? And if yes: how did you
> proceed?
> 
Unfortunately the world databases for the maps package are very old
(early 1990's) and have not been updated since, except for minor
typographical errors in country names.  Note that Czechoslovakia is
there (as is USSR!).  Also, it currently satisfies my needs, which is
why I don't have the resources to update it.

I suggest you start with the maptools package, but as I understand it,
you need to also find some appropriate shapefile data. [Other
respondents may improve upon this suggestion.]

Regards,
Ray Brownrigg
maps package maintainer

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] a question about box counting

2005-04-04 Thread Ray Brownrigg
I said:
> myfun <- function(x, y, ints) {
>   fx <- x %/% (1/ints)
>   fy <- y %/% (1/ints)
>   txy <- hist(fx + ints*fy+ 1, breaks=0:(ints*ints), plot=FALSE)$counts
>   dim(fxy) <- c(ints, ints)
^^^
>   return(txy)
> }
Of course it should be:
  dim(txy) <- c(ints, ints)
  ^^^

Sorry about that,
Ray

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] a question about box counting

2005-04-04 Thread Ray Brownrigg
> From: Deepayan Sarkar <[EMAIL PROTECTED]> Mon, 4 Apr 2005 13:52:48 -0500
> 
> On Monday 04 April 2005 13:22, Rajarshi Guha wrote:
> > Hi,
> >   I have a set of x,y data points and each data point lies between
> > (0,0) and (1,1). Of this set I have selected all those that lie in
> > the lower triangle (of the plot of these points).
> >
> > What I would like to do is to divide the region (0,0) to (1,1) into
> > cells of say, side = 0.01 and then count the number of cells that
> > contain a point.
> >
> > My first approach is to generate the coordinates of these cells and
> > then loop over the point list to see whether a point lies in a cell
> > or not.
> >
> > However this seems to be very inefficient esepcially since I will
> > have 1000's of points.
> >
> > Has anybody dealt with this type of problem and are there routines to
> > handle it?
> 
> A combination of cut and table/xtabs should do it, e.g.:
> 
> 
> x <- runif(3000)
> y <- runif(3000)
> 
> fx <- cut(x, breaks = seq(0, 1, length = 101))
> fy <- cut(y, breaks = seq(0, 1, length = 101))
> 
> txy <- xtabs(~ fx + fy)
> :

Another significantly faster way (but not generating row/column names)
is:
x <- runif(3000)
y <- runif(3000)
ints <- 100
myfun <- function(x, y, ints) {
  fx <- x %/% (1/ints)
  fy <- y %/% (1/ints)
  txy <- hist(fx + ints*fy+ 1, breaks=0:(ints*ints), plot=FALSE)$counts
  dim(fxy) <- c(ints, ints)
  return(txy)
}
myfun(x, y, ints)

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Recursive default argument reference

2005-01-25 Thread Ray Brownrigg
> From: boston knot <[EMAIL PROTECTED]>
> Date: Mon, 24 Jan 2005 20:50:53 -0800 (PST)

> "Recursive default argument reference" keeps appearing when I try to run a 
> haplo.score function in R for Windows. I'm new to using this program. Does 
> anyone know what this means? 
>  
I don't know exactly what it means, but I have seen it quite some time
ago when I was creating a package using delay().  It is easily
reproduced by:
> x <- delay({runtimeerror; 10})
> x
Error: Object "runtimeerror" not found
> x
Error: recursive default argument reference
> 

HTH,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] zero index and lazy evaluation in ifelse()

2005-01-13 Thread Ray Brownrigg
> Date: Thu, 13 Jan 2005 12:42:59 -0800
> From: [EMAIL PROTECTED]
> 
> I don't understand this behavior:
> 
> > a <- c(0, 1, 2, 3)
> > b <- c(1, 2, 3, 4)
> > ifelse (a == 0, 0, b[a])
> [1] 0 2 3 1
> 
> rather than the desired 0 1 2 3.  Thanks for any explanation.
> 
Look at b[a].

Ray

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] load object

2005-01-13 Thread Ray Brownrigg
> From: Douglas Bates <[EMAIL PROTECTED]> Fri Jan 14 08:35:33 2005
> 
> Weiwei Shi wrote:
> > Hi,
> > I happen to re-write my codes to save memory and my
> > approach is write my obj into file first and later I
> > load it.
> > 
> > However, it seems like:
> > load(filename) can load the object but the function
> > returns the name of the object instead of the
> > reference to it. For example, I have an object called
> > r0.prune, which is saved by
> > save(r0.prune, file='r0.prune')
> > 
> > and later, I want to load it by using:
> > load('r0.prune')
> > but I need to put the reference to the object r0.prune
> > into a var or a list. I tried:
> > t<-load('r0.prune'),
> > and class(t) gave me a char, which means t stores the
> > name of obj instead of the obj itself.
> > 
> > Sorry for the dumb question but please help...
> > 
> > Weiwei
> 
> I was going to suggest that you read the help page for load but when I 
> looked at it myself I found that it was not obvious what the effect of 
> calling load at the prompt is.  The help page is accurate but it is a 
> bit confusing if you don't know what the default environment is.
> 
> Anyway, when called from the prompt, load has the side effect of loading 
> the object into the global environment.  Try
> 
Well ?load does say:
Value:

 A character vector of the names of objects created, invisibly.

Note the plurals.  The point is that the file being loaded may contain
the definition of more than one R object.  You can say:
tt <- get(load('r0.prune')) # t is not a good example name to use
which will do what you want, but:
1) load() has a side-effect of also creating r0.prune (in your case)
2) if the file 'r0.prune' contains more than one object, only the first
is assigned to tt (but all of them are loaded into memory with their
original names).

Hope this helps,
Ray Brownrigg

__
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] counting numbers without replicates in a vector

2004-12-16 Thread Ray Brownrigg
> I am just wondering if there is an easy way to count
> in a numeric vector how many numbers don't have
> replicates. 
> For example, 
> a=c(1,1,2,2,3,4,5), how can I know there are three
> numbers (3, 4 and 5) without replicates?
> 
How about:
length(table(a)[table(a) == 1])

?
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Strange error from R CMD INSTALL

2004-12-08 Thread Ray Brownrigg
> Date: Wed, 08 Dec 2004 09:33:22 -0400
> From: Kjetil Brinchmann Halvorsen <[EMAIL PROTECTED]>
> 
> I am trying to install a local package and get this unexpected
> error:
> 
> -- Making package UMSA 
>   adding build stamp to DESCRIPTION
>   installing R files
>   installing data files
>   installing man source files
>   installing indices
> Error: couldn't find function "na.omit"
> Execution halted
> 
> na.omit of course is in package stats, and that is listed in the
> Depends field in the DESCRIPTION file.
> 
Check what is in your data directory.  Does something in there use
(implicitly) na.omit?

I have seen something like this (since 2.0.0) where a .R file in the
data/ directory used a function within the package, but that function
was not 'available' at the "installing indices" phase of the
check/build.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Manipulating contour plots

2004-12-01 Thread Ray Brownrigg
From: Caroline Houldcroft <[EMAIL PROTECTED]>
Date: 01 Dec 2004 15:16:55 +

> I would like to create a filled.contour plot with colour bar and overlay
> a map.
> 
> I need to create several of these representing variables with different
> spatial extents for comparison. When I try to fix the plot area size and
> axes limits as follows:
> 
> par(fig=c(0,1.0,0.3125,1.0)
> 
> filled.contour(x=seq((320/2)-180.5,(480/2)-180.5,len=nrow(oct94)),
> y=seq((190/2)-90.5,(240/2)-90.5,len=ncol(oct94)),
> z=as.matrix(oct94),
> color=heat.colors,
> levels=seq(0,30,2),
> plot.axes={axis(1,seq(-60,60,by=20));
> axis(2,seq(5,30,by=5))})
> 
> 
> my definitions were ignored producing a stretched image of the area with
> limits defined by the contour plot.  
> 
I cannot answer this part of the question because your example is not
reproducible (Error in nrow(oct94) : Object "oct94" not found) [and the
first line contains a syntax error.]
> 
> I have also had difficulties overlaying a map on the same axes. The
> following example, taken from the r-help archives, gave the error
> "couldn't find function "world""
> 
> library(maps)
> filled.contour(x=seq((320/2)-180.5,(480/2)-180.5,len=nrow(oct94)),
> y=seq((190/2)-90.5,(240/2)-90.5,len=ncol(oct94)),
> z=as.matrix(oct94),
> color=heat.colors,
> levels=seq(0,30,2),
> plot.axes={world(ylim=c(5,30),xlim=c(-30,50),add=T);
> axis(1);axis(2)})
> 
The world() function is not part of the maps package, it is in the
fields package.  If you want to use the maps package, try:
map("world", ylim=c(5,30), xlim=c(-30, 50), add=T)

Hope this helps,
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] library(fields) world shift function not working anymore

2004-11-29 Thread Ray Brownrigg
> I just upgraded to 2.01 on Mac OS 10.3.6.  I used to use the command 
> (on R 1.9.x):
> 
> world(ylim=c(-30,30), xlim = c(0,360), shift=TRUE, add=TRUE)
> 
> to draw a world outline over my image plots.  My data uses longitude 
> from (0, 360) so I need to use the shift function.  After I upgraded, I 
> get the following error:
> 
>  > world(ylim=c(-30,30), xlim = c(0,360), shift=TRUE, add=TRUE)
> Error in world(ylim = c(-30, 30), xlim = c(0, 360), shift = TRUE, add = 
> TRUE) :
>   NAs are not allowed in subscripted assignments
> 
> Does anyone know a workaround for this?
> 
Well, a workaround would be:
library(maps)
map("world2", ylim=c(-30,30), xlim = c(0,360), add = TRUE)

(at least until fields is updated).

CHANGES IN R VERSION 2.0.0

o   Subassignments involving NAs and with a replacement value of
length > 1 are now disallowed.  (They were handled
inconsistently in R < 2.0.0, see PR#7210.)  For data frames
they are disallowed altogether, even for logical matrix indices
(the only case which used to work).

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


RE: [R] IFELSE across large array?

2004-11-23 Thread Ray Brownrigg
> From: "Liaw, Andy" <[EMAIL PROTECTED]>
> Date: Tue, 23 Nov 2004 12:28:48 -0500
> 
> I'll give it half a crack:
> 
> Steps a through c can be done via nested ifelse(), treating A and M as
> vectors (as they really are).  Step d is the hard one.  I'd write a simple
> Fortran code and use .Fortran() for that.
> 
> I don't see how any of the *apply() functions can help here, as your
> operations are element-wise, not dimension-wise.
> 
> Andy
> 
The original message mentions that the value 10 is actually "NODATA",
and if one recodes 10 as NA, steps a) to c) become trivial, namely:

A[A == 10] <- NA
M[M == 10] <- NA
return(A + M - 1)

Then if step d) is performed first (i.e. appropriate values in A are
replaced by the 'most common neighbour' [perhaps using
round(mean(.., na.rm=T))] this still works, but would have to be repeated
for each replication (the third dimension).

Ray Brownrigg

> > From: Sander Oom
> > 
> > Dear all,
> > 
> > As our previous email did not get any response, we try again with a 
> > reformulated question!
> > 
> > We are trying to do something which needs an efficient loop 
> > over a huge 
> > array, possibly functions such as apply and related (tapply,
> > lapply...?), but can't really understand syntax and examples in 
> > practice...i.e. cant' make it work.
> > 
> > to be more specific:
> > we are trying to apply a mask to a 3D array.
> > By this I mean that when "overlaying" [i.e. comparing element 
> > by element] 
> > the mask on to the array the mask should change array 
> > elements according to 
> > the values of both array and mask elements
> > 
> > the mask has two values: 1 and 10.
> > 
> > the array elements have 3 values: 0, 1,  or 10
> > 
> > sticking for the moment to the single 2d array case
> > 
> > for example:
> > [A= array]  100 10 1  10  0
> >   1   10   1 0   0 10
> > 
> > [ M=mask] 1  10  10 1   1  1
> >  101   1  1 10 10
> > 
> > I would like the array elements to:
> > 
> > a) IF A(ij) !=10 and  Mij = 1
> >   leave A(ij) unchanged
> > 
> > b)  IF   A(ij) != 10 and M(ij) =10
> > change A(ij) to M(ij) i.e mask value (10)
> > 
> > c)IF A(ij) = 10 and M(ij) = 10
> >leave (Aij) unchanged
> > 
> > d) IF A(ij) = 10 and M(ij) !=10
> > replace A(ij) with the majority value in the 
> > 8-neighborhood
> > 
> > (or whatever if it is an edge element) BUT ignoring 10s in this 
> > neighborhood (i.e. with either 1 or 0, whichever is in majority)
> > 
> > because the array is 3d I would like to repeat the thing with 
> > all the k 
> > elements (2d sub-arrays) of the array in question, using the 
> > same mask for 
> > al k elements
> > 
> > Would you be able to suggest a strategy to do this?
> > 
> > thanks very much
> > 
> > Alessandro and Sander.

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Re: 3d Map with bars

2004-11-21 Thread Ray Brownrigg
> From: [EMAIL PROTECTED]
> Date: Fri, 19 Nov 2004 16:53:12 -0500
> 
> Thanks for reply. I need to first draw the map of USA a perspective plot. 
> I guess thats where my problem was.
> 

Try something like this:

library(maps)
states <- map("state", plot=F)
x1 <- rep(0, 3)
x2 <- rep(0, 3)
maxz <- 1
z <- matrix(c(0, 0, 0, 0, 0, 0, maxz, maxz, maxz), 3, 3)
x1[-2] <- states$range[1:2]
x2[-2] <- states$range[3:4]
x1[2] <- x1[3] - 1e-6
x2[2] <- x2[3] - 1e-6
pmat <- persp(x1, x2, z, xlab = "", ylab = "", zlab="", axes = F,
  theta=0, phi=20, d=10)
lines(trans3d(states$x, states$y, 0, pmat))
latitude <- -90
longitude <- 37
myz <- 0.6
lines(trans3d(rep(latitude, 2), rep(longitude, 2), c(0, myz), pmat), col=2)

where trans3d is as defined in the examples for persp().

Hope this helps,
Ray Brownrigg


> Uwe Ligges <[EMAIL PROTECTED]>
> 11/19/2004 04:33 PM
> 
> To: [EMAIL PROTECTED]
> cc: [EMAIL PROTECTED]
> Subject:Re: 3d Map with bars
> 
> 
> [EMAIL PROTECTED] wrote:
> 
> > Apologies in advance for the question. I am trying to draw a map of the 
> US
> > as a surface plot so that I would be able to drop bars on the different
> > states (something like Uwe Ligges' scatterplot3d example 4). I am not 
> sure
> > where to start looking for such a beast. If anyone has any pointers,
> > ideas, I will be grateful.
> >
> > TIA,
> > Partha
> 
> How to "drop bars" with persp() has been described on R-help yesterday
> or today, please check the mailing list's archives.
>

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] adjusting the map of France to 1830

2004-11-21 Thread Ray Brownrigg
> Date: Fri, 19 Nov 2004 15:59:25 -0500
> From: Michael Friendly <[EMAIL PROTECTED]>
> 
> Here's what I tried.   I can plot a selection of regions, but I
> can't seem to remove an arbitrary list of region numbers, unless I've 
> done something wrong
> by selecting the regions I want to plot with departements[-exclude].

I think here the problemis not using exact=T in the call to map(), see
below.

> I also get an error
> when I try to use map.text to label a map with only the regions I'm  
> selecting.
> 
>  > departements <- map('france',namesonly=T, plot=FALSE)
>  > # returns a vector of names of regions
>  >
>  > exclude <- c(47,  #Alpes-Maritimes
> + 66,  # Haute-Savoie
> + 76,  # Savoie
> + 95,  # Territore-de-Belfort
> + 109, 110, 111, # Var: Iles d'Hyeres
> + 49, 53, 54, 55, # Moribhan: Isles
> + 62, 64,# Vendee: Isles
> + 72, 75 # Charente-Maritime: Isles
> + )
>  >
>  > depts <- departements[-exclude]
>  > gfrance <-map('france', regions=depts)
>  > labels <- (as.character(1:length(departements)))[-exclude]
>  > gfrance <-map.text('france', regions=depts, add=FALSE, labels=labels)
> Error in map.text("france", regions = depts, add = FALSE, labels = labels) :
> map object must have polygons (fill=TRUE)
> 
That error message is issued when regions= specifies less than the whole
database, in which case the alternate "list of 'x', 'y', and 'names'
obtained from a previous call to 'map'" is required as a first parameter.

So to do what you want (notwithstanding the next problem you raise), try
the following as a demonstration:

gfrance <- map('france', regions=depts, exact=T, fill=T, plot=F)
map.text(gfrance, regions=depts, labels=labels)  
map('france', regions=departements[exclude], fill=T, col=1, add=T)

> Another problem, potentially more difficult for mapping data on the map 
> of France is that
> the "departements" are actually just the polygons in the map, 
> arbitrarily numbered from
> east to west, and from north to south --- they don't correspond to the 
> 'official' administrative
> region numbers.  As well, the departement names don't always match 
> exactly (ignoring
> accents, e.g., Val-d'Oise vs. Val-Doise) so it would be another 
> challenge to plot my
> historical data on the map of France.
> 
Well, maps is a source package! [:-)].

You are most welcome to modify the source files
maps/src/france.{gon,line,name} to reorder the polygons (and correct
errors in the names).  If the relationship between those 3 files is not
obvious, contact me for further details.  Also, I am happy to fold your
changes back into the original maps package.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] adjusting the map of France to 1830

2004-11-18 Thread Ray Brownrigg
At 16:29 18/11/2004, Michael Friendly wrote:
> I'm doing some analyses of historical data from France in 1830 on 'moral 
> statistics' that I'd like to
> show on a map.  I've done most of my analyses in SAS, but a few things 
> would work better in R.
> To do this, I have to adjust the modern map,
> 
> library(maps)
> map('france')
> 
> to adjust for changes in departments (86 in 1830, to 97 now).  I've read 
> the documentation
> for the maps and maptools package, but there seems to be no functions to 
> allow this, and
> I can't find information on the exact structure of map datasets, but I 
> understand them to
> be delimited lists of polygon coordinates.
> 
> In SAS, all maps have (one or more) ID variables representing the 
> geographical region,
> and there is also a proc gremove that can remove internal boundaries 
> inside the polygons
> for regions with the same ID.  Is there some way I can do this in R?
>
Unfortunately not with the current implementation of several of the
'extra' databases in the mapdata package.  The map() function does have
the interior=FALSE option, which would normally do what you want, but
only when the data has been manipulated to allow it.  Currently this
option is only useful with the world and usa maps (and their
derivatives, such as world2 and state).

Currently every department is a complete polygon, and so every interior
line segment occurs twice in the data.  What has to happen to the data
is for it to be split up into non-overlapping line segments, and each
polygon reconstructed from a list of these line segments (with
direction being important).

If you are prepared to perform this somewhat tedious process, I am
happy to assist you with further details.

However even with the interior= option functioning, it would still not
be easy to produce the map you would require, since you would have to
build it up from many components (namely each of the 'combined'
departments, plus 'all the rest').

HTH,
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


RE: [R] Seg. faults in mapthin (in package maps!)

2004-07-18 Thread Ray Brownrigg
> From: "Kevin Bartz" <[EMAIL PROTECTED]>
> Date: Thu, 1 Jul 2004 10:18:36 -0700
> 
> Dr. Ripley, you are the master. That fix worked like a charm! All the way to
> 50, with no problems. Thanks again,
> 
Hear, hear!

> From: Prof Brian Ripley [mailto:[EMAIL PROTECTED] 
> Sent: Thursday, July 01, 2004 10:07 AM
> 
> There are a few other quirks, and I would check some of the other uses of 
> `long' in the C sources.
> 
OK, I've changed most of the 'long's to 'int's in the next version
(maps_2.0-22), soon to appear on CRAN.  It passes R CMD check on my 32-bit
systems, but please test it thoroughly on a 64-bit system.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


RE: [R] Seg. faults in mapthin (in package maps!)

2004-07-01 Thread Ray Brownrigg
> From: "Kevin Bartz" <[EMAIL PROTECTED]>
> Date: Thu, 1 Jul 2004 10:18:36 -0700
> 
> Dr. Ripley, you are the master. That fix worked like a charm! All the way to
> 50, with no problems. Thanks again,
> 
And all while I was asleep!   I'll submit a corrected maps package in a
couple of weeks, when I return from China (no guaranteed connection
where I'll be).

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Greater than 1 or less than 1?

2004-06-20 Thread Ray Brownrigg
> Date: Sun, 20 Jun 2004 18:32:46 -0400
> From: "Zhen Chen" <[EMAIL PROTECTED]>
> 
> I have problem evaluating the expression h = exp(x)/(exp(exp(x))-1) for 
> large negative x. This expression is actually the probability
> that y = 1 when y is a Poisson random variable truncated at 0, hence 
> must satisfy 0 <= h <= 1. However, when
> x < -18, I may get an h value that is  larger than 1 while the true 
> value should be a quantity that is smaller than but very close to 1.
> For example when x = -19, h = 1.031174. I tried to use 
> different form of h and none of
> them give me an h value that is less than or equal to 1. I also tried 
> to find patterns, but discovered none: for some x < -18, h is below 1; 
> for others
> h is greater than 1. Is there any trick that enables me to obtain 
> theoretically correct results from this expression?
> 
exp(x)/(exp(exp(x)) - 1) is less than 1 for me for x = -18 with R-1.9.0
on a Solaris 9 system, a NetBSD 2.0C system and a Windows XP system.

You could try exp(x)/expm1(exp(x)), but if the longer expression doesn't
work, then this may not either, though it does give slightly different
results on my systems.

E.g.:
> x = -19
> exp(x)/(exp(exp(x)) - 1) - 1
[1] -2.047911e-09
> exp(x)/(expm1(exp(x))) - 1
[1] -2.801398e-09

log(1 + x) and exp(x) - 1 are often treated specially in math libraries
for when x << 0.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] R on Tru64 OSF 5.1

2004-04-02 Thread Ray Brownrigg
> From: =?iso-8859-1?q?J=20Swinton?= <[EMAIL PROTECTED]>
> 
> Has anyone achieved a current successful build of R
> 1.8 or 1.9 or earlier for HPUX Tru 64 OSF 5.1?
> 
> There is no binary version of R for HPUX Tru 64 OSF
> 5.1. The R admin manual mentions that the native make
> fails on "Alpha/OSF (aka Tru 64)" and gnu make must be
> used instead. There are problems reported with
> building R versions around 1.4 from 2002 on R-devel,
> some of which were fixed with using gmake instead of
> make, some of which were fixed by changes to R
> sources, but it is unclear if current R compiles on
> the current OSF 5.1, which is what I need to know
> before trying to get our local IS admins to try again
> after their first three failures.
> 
> Jonathan
> 
It depends what you mean by "successful build".  I succeeded in getting
R-1.8.1 to "./configure" and "make" on a system:
Compaq Tru64 UNIX V5.1A (Rev. 1885); Wed Nov 20 00:37:24 NZDT 2002

However it does not pass  "make check", with the following message:
496306:/vol/R/src/osf1/R-1.8.1/bin/R.bin: /sbin/loader: Fatal Error: call to 
unresolved symbol from /vol/R/src/osf1/R-1.8.1/modules/lapack.so (pc=0x3ffbfde3308)

during reg-tests-1.

This is possibly a consequence of the following (Info) warnings issued
at make time:

f77 -fpe3   -g -c dlapack0.f -o dlapack0.lo
f90: Info: dlapack0.f, line 1721: This directive is not supported in this platform.
CDEC$  NOVECTOR
---^
f77 -fpe3   -g -c dlapack1.f -o dlapack1.lo
f77 -fpe3   -g -c dlapack2.f -o dlapack2.lo
f77 -fpe3   -g -c dlapack3.f -o dlapack3.lo
f90: Info: dlapack3.f, line 21481: This directive is not supported in this platform.
CDEC$  NOVECTOR
---^
f90: Info: dlapack3.f, line 21495: This directive is not supported in this platform.
CDEC$ NOVECTOR
--^
f90: Info: dlapack3.f, line 21522: This directive is not supported in this platform.
CDEC$  NOVECTOR
---^

The weird thing is that the problem occurs during the following code
fragment:

## predict.arima0 needed a matrix newxreg: Roger Koenker, 2001-09-27
u <- rnorm(120)
s <- 1:120
y <- 0.3*s + 5*filter(u, c(.95,-.1), "recursive", init=rnorm(2))
fit0 <- arima0(y,order=c(2,0,0), xreg=s)

which runs perfectly OK if pasted into an interactive session.

I did spend some time trying to isolate/reproduce the problem without
success.  Unfortunately, this is the only Alpha/OSF1 system we have and
it is not used much (particularly not for R), so the installation has
not been well tested for 'real use'.

Let me know if you want the necessary configuration changes.

Hope this helps,
Ray Brownrigg <[EMAIL PROTECTED]>   http://www.mcs.vuw.ac.nz/~ray

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] writing text on graphics' window

2004-03-21 Thread Ray Brownrigg
> From: Duncan Murdoch <[EMAIL PROTECTED]>
> Date: Sun, 21 Mar 2004 06:58:03 -0500
> 
> On Sun, 21 Mar 2004 11:07:01 +, you [EMAIL PROTECTED] wrote:
> 
> >Does anyone know of a method for writing text to the graphics window, 
> >where there is *no* plot? Basically, I have developed a 'significance 
> >test' and I would like the output on the graphics window to say 
> >something about the input parameters and the stats of the significance test.
> 
> You need to make sure a graphics device is active and establish a
> coordinate system there.  The easiest way to do that is to make a call
> to plot() with everything turned off:
> 
>  plot(0:100,0:100,type='n',axes=FALSE,xlab="",ylab="")
> 
> You may also want to reduce the "margins" if you want your output to
> take up the full frame, e.g.
> 
>  oldmargins <- par(mar=c(0,0,0,0))
>  plot(0:100,0:100,type='n',axes=FALSE,xlab="",ylab="")
> 
An easier way to activate a graphics device and establish a coordinate
system is to call plot.new().


Try:
> plot.new()
> text(0, 0, "ABC")
> par("usr")
[1] -0.04  1.04 -0.04  1.04

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] R equiv to proc gremove in maps package

2004-03-15 Thread Ray Brownrigg
> Is there an R equivalent to SAS's proc gremove?  You would use this procedure to 
> combine the units on an existing map, for example to build a map of Metropolitan 
> Statistical Areas (MSAs) from the [US] counties dataset where the internal boundries 
> surround the MSAs (which are groups of counties) rather than the individual 
> counties.  I can imagine the mechanism would be to find and erase the boundries 
> within each group.  Is there a way to do that in R?  I looked at the S reference on 
> building geographical databases, and I could probably fumble through the C code to 
> build the binary files of lines and polygons, but have no idea how to implement the 
> shell program to extract the datapoints.  Thanks.
> 
It's already built in to the maps package.

Try e.g.:
map("county", c("arizona,gila", "arizona,maricopa", "arizona,yavapai"), interior=F)

[I wouldn't have believed that it worked until I tried!]

Hope this helps,
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Map of British Colonial America 1775

2004-03-10 Thread Ray Brownrigg
>   I've reviewed the "map" documentation I've found so far with 
> S-Plus 6.2 and R including "www.r-project.org" -> search -> "R site 
> search" -> "world map" and Becker & Wilkes (1993) "Maps in S" cited in 
> the "map" documentation.  I've noticed that, for example, 
> 'library(maps);  map("world", "UK:Bermuda")' produces an 83-point 
> polygon sketch in S-Plus 6.2 but only a degenerate 3-point polygon (with 
> the third point = the first) in R 1.8.1. 
> 
You need the mapdata package.  The default "world" database in S-Plus is
the high-resolution one (with also "world.thin").  In R the high
resolution one is "worldHires" from mapdata (mainly so that the
standard maps package is a manageable size).

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] drawing filled countries according to data using map('world')?

2004-03-07 Thread Ray Brownrigg
> My question then is, how to compute a similar procedure using the 'world'
> database. Specifically, how can I access the country names in the 'world'
> database to accomplish the translation to the country names in my dataset?
> Is there any way to unpack the 'world' database to do the matching in an
> external program? And does anybody now of other (more recent) world maps
> that I could use?
> 
The short answer is that the file world.N in .../library/maps/mapdata/
contains the mapping from polygon numbers to names.

I think the only other thing you need to know is that the colours are
allocated sequentially in increasing numerical order of polygon number.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] ISO country codes available in R's {map}?

2004-03-04 Thread Ray Brownrigg
> Does any one know wether it is possible to get the Country codes from the
> Map package - I've not found anything in the docs and only seen a post re.
> how to get country names.  I'd like to programmatically marry some data
> indexed by ISO country code to the map package contry names.  If not then
> I'll need to hunt for a listing and try and marry that to the {map}
> names.
> 
No part of the maps package knows anything about ISO country codes.  I
guess if there was a database of ISO country codes somewhere, it
wouldn't be too difficult to provide some add-on functionality to merge
the internal polygon numbers with country code numbers.  The hardest
part would be what to do with mismatches.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] Distance and Aggregate Data - Again...

2004-02-26 Thread Ray Brownrigg
[EMAIL PROTECTED] wrote:

> I appreciate the help I've been given so far.  The issue I face is
> that the data I'm working with has 53000 rows, so in calculating
> distance, finding all recids that fall within 2km and summing the
> population, etc. - a) takes too long and b) have no sense of progress.
>
Well, there are ways to speed this up.

> Below is a loop that reads each recid one at a time, calculates the
> distance and identifies the recids that fall within 2 km.  It iterates
> through all records successfully.
>
But you don't need to subset d[, 2:3] every time, e.g.  Also, my
experience is that writing a single line appended to a file every time
around the loop is very inefficient.

> Where I'm stuck is how to get the sum of population and dwellings and
> the mean age for the records that are selected.  Also, the desired
> output should have the following fields:  recid, sum(pop), sum(dwell),
> mean(age).  I don't know how to write only those fields out to the
> file.
>
Well, that part is easy, essentially, you have to save the indices, then
extract the relevant rows of pop, dwell and age. So I would modify your
code as follows:

library(fields)
d <- as.matrix( read.csv("filein.csv") )

lonlat2 <- d[,2:3]
results <- matrix(0, nrow=nrow(d), ncol=4)
for(i in 1:nrow(d)){
lonlat1 <- d[i, 2:3]
whichdist <- which(rdist.earth(t(as.matrix(lonlat1)),
as.matrix(lonlat2), miles=F) < 2)
distval <- d[, 1][whichdist]
sumpop <- sum(data[, "pop"][whichdist])
sumdwell <- sum(data[, "dwell"][whichdist])
meanage <- mean(data[, "age"][whichdist])
results[i, ] <- c(d[i, "recid"], sumpop, sumdwell, meanage)
}
write(t(results), file="C:\\outfile.out", ncol=ncol(results))

Then there is a trick you can use to speed up the process.  Essentially
you reduce the 'inner loop' which is inside rdist.earth().  This is
achieved by initially computing a single distance vector with reference
to a fixed point [(0, 0) seems reasonable since it is in the Atlantic
Ocean].  Then you select a subset of your points that are the same
distance is your circle centre from the fixed point plus or minus the
radius you want (and a bit of tolerance).  This will generate an
annulus of points rather than a circle, but in particular all the
points within the circle of interest will also be within this annulus.
Then you apply rdist.earth() to find the distance of each of these
points from your circle centre.

This is what I have used to process a 52000 length matrix in about 40
minutes (and ~50MB memory):

DistAgg <-
function(data, dist=2) {
  results <- matrix(0, nrow=nrow(data), ncol=4)
  colnames(results) <- c("recid", "sumpop", "sumdwell", "meanage")
  lonlat2 <- data[, 2:3]
  basedist <- rdist.earth(t(as.matrix(c(0, 0))), as.matrix(lonlat2))
  for(i in 1:nrow(data)){
lonlat1 <- data[i, 2:3]
approxval <- which(abs(basedist - basedist[i]) < dist*1.001)
if (length(approxval) > 1) {
  whichval <- approxval[which(rdist.earth(t(as.matrix(lonlat1)),
as.matrix(lonlat2[approxval, ]), miles=F) < dist)]
} else {
  whichval <- approxval
}
sumpop <- sum(data[, "pop"][whichval])
sumdwell <- sum(data[, "dwell"][whichval])
meanage <- mean(data[, "age"][whichval])
results[i, ] <- c(data[i, "recid"], sumpop, sumdwell, meanage)
  }
  write(t(results), file="C:\\outfile.out", ncol=ncol(results))
}

The data I used was based on yours, but randomised latitude and
longitude (in restricted ranges).

Hope this helps,
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] parse error in GLMM function

2004-02-18 Thread Ray Brownrigg
"CENDOYA, Gabriela" <[EMAIL PROTECTED]> writes:
> I'm trying to use the function GLMM from lme4 package, (R-1.8.1, Windows
> 98),and I get the following error:
> > pd5 = GLMM(nplant~sitio+
> + fert+
> + remo+
> + sitio:fert+
> + remo:sitio+
> + remo:fert+
> + remo:fert:sitio
> + data=datos,
> + family=binomial,
> + random=~repe:sitio)
> Error in parse(file, n, text, prompt) : parse error

Isn't the problem that you really do have a syntax error?

Note there is no comma before the data=.

Ray

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] large fonts on plots

2004-02-11 Thread Ray Brownrigg
On 11 Feb 2004 at 13:30, Piet van Remortel wrote:

> Hi all,
> 
> I need to enlarge te fonts used oo R-plots (plots, histograms, ...) in
> labels and titles etc.

Nobody seems to have mentioned another way of achieving this, which is
to scale the original plot to suit the intended final size in the
document.  E.g. if you use:
> postscript(file="plot1.ps", width=4, height=4)
then the 4"x4" figure will 'automatically' have relatively larger
text and plot symbols.  This works well with x11() which is
subsequently printed, and pdf() also.

HTH,
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


Re: [R] filled maps

2004-02-03 Thread Ray Brownrigg
> Date: Tue, 3 Feb 2004 11:34:06 +0100 (CET)
> From: Roger Bivand <[EMAIL PROTECTED]>
> 
> On Tue, 3 Feb 2004, Janus Larsen wrote:
> 
> > Hi R-Help,
> > 
> > I would like to make filled contour maps of ocean data overlaid by
> > costlines from the map package.
> > I can draw the filled contours and the coastlines om the same plot, but
> > the filled contour also covers part of the land. To get rid of that I
> > tried to draw a filled coastline map on top of the filled contour, but
> > the filled map only draws the closed contours - so most of the land is
> > missing.
> > Example:
> > map("worldHires",xlim=c(0,15),ylim=c(50,60)) #Draw relevant region
> > (North Sea and Denmark waters)
> > map("worldHires",xlim=c(0,15),ylim=c(50,60),fill=TRUE) # This only draws
> > Denmark and Holland (Sweden, uk, Germany etc. disappears because they
> > are not closed polygons).
> 
> Looks as though you can use the regions= argument:
> 
> > res <- map("worldHires",xlim=c(0,15),ylim=c(50,60), plot=FALSE, 
> + namesonly=TRUE)
> > map("worldHires",xlim=c(0,15),ylim=c(50,60))
> > map("worldHires",regions=res, fill=TRUE, add=TRUE)
> 
Thanks Roger, I hadn't got around to figuring out that workaround.
Certainly the mapgetg code explicitly excludes incomplete polygons if 
fill=TRUE, but since the polygon() function is happy to clip what it is
filling, your workaround takes advantage of that.

I'll have a look at getting mapgetg to do this itself - it would have
to always provide complete polygons if fill=TRUE instead of just
'visible' polylines.

> fixes it for me. Curiously, the longitudes add 2.5 degrees on each side.
> 
Actually this is to preserve aspect ratio.  An earlier version of the
maps package adjusted the shape of the figure to conform to both aspect
ratio and imposed limits.  I didn't deliberately change this, but I
think I know why it has happened, so I'll look at this also.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html


RE: [R] beginner programming question

2003-12-17 Thread Ray Brownrigg
> From: "Gabor Grothendieck" <[EMAIL PROTECTED]>
> Date: Wed, 17 Dec 2003 15:02:49 -0500 (EST)
> 
> Define function f to take a vector as input representing
> a single input row.   f should (1) transform this to a vector 
> representing the required row of output or else (2) produce 
> NULL if no row is to be output for that input row.
> 
> Then use this code where z is your input matrix:
> 
> t( matrix( unlist( apply( z, 1, f ) ), 2) )
> 
But as has been pointed out recently, apply really is still just a for
loop.

> > From: Adrian Dusa <[EMAIL PROTECTED]>
> > Date: Wed, 17 Dec 2003 21:28:05 +0200 
> > 
> > I have a (rather theoretical) programming problem for which I have found
> > a solution, but I feel it is a rather poor one. I wonder if there's some
> > other (more clever) solution, using (maybe?) vectorization or
> > subscripting.

Here is a subscripting solution, where (for consistency with above) z is
your data [from read.table(filename, header=T)]:

> z
  rel1 rel2 rel3 age0 age1 age2 age3 sex0 sex1 sex2 sex3
113   NA   25   232   NA121   NA
2413   35   67   34   102212
3144   39   40   59   601221
44   NA   NA   45   70   NA   NA22   NA   NA
> res <- matrix(NA, nrow=length(z[, 1]), ncol=2,
 dimnames=list(rownames=rownames(z), colnames=c("ageh", "agew")))
> w <- w0 <- w1 <- w2 <- which(z[, c("rel1", "rel2", "rel3")] == 1, T)
# find spouse entries
> w0[, 2] <- z[, "sex0"][w[, 1]]# indices for respondent's age
> w1[, 2] <- 3 - w0[, 2]# indices for spouse's age
> w2[, 2] <- 4 + w[, 2] # indices of spouse's age
> res[w0] <- z[, "age0"][w[, 1]]# set respondent's age
> res[w1] <- z[w2]  # set spouse's age
> res
colnames
rownames ageh agew
   1   25   23
   2   34   35
   3   39   40
   4   NA   NA
>
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] Recoding problem

2003-12-08 Thread Ray Brownrigg
> I have the following variables, all of which are logicals
> 
> fmar15   fcoc15   fher15fcrk15fidu15
> 
> what I would like is a variable drug15 which equals 
> idu if fidu15 is T; crk if fidu15 is F but fcrk is T, her if fher15 is
> T but fcrk15 and fidu15 are F and so on
> 
> What's the best way to do this?  
> 
I don't know about the best way, but if I understand your question, the
following achieves what you want:

> tab <- cbind(fidu15, fcrk15, fher15, fcoc15, fmar15)
> substring(colnames(tab), 2, 4)[apply(tab, 1, match, x = T)]

HTH
Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] ?for

2003-11-20 Thread Ray Brownrigg
> I thought I was getting to grips with it, until last reply from Prof Brian
> Ripley's. At least in my Windows and Linux boxes

Actually it was my reply that suggested ?paste("help")

> ? paste("help")
> opens the help for help (!?)

Not on my systems! [Solaris, NetBSD, Windows 98].  I believe it should
open the help for paste().

The point is that it requires a syntactically correct statement after
the ? so that the command parses correctly (with ? as a unary
operator), but the result is to provide the help for the first
syntactic element of the parsed expression, *without* evaluating the expression.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] file not found?

2003-11-19 Thread Ray Brownrigg
Thomas W Blackwell <[EMAIL PROTECTED]> wrote:

> I recall a message on the [r-pkgs] list from Ray Brownrigg
> on November 1 this year that may be relevant.  At that time
> he was announcing the availability of  maps_2.0-8.zip  and
> some other packages.
> 
> Therefore,  2.0-10 must be *very* recent and might still be
> going through the CRAN propagation process.  (I don't know
> anything about whther it is or not.)  I would wait a couple
> of days and try again, if I were you.
> 
Well, actually, 2.0-10 in source form is on CRAN since last week, but
for some reason only in the 1.7/ directory in the Windows area.  I just
this morning noticed this and sent a message to the Windows CRAN
maintainer.

In the meantime, 2.0-11 has been uploaded, probably to appear in a few
days (no noticeable changes).

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] Problem with parser and if/else

2003-11-14 Thread Ray Brownrigg
Prof Brian Ripley <[EMAIL PROTECTED]>

> On Fri, 14 Nov 2003, Paul Lemmens wrote:
> 
> > I'm trying to grasp this: if you're saying (or are you saying) that the 
> > only way to have if() know that an else will be present is to put it on the 
> > same line as the closing curly brace of the if() (compound) statement. But 
> > if I look at some code from, e.g., aov and lm, I see plenty violations of 
> > that rule.
> 
> The actual rule is given in my reference (the one that Ben Bolker did not 
> bother to look up) earlier in this thread.
> 
> You need to ensure that the code is syntactically incomplete when `else'
> is encountered.  That will always be true inside a braced expression such 
> as the bodies of the functions you quote.  But at top-level, you do need
> to write
> 
> if(condition) something else something_else
> 
> or
> 
> if(condition) {
> } else {
> }
> 
> since
> 
> if(condition) {
> } 
> else {
> }
> 
> fails the test.
> 
BUT, just to make it even clearer, what does succeed interactively is
(what is effectively the same as a function body):

{
if(condition) {
"TRUE"
}
else {
"FALSE"
}
}
because the parser has to find the final closing brace before the
syntactic item is complete.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] stop further sourcing of an R file

2003-11-13 Thread Ray Brownrigg
> Thankyou, but i *donot* want the R session to quit.
> Someone suggested writing functions with return kept at
> the locations.
> I wanted to ask if there was another way.
> 
Just have a single
}
as your last line, and put:
if (FALSE) {

just after where you want to stop.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] Test for new page

2003-11-11 Thread Ray Brownrigg
Marc Schwartz <[EMAIL PROTECTED]> wrote:
> On Tue, 2003-11-11 at 16:49, Paul Murrell wrote:
> > Hi
> > 
> > 
> > Aleksey Naumov wrote:
> > > Dear R experts,
> > > 
> > > I am writing a multi-page PDF file and would like to put a header on each 
> > > page. Is there a way to test a graphic device to see if a new page is started 
> > > (so that I know when to write the header)?
> > 
> > 
> > Sorry.  Not that I can think of.
> > 
> > Paul
> > 
> > 
> > > I could simply count the plots made (each page has the same number of plots), 
> > > but wanted to see if a more general solution is available.
> 
> 
> I was trying to think of a way but could not either.
> 
It's not clear what you want as the "header", but if it is the same on
each page, then using mtext("Header", outer=TRUE) (after a suitable
par(oma=c(0, 0, 1, 0))), you can just write the header after every plot.

Ray

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] map does not display maps, MacOSX

2003-11-05 Thread Ray Brownrigg
Philippe Glaziou <[EMAIL PROTECTED]> wrote:
> 
> I installed the maps and mapdata libraries on my R-1.8.0 on
> MacOSX 10.2.8 (jaguar on a powerbook G4), and failed to make the
> map function work properly:
> 
> 
> R : Copyright 2003, The R Development Core Team
> Version 1.8.0  (2003-10-08)
> [...]
> 
> > library(maps)
> > map()
> Error in file(file, "r") : unable to open connection
> In addition: Warning message: 
> cannot open file `/Users/glaziou/Library/R/maps/mapdata//world.N' 
> > map('usa')
> Error in file(file, "r") : unable to open connection
> In addition: Warning message: 
> cannot open file `/Users/glaziou/Library/R/maps/mapdata//usa.N' 
> > system('ls -l /Users/glaziou/Library/R/maps/mapdata')
> total 1796
> -rw-r--r--   1 root staff  143902 Oct 14 11:30 county.G
> -rw-r--r--   1 root staff  690260 Oct 14 11:30 county.L
> -rw-r--r--   1 root staff 618 Oct 14 11:30 nz.G
> -rw-r--r--   1 root staff   13040 Oct 14 11:30 nz.L
> -rw-r--r--   1 root staff2642 Oct 14 11:30 state.G
> -rw-r--r--   1 root staff   96892 Oct 14 11:30 state.L
> -rw-r--r--   1 root staff 282 Oct 14 11:30 usa.G
> -rw-r--r--   1 root staff   58232 Oct 14 11:30 usa.L
> -rw-r--r--   1 root staff   74434 Oct 14 11:30 world.G
> -rw-r--r--   1 root staff  295152 Oct 14 11:30 world.L
> -rw-r--r--   1 root staff   74434 Oct 14 11:30 world2.G
> -rw-r--r--   1 root staff  295152 Oct 14 11:30 world2.L
> -rw-r--r--   1 root staff   54832 Oct 14 11:30 world2.N
> 
Note the actual file requested (usa.N) doesn't actually exist!

> :
> 
> Any hint appreciated,
> 
I had this sort of problem on a Windows system, and had to add an extra
command into Makefile.win.  [Is there recognised such a thing as
Makefile.mac? I seem to recall there used to be, but it is not mentioned
in the latest R-exts.]

What you need to do is add the line:

$(CP) ${*}.n ../inst/mapdata/${*}.N     # need this here for Mac

as the second line of the .gon.g: target in maps/src/Makefile.
[note the first character is a tab]

I'll fix this (and the //) in the next version.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] USA map

2003-11-05 Thread Ray Brownrigg
> Ivar Herfindal <[EMAIL PROTECTED]> wrote:
> 
> On Wed,  5 Nov 2003 11:28:42 +0100 (MET), Mathieu Ros 
> <[EMAIL PROTECTED]> wrote:
> 
> >>>>>> "k" == kjetil  <[EMAIL PROTECTED]> disait:
> >
> > 
> > k> I also tried
> >
> > k>> map("worldHires","sweden")
> > k>> map("worldHires","denmark") # which comes out very small since it
> > k>  # includes the Faroe k>   
> > # islands properly faraway
> >
> > and, just to know, how would you do to plot *only* continental
> > denmark? The same applies for france, UK, ...
> >
> Hello
> 
> One simple way of doing it is to specify the xlim and ylim of your map, 
> e.g. library(maps)
> map('world', 'Norway', xlim=c(5, 33), ylim=c(55, 75))
> 
But the 'best' way is RTFM!
map("world", "Norway", exact=T)

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] USA map

2003-11-04 Thread Ray Brownrigg
> From [EMAIL PROTECTED] Wed Nov  5 11:31:52 2003
> 
> but I cannot find any function to find the avalable region names in
> a map database
> 
They are in the .N file for the appropriate map.

Try:
library(maps)
data(worldMapEnv)
file.show(paste(Sys.getenv(worldMapEnv), "world.N", sep="/"))

These ASCII files are in the form:
region[:subregion] polygonNumber

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] USA map

2003-11-04 Thread Ray Brownrigg
Prof Brian Ripley <[EMAIL PROTECTED]> wrote:
> 
> On Tue, 4 Nov 2003 [EMAIL PROTECTED] wrote:
> 
> > On 3 Nov 2003 at 15:04, [EMAIL PROTECTED] wrote:
> > 
> > > 
> > > I believe in R you need
> > > 
> > > library(maps)
> > > 
> > > before using usa().  Of course, the library has to be installed first.
> > > 
> > 
> > Yes, I new that. When I did help.search("usa") I HAD maps
> > attached. It dosen't have any function usa(). R is more international 
> > minded:
> > 
> > > library(maps)
> > > usa()
> > Error: couldn't find function "usa"
> > > map('usa')
> > > 
> > map('china')
> > map('bolivia')
> > 
> > and so on also works!
> 
> You need library(mapdata) for those.
> 
Except the Bolivia map must be a local one (from [EMAIL PROTECTED]);
it is not in library(mapdata).

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] USA map

2003-11-03 Thread Ray Brownrigg
> > I don't know the usa() function in S but usa() in R gives a map of the US
> 
> Where did you find that function? In rw1080 on windows XP I get:
> 
Some earlier versions of the maps package (which was then available for
"Unix" only), did have individual functions for each country in the
world database.

In the current version, you need to type:
library(maps)
map("usa") 
#or:
map("world2", "USA")

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


[R] [R-pkgs] R maps, mapdata and mapproj for Windows

2003-11-02 Thread Ray Brownrigg
Nearly a month ago I announced the availability of a Windows port of
the maps for S package but for various reasons the Windows binaries
did not appear on CRAN until yesterday.  Now I am pleased to announce
that the three packages maps_2.0-8 (the base maps package with
low-resolution databases), mapdata_2.0-6 (high-resolution and other
databases) and mapproj_1.1 (projection add-on for maps, maintained by
Tom Minka) are all now available on CRAN both in source and Windows
binary forms.

I apologise for the delays.
Ray Brownrigg

___
R-packages mailing list
[EMAIL PROTECTED]
https://www.stat.math.ethz.ch/mailman/listinfo/r-packages

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


[R] [R-pkgs] R maps now for Windows

2003-10-05 Thread Ray Brownrigg
The original S maps package, previously only available in R on Unix
systems is now available in R for Windows.

maps_2.0-2 and mapdata_2.0-2 (see below) in both source and Windows
binary form (for rw1071) have been uploaded to CRAN, hopefully soon to
appear.

The key aspect of this package which differentiates it from other
GIS related R packages under development is that it is standalone
within R, i.e. it does not rely on the presence of any other
software.

Features are:
1) Built-in databases for usa, state (states of the USA), county
(counties of the USA), world, world2 (Pacific-centric world) and nz.
2) Databases supplied in source form with accompanying software to
generate the binary forms, allowing users to edit the data or add new
databases (basic format is latitude/longitude pairs).
3) A database of world cities (generally >~100,000 population).
4) A companion package called mapdata provides high-resolution versions
of world, world2 and nz, as well as china (provinces) and rivers.
5) A companion package called mapproj provides a variety of map 
   projections (from Doug McIlroy).
6) While being backwards-compatible with the original 'S maps library',
it also has many new features, including:
   - displaying region names on a map,
   - finding the map region containing a given point (point-in-polygon),
   - computing areas of map regions, and
   - interpolating data aggregated over map regions (for statistical maps).

The maps and mapdata packages are maintained by Ray Brownrigg with
generous assistance from Tom Minka.

The mapproj package is available at 
   http://www.stat.cmu.edu/~minka/software/maps/

Ray Brownrigg
Tom Minka
October 2003

___
R-packages mailing list
[EMAIL PROTECTED]
https://www.stat.math.ethz.ch/mailman/listinfo/r-packages

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] trouble with maps

2003-07-23 Thread Ray Brownrigg
[EMAIL PROTECTED] wrote:
> On Wed, 23 Jul 2003, Deborah Swayne wrote:
> 
> > Has anyone else seen this behavior from the "maps" package?
> > 
> > map('state', fill=TRUE)
> > 
> > results in a lively mix of overlapping polygons inside a
> > map of the US, but they have no obvious relationship to
> > state boundaries.  (See attached jpeg.)
> 
Ah, this was a 'known problem' :-) first noticed in January this year,
by Ott Toomet <[EMAIL PROTECTED]>.

I suspect it is a feature of the port to R which never worked
properly.  The solution is to reverse the sign of every line number
making up the polygons in the .gon file (or alternately, reverse the
order of the line numbers there).

A 'fixed' maps package (Unix only) is available at:
ftp://ftp.mcs.vuw.ac.nz/pub/statistics/map/maps_1.1-2.tar.gz

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] maps library for R?

2003-06-13 Thread Ray Brownrigg
> Does it already exist a library in R to draw maps (something like a Generic 
> Mapping Toolbox, http://gmt.soest.hawaii.edu/ port).
> 
> I've seen in an old R-help that Ross Ihaka once tried to port the S-plus map 
> library (http://maths.newcastle.edu.au/~rking/R/help/99b/0832.html ).
> Anybody know if this package is available somewhere or if there is somebody 
> developing a Mapping package for R. I'll might then try to help.

There are two, but AFAIK neither of them will work on Windows.  My port
of the basic S maps package (no projections) has recently been modified
to provide a small (700KB) 'base' package with limited maps data (just
usa, state, counties, nz and world.thin), plus an add-on package called
mapdata (24MB) with the high-resolution maps. [This split has not yet
been fully tested, and requires maps to be *installed* before mapdata,
and mapdata must be installed in the same library directory as maps.]

These are available from
ftp://ftp.mcs.vuw.ac.nz/pub/statistics/map/map*_1.1-0*

The other solution in development is called Rmap from Barry Rowlingson
(http://www.maths.lancs.ac.uk/Software/Rmap/), but as I understand it,
this requires quite a lot of "third-party" stuff to be installed first.

Ray Brownrigg

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] Combining the components of a character vector

2003-04-04 Thread Ray Brownrigg
>Suppose I have a character vector.
> 
> x <- c("Bob", "loves", "Sally")
> 
> I want to combine it into a single string:  "Bob loves Sally" .
> paste(x) yields:
> paste(x)
> [1] "Bob"   "loves" "Sally"
> 
> Is there a more natural (no loop) way to do this in R?
> 
RTFM:
> paste(x, collapse=" ")
[1] "Bob loves Sally"

Ray

__
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help


Re: [R] inserting elements in a list

2003-02-19 Thread Ray Brownrigg
> From: Peter Wolf <[EMAIL PROTECTED]>
> 
> "insert.values" <- function(x,pos.insert,x.insert){
> :

Which when reduced to its bare minimum to solve the stated problem,
becomes:
  threes <- which(a==3)
  a <- c(a, rep(7, length(threes)))[order(c(seq(a), threes))]

This is by far the overall fastest solution so far proposed, and
certainly is reasonably easy to understand.  However an even faster
solution is based on rep():

 N <- which(threes <- a == 3)   # which ones
 a <- rep(a, times = 1 + threes)# duplicate them
 a[N + 1:length(N)] <- 7# replace the duplicates

Relative timing for vector lengths n=1000 and n=1 with 1000 repeats
(on a 1.7GHz P4 running R-1.6.1 under NetBSD) is as follows:

a <- sample(1:9, n, TRUE)   n=1000  n=1
lapply (Stephen Upton and Peter Dalgaard) 9.2s   117.8s
recursion (John Fox) [options(expressions=n)]59.9s   stack overflow
offset <- (Thomas Lumley) 1.4s   12.7s
offset <- (Peter Dalgaard)1.3s   10.9s
order() (Peter Wolf)  1.1s8.7s
rep() (me)        0.9s7.6s

Hope this helps,
Ray Brownrigg <[EMAIL PROTECTED]>   http://www.mcs.vuw.ac.nz/~ray

__
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help



Re: [R] Interpolation

2003-02-12 Thread Ray Brownrigg
> I have two vectors and connect the points by line segments:
> 
> x <- c(1990,1994,1995,1997)
> y <- c(30,40,80,20)
> plot(x,y,type="l")
> 
> I want to get the values of y's on these lines at missing x's,
> i.e., at 1991,1992,1993,and 1996. Is there a simple way to do this?
> 
Easy:
approx(x, y, c(1991:1993, 1996))$y
while
approx(x, y, 1990:1997)$y
will give you all the y values.

Hope this helps,
Ray Brownrigg <[EMAIL PROTECTED]>   http://www.mcs.vuw.ac.nz/~ray

__
[EMAIL PROTECTED] mailing list
http://www.stat.math.ethz.ch/mailman/listinfo/r-help