Re: [R] help with grouping data and calculating the means

2018-11-15 Thread Anthoni, Peter (IMK)
Hi Sasa,

Those latitude look equidistant with a separation of 0.05.
I guess you want to calculate the zonal mean along the latitude, right?

#estimate the lower and upper latitude for the cut
lat.dist=0.05 #equidistant spacing along latitude
lat.min=min(df$LAT,na.rm=T)-lat.dist/2
lat.max=max(df$LAT,na.rm=T)+lat.dist/2
cat.lat=cut(df$LAT,breaks=seq(lat.min,lat.max,by=lat.dist));cat.lat

#just show which indices are grouped
tapply(df$TK.QUADRANT,cat.lat, paste,collapse=",")

#calculate the mean of whatever column. The lat.mean will have NA for any 
latitude cell where the df column has no data
lat.mean=tapply(df$TK.QUADRANT,cat.lat, mean)

#if you need to remove any potential NAs
lat.mean[!is.na(lat.mean)]

cheers/beste Grüße
Peter

On 15. Nov 2018, at 17:48, sasa kosanic 
mailto:sasa.kosa...@gmail.com>> wrote:

`TK-QUADRANT` <- c(9161,9162,9163,9164,10152,10154,10161,10163)
LAT <- c(55.07496,55.07496,55.02495,55.02496
,54.97496,54.92495,54.97496,54.92496)
LON <-
c(8.37477,8.458109,8.37477,8.45811,8.291435,8.291437,8.374774,8.374774)
df <- data.frame(`TK-QUADRANT`=`TK-QUADRANT`,LAT=LAT,LON=LON)


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Help needed with aggregate or other solution

2017-10-26 Thread Anthoni, Peter (IMK)
Hi Thomas,

Would this work:
res1=aggregate(dta[,"fcst"],by=list(basistime=dta[,"basistime"]),FUN=max)
mm=match(paste(res1[,"basistime"],res1[,"x"]),paste(dta[,"basistime"],dta[,"fcst"]))
dta[mm,]
> dta[mm,]
  date   basistime fcst usgs
20 2012-01-30 12:00:00 2012-01-25 15:02:00 47.9 44.5
40 2012-01-31 12:00:00 2012-01-26 15:11:00 50.4 44.2
50 2012-01-29 12:00:00 2012-01-27 01:41:00 43.8 43.8

cheers
Peter



> On 26. Oct 2017, at 17:07, Jeff Newmiller  wrote:
> 
> On Thu, 26 Oct 2017, Thomas Adams wrote:
> 
>> Hi Jeff,
>> Thank you for the suggestions -- I appreciate your help. Unfortunately, the
>> result2 has two problems...
>> (1) there are now 3 date columns (it looks like 2 cols are merged into 1
>> col)
> 
> No, there are two date columns. Result2 includes the grouping value as a row 
> name (pulled from the names of the dta2list items by rbind).
> 
>> (2) the output rows should not have any of the basistime dates repeated
>> (maybe I misstated the problem); I need the max fcst value by basistime, but
>> also list the date value for that row; for example:
> 
> Then try out my code replacing
> 
> dta2list <- split( dta2, dta2$date )
> 
> with
> 
> dta2list <- split( dta2, dta2$basistime )
> 
> and
> 
> %>% group_by( date )
> 
> with
> 
> %>% group_by( basistime )
> 
> Please study how the code works and ask questions based on gaps in that 
> knowledge rather than how the results are not what you expected. This mailing 
> list is not a do-your-work-for-you coding service. Some help pages you should 
> look at include:
> 
> ?rownames
> ?split
> ?lapply
> ?do.call
> ?rbind
> ?group_by
> ?do
> 
> You might also find [1] helpful in general, and [2] helpful for understanding 
> dplyr.
> 
> [1] H. Wickham, The Split-Apply-Combine Strategy for Data Analysis, Journal 
> of Statistical Software, vol. 40, no. 1, Apr. 2011.
> 
> [2] H. Wickham and G. Grolemund, R for Data Science. OReilly UK Ltd, 2017. 
> URL: https://r4ds.had.co.nz.
> 
>>  basistime fcst
>> 1   2012-01-25 15:02:00 47.9
>> 2   2012-01-26 15:11:00 50.4
>> 3   2012-01-27 01:41:00 46.0
>> 4   2012-01-27 10:15:00 47.3
>> 5   2012-01-27 15:15:00 47.3
>> 6   2012-01-28 14:22:00 46.2
>> 7   2012-01-29 13:33:00 45.8
>> 8   2012-01-30 14:11:00 44.8
>> 9   2012-01-31 14:24:00 43.9
>> 10  2012-02-01 14:55:00 41.1
>> 11  2012-02-02 14:56:00 38.1
>> 12  2012-02-03 14:40:00 36.2
>> 13  2012-02-04 15:01:00 34.7
>> 14  2012-02-05 15:04:00 33.1
>> 15  2012-02-06 14:37:00 32.2
>> This is very close to what I need. The basistime dates are all unique, with
>> the max fcst value for the available basistime dates; but I additionally
>> need the corresponding 'date' value.
>> Best,
>> Tom
>> On Thu, Oct 26, 2017 at 1:28 AM, Jeff Newmiller 
>> wrote:
>>  Thanks for the dput...
>> 
>>   reproducible example of split-apply-combine ###
>> 
>>  dta <- structure(list(date = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
>>  7L,
>>  8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L,
>>  5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L,
>>  19L, 20L, 21L, 22L, 23L, 24L, 7L, 8L, 9L, 10L, 11L, 12L, 13L,
>>  14L, 15L, 16L), .Label = c("2012-01-25 18:00:00", "2012-01-26
>>  00:00:00",
>>  "2012-01-26 06:00:00", "2012-01-26 12:00:00", "2012-01-26
>>  18:00:00",
>>  "2012-01-27 00:00:00", "2012-01-27 06:00:00", "2012-01-27
>>  12:00:00",
>>  "2012-01-27 18:00:00", "2012-01-28 00:00:00", "2012-01-28
>>  06:00:00",
>>  "2012-01-28 12:00:00", "2012-01-28 18:00:00", "2012-01-29
>>  00:00:00",
>>  "2012-01-29 06:00:00", "2012-01-29 12:00:00", "2012-01-29
>>  18:00:00",
>>  "2012-01-30 00:00:00", "2012-01-30 06:00:00", "2012-01-30
>>  12:00:00",
>>  "2012-01-30 18:00:00", "2012-01-31 00:00:00", "2012-01-31
>>  06:00:00",
>>  "2012-01-31 12:00:00", "2012-01-31 13:00:00", "2012-01-31
>>  18:00:00",
>>  "2012-02-01 00:00:00", "2012-02-01 06:00:00", "2012-02-01
>>  12:00:00",
>>  "2012-02-01 18:00:00", "2012-02-02 00:00:00", "2012-02-02
>>  06:00:00",
>>  "2012-02-02 12:00:00", "2012-02-02 18:00:00", "2012-02-03
>>  00:00:00",
>>  "2012-02-03 06:00:00", "2012-02-03 12:00:00", "2012-02-03
>>  18:00:00",
>>  "2012-02-04 00:00:00", "2012-02-04 06:00:00", "2012-02-04
>>  12:00:00",
>>  "2012-02-04 18:00:00", "2012-02-05 00:00:00", "2012-02-05
>>  06:00:00",
>>  "2012-02-05 12:00:00", "2012-02-05 18:00:00", "2012-02-06
>>  00:00:00",
>>  "2012-02-06 06:00:00", "2012-02-06 12:00:00", "2012-02-06
>>  18:00:00",
>>  "2012-02-07 00:00:00", "2012-02-07 06:00:00", "2012-02-07
>>  12:00:00",
>>  "2012-02-07 18:00:00", "2012-02-08 00:00:00", "2012-02-08
>>  06:00:00",
>>  "2012-02-08 12:00:00", "2012-02-08 18:00:00", "2012-02-09
>>  00:00:00",
>>  "2012-02-09 06:00:00", "2012-02-09 12:00:00", "2012-02-09
>>  18:00:00",
>>  "2012-02-10 00:00:00", "2012-0

Re: [R] submitting R scripts with command_line_arguments to PBS HPC clusters

2017-07-11 Thread Anthoni, Peter (IMK)
Hi,

The problem is most likely, you need to call a R CMD BATCH with your arguments 
and the R-script inside of a shell script that you submit to your qsub.
Unfortunately we don't use qsub anymore so can't test it, but it should be as 
follows:

R-script eg. test.R:
> ##First read in the arguments listed at the command line
> args=(commandArgs(TRUE))
> 
> ##args is now a list of character vectors
> ## First check to see if arguments are passed.
> if(length(args)==0){
>   stop("no args specified")
> }
> ## Then cycle through each element of the list and evaluate the expressions.
> for(i in 1:length(args)){
>   print(args[[i]])
>   eval(parse(text=args[[i]]))
> }
> print(TUMOR)
> print(GERMLINE)
> print(CHR)


qsub shell script test.sh:  
> #!/bin/bash
> 
> #Note: the single quote '...' around the --args ... "..." "..." is important!
> R CMD BATCH --no-save --no-restore '--args TUMOR="tumor.bam" 
> GERMLINE="germline.bam" CHR="chr22"' test.R test.Rout

then you submit with a qsub with all the options you specified the test.sh
qsub  test.sh

cheers
Peter



> On 12. Jul 2017, at 03:01, Jeff Newmiller  wrote:
> 
> This sounds like an operating system specific question, in that "submit the R 
> script to a PBS HPC scheduler" would be the kind of action that would run R 
> with very different environment variables and possibly different access 
> credentials than your usual interactive terminal.  A thorough reading of the 
> "Installation and Administration Guide" and some study of your HPC 
> documentation are in order. 
> -- 
> Sent from my phone. Please excuse my brevity.
> 
> On July 11, 2017 5:25:20 PM PDT, Bogdan Tanasa  wrote:
>> Dear all,
>> 
>> please could you advise me on the following : I've written a R script
>> that
>> reads 3 arguments from the command line, i.e. :
>> 
>> " args <- commandArgs(TRUE)
>> TUMOR <- args[1]
>> GERMLINE <- args[2]
>> CHR <- args[3] ".
>> 
>> when I submit the R script to a PBS HPC scheduler, I do the following
>> (below), but ... I am getting an error message.
>> (I am not posting the error message, because the R script I wrote works
>> fine when it is run from a regular terminal ..)
>> 
>> Please may I ask, how do you usually submit the R scripts with command
>> line
>> arguments to PBS HPC schedulers ?
>> 
>> qsub -d $PWD -l nodes=1:ppn=4 -l vmem=10gb -m bea -M tan...@gmail.com \
>> -v TUMOR="tumor.bam",GERMLINE="germline.bam",CHR="chr22" \
>> -e script.efile.chr22 \
>> -o script.ofile.chr22 \
>> script.R
>> 
>> Thank you very very much  !
>> 
>> -- bogdan
>> 
>>  [[alternative HTML version deleted]]
>> 
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Efficient swapping

2017-07-07 Thread Anthoni, Peter (IMK)
how about gdata functions?

set.seed(1)
(tmp <- data.frame(x = 1:10, R1 = sample(LETTERS[1:5], 10, replace = TRUE), R2 
= sample(LETTERS[2:6], 10, replace = TRUE)))
tmp.orig=tmp

library(gdata)
bigMap=mapLevels(x=list(factor(tmp[,"R1"]),factor(tmp[,"R2"])),combine = 
T,codes=FALSE)
mapLevels(tmp[,"R1"]) = bigMap
mapLevels(tmp[,"R2"]) = bigMap
dd=as.integer(tmp[,"R2"])-as.integer(tmp[,"R1"])
odd=dd%%2!=0
swap=tmp[odd,"R2"]
tmp[odd,"R2"]=tmp[odd,"R1"]
tmp[odd,"R1"]=swap

cbind(tmp,odd,tmp.orig)

cheers
Peter



> On 7. Jul 2017, at 00:54, Gang Chen  wrote:
> 
> Thanks a lot, Ista! I really appreciate it.
> 
> How about a slightly different case as the following:
> 
> set.seed(1)
> (tmp <- data.frame(x = 1:10, R1 = sample(LETTERS[1:5], 10, replace =
> TRUE), R2 = sample(LETTERS[2:6], 10, replace = TRUE)))
> 
>   x R1 R2
>   1  C  B
>   2  B  B
>   3  C  E
>   4  E  C
>   5  E  B
>   6  D  E
>   7  E  E
>   8  D  F
>   9  C  D
>  10  A  E
> 
> Notice that the factor levels between the two factors, R1 and R2,
> slide by one level; that is, factor R1 does not have level F while
> factor R2 does not have level A. I want to swap the factor levels
> based on the combined levels of the two factors as shown below:
> 
> tl <- unique(c(levels(tmp$R1), levels(tmp$R2)))
> for(ii in 1:dim(tmp)[1]) {
>   kk <- which(tl %in% tmp[ii,'R2'], arr.ind = TRUE) - which(tl %in%
>  tmp[ii,'R1'], arr.ind = TRUE)
>   if(kk%%2!=0) { # swap the their levels between the two factors
>  qq <- tmp[ii,]$R1
>  tmp[ii,]$R1 <- tmp[ii,]$R2
>  tmp[ii,]$R2 <- qq
>  }
> }
> 
> How to go about this case? Thanks!
> 
> 
> On Thu, Jul 6, 2017 at 5:16 PM, Ista Zahn  wrote:
>> How about
>> 
>> foo <- with(list(r1 = tmp$R1,
>> r2 = tmp$R2,
>> swapme = (as.numeric(tmp$R1) - as.numeric(tmp$R2)) %% 2 != 
>> 0),
>> {
>>tmp[swapme, "R1"] <- r2[swapme]
>>tmp[swapme, "R2"] <- r1[swapme]
>>tmp
>> })
>> 
>> Best,
>> Ista
>> 
>> On Thu, Jul 6, 2017 at 4:06 PM, Gang Chen  wrote:
>>> Suppose that we have the following dataframe:
>>> 
>>> set.seed(1)
>>> (tmp <- data.frame(x = 1:10, R1 = sample(LETTERS[1:5], 10, replace =
>>> TRUE), R2 = sample(LETTERS[1:5], 10, replace = TRUE)))
>>> 
>>>x R1 R2
>>> 1   1  B  B
>>> 2   2  B  A
>>> 3   3  C  D
>>> 4   4  E  B
>>> 5   5  B  D
>>> 6   6  E  C
>>> 7   7  E  D
>>> 8   8  D  E
>>> 9   9  D  B
>>> 10 10  A  D
>>> 
>>> I want to do the following: if the difference between the level index
>>> of factor R1 and that of factor R2 is an odd number, the levels of the
>>> two factors need to be switched between them, which can be performed
>>> through the following code:
>>> 
>>> for(ii in 1:dim(tmp)[1]) {
>>>   kk <- which(levels(tmp$R2) %in% tmp[ii,'R2'], arr.ind = TRUE) -
>>> which(levels(tmp$R1) %in% tmp[ii,'R1'], arr.ind = TRUE)
>>>   if(kk%%2!=0) { # swap the their levels between the two factors
>>>  qq <- tmp[ii,]$R1
>>>  tmp[ii,]$R1 <- tmp[ii,]$R2
>>>  tmp[ii,]$R2 <- qq
>>>  }
>>> }
>>> 
>>> More concise and efficient way to do this?
>>> 
>>> Thanks,
>>> Gang
>>> 
>>> __
>>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>>> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] expand gridded matrix to higher resolution

2017-07-05 Thread Anthoni, Peter (IMK)
Hi Jeff,

thanks, the raster package disaggregate will do the trick as well.

library(raster)
rmm <- raster(ncols=5, nrows=3)
rmm[] <- matrix(1:15,nrow=3,byrow = T)
xrmm <- disaggregate(rmm, fact=c(3, 3))
> > as.matrix(rmm)
>  [,1] [,2] [,3] [,4] [,5]
> [1,]12345
> [2,]6789   10
> [3,]   11   12   13   14   15
> > as.matrix(xrmm)
>   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] 
> [,14] [,15]
>  [1,]111222333 4 4 4 5
>  5 5
>  [2,]111222333 4 4 4 5
>  5 5
>  [3,]111222333 4 4 4 5
>  5 5
>  [4,]666777888 9 9 910
> 1010
>  [5,]666777888 9 9 910
> 1010
>  [6,]666777888 9 9 910
> 1010
>  [7,]   11   11   11   12   12   12   13   13   1314141415
> 1515
>  [8,]   11   11   11   12   12   12   13   13   1314141415
> 1515
>  [9,]   11   11   11   12   12   12   13   13   1314141415
> 1515


the disaggregate as a bit faster than the tapply.

mmb=matrix(1:259200,nrow=720,ncol=360)
rmmb <- raster(ncols=360, nrows=720)
rmmb[] <- mmb[]
system.time(for(i in 1:10) {xmm=matrix(NA,nrow=nrow(mmb)*3,ncol=ncol(mmb)*3)
for(icol in 1:ncol(mmb)) {
  for(irow in 1:nrow(mmb)) {
xicol=(icol-1)*3 +c(1:3)
xirow=(irow-1)*3 +c(1:3)
xmm[xirow,xicol]=mmb[irow,icol]
  }
}
})
system.time(for(i in 1:10) {apply(t(apply(mmb,1,rep,each=3)),2,rep,each=3)}) 
#ca. 10x faster
system.time(for(i in 1:10) {xrmmb <- disaggregate(rmmb, fact=c(3, 3))}) 

> > system.time(for(i in 1:10) {xmm=matrix(NA,nrow=nrow(mmb)*3,ncol=ncol(mmb)*3)
> + for(icol in 1:ncol(mmb)) {
> +   for(irow in 1:nrow(mmb)) {
> + xicol=(icol-1)*3 +c(1:3)
> + xirow=(irow-1)*3 +c(1:3)
> + xmm[xirow,xicol]=mmb[irow,icol]
> +   }
> + }
> + })
>user  system elapsed 
>   8.297   0.048   8.364 
> > system.time(for(i in 1:10) 
> > {apply(t(apply(mmb,1,rep,each=3)),2,rep,each=3)}) #ca. 10x faster
>user  system elapsed 
>   0.785   0.093   0.881 
> > system.time(for(i in 1:10) {xrmmb <- disaggregate(rmmb, fact=c(3, 3))})
>user  system elapsed 
>   0.583   0.147   0.731 

cheers
Peter



> On 5. Jul 2017, at 16:57, Jeff Newmiller  wrote:
> 
> You probably ought to be using the raster package. See the CRAN Spatial Task 
> View.
> -- 
> Sent from my phone. Please excuse my brevity.
> 
> On July 5, 2017 12:20:28 AM PDT, "Anthoni, Peter (IMK)" 
>  wrote:
>> Hi all,
>> (if me email goes out as html, than my email client don't do as told,
>> and I apologies already.)
>> 
>> We need to downscale climate data and therefore first need to expand
>> the climate from 0.5deg to the higher resolution 10min, before we can
>> add high resolution deviations. We basically need to have the original
>> data at each gridcell replicated into 3x3 gridcells. 
>> A simple for loop can do this, but I could need a faster procedure.
>> Anybody know a faster way? Is there package than can do what we need
>> already?
>> I tried matrix with rep, but I am missing some magic there, since it
>> doesn't do what we need. 
>> replicate might be promising, but then still need to rearrange the
>> output into the column and row format we need. 
>> 
>> A simple example:
>> mm=matrix(1:15,nrow=3,byrow = T)
>> xmm=matrix(NA,nrow=nrow(mm)*3,ncol=ncol(mm)*3)
>> for(icol in 1:ncol(mm)) {
>> for(irow in 1:nrow(mm)) {
>>   xicol=(icol-1)*3 +c(1:3)
>>   xirow=(irow-1)*3 +c(1:3)
>>   xmm[xirow,xicol]=mm[irow,icol]
>> }
>> }
>> mm
>>>> mm
>>> [,1] [,2] [,3] [,4] [,5]
>>> [1,]12345
>>> [2,]6789   10
>>> [3,]   11   12   13   14   15
>>> 
>> xmm
>>>> xmm
>>>  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
>> [,13] [,14] [,15]
>>> [1,]111222333 4 4 4 
>>  5 5 5
>>> [2,]111222333 4 4 4 
>>  5 5 5
>>> [3,]111222333 4 4 4 
>>  5 5 5
>>> [4,]666777888 9 9 9 
>> 101010
>>> [5,]6667778   

Re: [R] expand gridded matrix to higher resolution

2017-07-05 Thread Anthoni, Peter (IMK)
Hi Jim,

thanks that works like a charm.  

cheers
Peter



> On 5. Jul 2017, at 12:01, Jim Lemon  wrote:
> 
> Hi Peter,
> 
> apply(t(apply(mm,1,rep,each=3)),2,rep,each=3)
> 
> Jim
> 
> On Wed, Jul 5, 2017 at 5:20 PM, Anthoni, Peter (IMK)
>  wrote:
>> Hi all,
>> (if me email goes out as html, than my email client don't do as told, and I 
>> apologies already.)
>> 
>> We need to downscale climate data and therefore first need to expand the 
>> climate from 0.5deg to the higher resolution 10min, before we can add high 
>> resolution deviations. We basically need to have the original data at each 
>> gridcell replicated into 3x3 gridcells.
>> A simple for loop can do this, but I could need a faster procedure. Anybody 
>> know a faster way? Is there package than can do what we need already?
>> I tried matrix with rep, but I am missing some magic there, since it doesn't 
>> do what we need.
>> replicate might be promising, but then still need to rearrange the output 
>> into the column and row format we need.
>> 
>> A simple example:
>> mm=matrix(1:15,nrow=3,byrow = T)
>> xmm=matrix(NA,nrow=nrow(mm)*3,ncol=ncol(mm)*3)
>> for(icol in 1:ncol(mm)) {
>>  for(irow in 1:nrow(mm)) {
>>xicol=(icol-1)*3 +c(1:3)
>>xirow=(irow-1)*3 +c(1:3)
>>xmm[xirow,xicol]=mm[irow,icol]
>>  }
>> }
>> mm
>>>> mm
>>> [,1] [,2] [,3] [,4] [,5]
>>> [1,]12345
>>> [2,]6789   10
>>> [3,]   11   12   13   14   15
>>> 
>> xmm
>>>> xmm
>>>  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] 
>>> [,14] [,15]
>>> [1,]111222333 4 4 4 5   
>>>   5 5
>>> [2,]111222333 4 4 4 5   
>>>   5 5
>>> [3,]111222333 4 4 4 5   
>>>   5 5
>>> [4,]666777888 9 9 910   
>>>  1010
>>> [5,]666777888 9 9 910   
>>>  1010
>>> [6,]666777888 9 9 910   
>>>  1010
>>> [7,]   11   11   11   12   12   12   13   13   1314141415   
>>>  1515
>>> [8,]   11   11   11   12   12   12   13   13   1314141415   
>>>  1515
>>> [9,]   11   11   11   12   12   12   13   13   1314141415   
>>>  1515
>> 
>> I tried various rep with matrix, but don't get the right result.
>> xmm2=matrix(rep(rep(mm,each=3),times=3),nrow=nrow(mm)*3,ncol=ncol(mm)*3,byrow
>>  = F)
>>> identical(xmm,xmm2)
>> [1] FALSE
>> 
>> rr=replicate(3,rep(t(mm),each=3))
>> rr
>>>> rr
>>>  [,1] [,2] [,3]
>>> [1,]111
>>> [2,]111
>>> [3,]111
>>> [4,]222
>>> [5,]222
>>> [6,]222
>>> [7,]333
>>> ...
>> identical(xmm,matrix(rr,ncol=15,nrow=9,byrow=T))
>>>> identical(xmm,matrix(rr,ncol=15,nrow=9,byrow=T))
>>> [1] FALSE
>> 
>> Many thanks for any advice.
>> 
>> cheers
>> Peter
>> 
>> __
>> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
>> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] expand gridded matrix to higher resolution

2017-07-05 Thread Anthoni, Peter (IMK)
Hi all,
(if me email goes out as html, than my email client don't do as told, and I 
apologies already.)

We need to downscale climate data and therefore first need to expand the 
climate from 0.5deg to the higher resolution 10min, before we can add high 
resolution deviations. We basically need to have the original data at each 
gridcell replicated into 3x3 gridcells. 
A simple for loop can do this, but I could need a faster procedure. Anybody 
know a faster way? Is there package than can do what we need already?
I tried matrix with rep, but I am missing some magic there, since it doesn't do 
what we need. 
replicate might be promising, but then still need to rearrange the output into 
the column and row format we need. 

A simple example:
mm=matrix(1:15,nrow=3,byrow = T)
xmm=matrix(NA,nrow=nrow(mm)*3,ncol=ncol(mm)*3)
for(icol in 1:ncol(mm)) {
  for(irow in 1:nrow(mm)) {
xicol=(icol-1)*3 +c(1:3)
xirow=(irow-1)*3 +c(1:3)
xmm[xirow,xicol]=mm[irow,icol]
  }
}
mm
> > mm
>  [,1] [,2] [,3] [,4] [,5]
> [1,]12345
> [2,]6789   10
> [3,]   11   12   13   14   15
> 
xmm
> > xmm
>   [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] 
> [,14] [,15]
>  [1,]111222333 4 4 4 5
>  5 5
>  [2,]111222333 4 4 4 5
>  5 5
>  [3,]111222333 4 4 4 5
>  5 5
>  [4,]666777888 9 9 910
> 1010
>  [5,]666777888 9 9 910
> 1010
>  [6,]666777888 9 9 910
> 1010
>  [7,]   11   11   11   12   12   12   13   13   1314141415
> 1515
>  [8,]   11   11   11   12   12   12   13   13   1314141415
> 1515
>  [9,]   11   11   11   12   12   12   13   13   1314141415
> 1515

I tried various rep with matrix, but don't get the right result.
xmm2=matrix(rep(rep(mm,each=3),times=3),nrow=nrow(mm)*3,ncol=ncol(mm)*3,byrow = 
F)
> identical(xmm,xmm2)
[1] FALSE

rr=replicate(3,rep(t(mm),each=3))
rr
> > rr
>   [,1] [,2] [,3]
>  [1,]111
>  [2,]111
>  [3,]111
>  [4,]222
>  [5,]222
>  [6,]222
>  [7,]333
> ...
identical(xmm,matrix(rr,ncol=15,nrow=9,byrow=T))
> > identical(xmm,matrix(rr,ncol=15,nrow=9,byrow=T))
> [1] FALSE
 
Many thanks for any advice.

cheers
Peter

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Help: ifelse selection for x,y coordinates

2017-06-22 Thread Anthoni, Peter (IMK)
Hi Celine,

what about removing the unwanted after you made the x and y
x<-x[x>0]  # or x<-x[x>0&&y>0], ditto for y, x[x!=""] in your ifelse (... ,"") 
case

if x and y will not have the same length afterwards you need to make that list 
thingy.

cheers
Peter



On 23. Jun 2017, at 07:30, Céline Lüscher 
mailto:c-luesc...@hispeed.ch>> wrote:

Hi Jim,

Thank you very much for the answer ! The result is really better with this 😊

Here is the code :
kk<- function(x.Koordinate, y.Koordinate, data=data)
+ {
+   coordx<-data$x.Koordinate[data$G==24]
+   coordy<-data$y.Koordinate[data$G==24]
+   x <- ifelse(data$x.Koordinate>coordx-51 & data$G>15,data$x.Koordinate," ")
+   y<-ifelse(data$y.Koordinate>coordy-51 & data$G>15,data$y.Koordinate," ")
+   xy<-as.data.frame(list(x,y))
+   names(xy)<-c("x","y")
+   return(xy)
+ }
kk(x.Koordinate, y.Koordinate, data=data)
   x  y
1
2
3
4
5
6  205550
7  205550 604100
8
9  205600 604150
10 205600 604100

Best regards,
C.


Gesendet von Mail für Windows 10

Von: Jim Lemon
Gesendet: vendredi, 23 juin 2017 05:28
An: Céline Lüscher
Cc: r-help@r-project.org
Betreff: Re: [R] Help: ifelse selection for x,y coordinates

Hi Celine,
Perhaps if you modify your return value like this:

xy<-as.data.frame(list(x,y))
names(xy)<-c("x","y")
return(xy)

Jim

On Thu, Jun 22, 2017 at 6:48 PM, Céline Lüscher 
mailto:c-luesc...@hispeed.ch>> wrote:
Hi everyone,
My database has 3 columns : the x coordinates, the y coordinates and the value 
of G for every individual (20 in total). I would like to have the x,y 
coordinates of individuals, Under these conditions : the distance to the 
reference coordinates must be under or equal to 50 meters + the value of G must 
be bigger or equal to 16.

Here’s what I’ve done first :

kk<- function(x, y)
+ {
+   coordx<-data$x.Koordinate[data$G==24]
+   coordy<-data$y.Koordinate[data$G==24]
+   x <- ifelse(data$x.Koordinate>coordx-51 & data$G>15,data$x.Koordinate,0)
+   y<-ifelse(data$y.Koordinate>coordy-51 & data$G>15,data$y.Koordinate,0)
+   return(c(x,y))
+ }
kk(data$x.Koordinate, data$y.Koordinate)
[1]  0  0  0  0  0 205550 205550  0 205600 205600  
0  0  0  0  0  0  0
[18] 604100  0 604150 604100  0

The problem here is that we can not clearly see the difference between the 
coordinates for x and the ones for y.

Then I tried this :
kk<- function(x, y)
+ {
+   coordx<-data$x.Koordinate[data$G==24]
+   coordy<-data$y.Koordinate[data$G==24]
+   x <- ifelse(data$x.Koordinate>coordx-51 & data$G>15,data$x.Koordinate," ")
+   y<-ifelse(data$y.Koordinate>coordy-51 & data$G>15,data$y.Koordinate," ")
+   return(list(x,y))
+ }
kk(data$x.Koordinate, data$y.Koordinate)
[[1]]
[1] " "  " "  " "  " "  " "  "205550" "205550" " "  
"205600" "205600" " "

[[2]]
[1] " "  " "  " "  " "  " "  " "  "604100" " "  
"604150" "604100" " "



Where we can see better the two levels related to the x and y coordinates.

My question is simple : Is it possible for this function to return the values 
in a form like x,y or x y ? (without any 0, or « », or space) Or should I use 
another R function to obtain this result ?

Thank you for your help, and sorry for the English mistakes,
C.


Gesendet von Mail für Windows 10


   [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To 
UNSUBSCRIBE and more, see
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.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To 
UNSUBSCRIBE and more, see
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.


[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] memory problem

2017-05-02 Thread Anthoni, Peter (IMK)
Hi Amit,

Is the file gzipped or extracted?
if you read the plain text file, try to gzip it and make a read.table on the 
gzipped file, the read.table can handle gzipped files at least on linux and mac 
OS, not sure about windows.

cheers
Peter



> On 2. May 2017, at 18:59, Amit Sengupta via R-help  
> wrote:
> 
> Hi,I was unable to read a 2.4 gig file into an R object using read.table in 
> 64 bit R environment. Please let me have your suggestions.Amit Sengupta
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] asking for help

2017-04-25 Thread Anthoni, Peter (IMK)
Hi,

the cut function might be helpful.

vec=1: 163863
fcut=cut(vec,seq(1, 163863+1,by= 6069),include.lowest = T,right=F)
aggregate(vec,by=list(fcut),min)
aggregate(vec,by=list(fcut),max)

cheers
Peter



On 25. Apr 2017, at 14:33, PIKAL Petr 
mailto:petr.pi...@precheza.cz>> wrote:

Hi

-Original Message-
From: R-help [mailto:r-help-boun...@r-project.org] On Behalf Of Thomas
Mailund
Sent: Tuesday, April 25, 2017 11:20 AM
To: r-help@r-project.org; Saifuddin, Miah Mohammad
mailto:miah.mohammadsaif...@mavs.uta.edu>>
Subject: Re: [R] asking for help

If you write something like

indices <- rep(1:(163863/6069), each = 6069)

You can get similar result by

indices2 <- 0:163862%/%6069

but starting with zero.

or the same with
indices2 <- (0:163862%/%6069)+1

Cheers
Petr


you can get the i’th block of rows with

table[indices == i,]

It looks like a little more work than a loop would be since you have to run
through all rows for each block, but the implicit loop in this approach is 
likely
to be faster than an explicit for-loop.

Cheers
 Thomas

On 25 Apr 2017, 07.01 +0200, Saifuddin, Miah Mohammad
mailto:miah.mohammadsaif...@mavs.uta.edu>>, 
wrote:
I have a data frame having 163863 values. I want to subset it so that each set
has 6069 values in it. for example 1:6069 as first, 6070: 6070+6068 as second.
how can I do that, preferably in a loop.


TIA

[[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To 
UNSUBSCRIBE and more, see
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.

 [[alternative HTML version deleted]]

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.


Tento e-mail a jakékoliv k němu připojené dokumenty jsou důvěrné a jsou určeny 
pouze jeho adresátům.
Jestliže jste obdržel(a) tento e-mail omylem, informujte laskavě neprodleně 
jeho odesílatele. Obsah tohoto emailu i s přílohami a jeho kopie vymažte ze 
svého systému.
Nejste-li zamýšleným adresátem tohoto emailu, nejste oprávněni tento email 
jakkoliv užívat, rozšiřovat, kopírovat či zveřejňovat.
Odesílatel e-mailu neodpovídá za eventuální škodu způsobenou modifikacemi či 
zpožděním přenosu e-mailu.

V případě, že je tento e-mail součástí obchodního jednání:
- vyhrazuje si odesílatel právo ukončit kdykoliv jednání o uzavření smlouvy, a 
to z jakéhokoliv důvodu i bez uvedení důvodu.
- a obsahuje-li nabídku, je adresát oprávněn nabídku bezodkladně přijmout; 
Odesílatel tohoto e-mailu (nabídky) vylučuje přijetí nabídky ze strany příjemce 
s dodatkem či odchylkou.
- trvá odesílatel na tom, že příslušná smlouva je uzavřena teprve výslovným 
dosažením shody na všech jejích náležitostech.
- odesílatel tohoto emailu informuje, že není oprávněn uzavírat za společnost 
žádné smlouvy s výjimkou případů, kdy k tomu byl písemně zmocněn nebo písemně 
pověřen a takové pověření nebo plná moc byly adresátovi tohoto emailu případně 
osobě, kterou adresát zastupuje, předloženy nebo jejich existence je adresátovi 
či osobě jím zastoupené známá.

This e-mail and any documents attached to it may be confidential and are 
intended only for its intended recipients.
If you received this e-mail by mistake, please immediately inform its sender. 
Delete the contents of this e-mail with all attachments and its copies from 
your system.
If you are not the intended recipient of this e-mail, you are not authorized to 
use, disseminate, copy or disclose this e-mail in any manner.
The sender of this e-mail shall not be liable for any possible damage caused by 
modifications of the e-mail or by delay with transfer of the email.

In case that this e-mail forms part of business dealings:
- the sender reserves the right to end negotiations about entering into a 
contract in any time, for any reason, and without stating any reasoning.
- if the e-mail contains an offer, the recipient is entitled to immediately 
accept such offer; The sender of this e-mail (offer) excludes any acceptance of 
the offer on the part of the recipient containing any amendment or variation.
- the sender insists on that the respective contract is concluded only upon an 
express mutual agreement on all its aspects.
- the sender of this e-mail informs that he/she is not authorized to enter into 
any contracts on behalf of the company except for cases in which he/she is 
expressly authorized to do so in writing, and such authorization or power of 
attorney is submitted to the recipient or the person represented by the 
recipient, or the existence of such au

Re: [R] Looping Through DataFrames with Differing Lenghts

2017-03-27 Thread Anthoni, Peter (IMK)
Hi Paul,

match might help, but without a real data sample, it is hard to check if the 
following might work.

mm=match(df.col378[,"Date"],df.col362[,"Date"])
#mm will have NAs, where there is no matching date in df.col362
#and have the index of the match, where the two dates match
new.df=cbind(df.col378,"transits.col362"=df.col362[mm,"transits"])

cheers
Peter



> On 27 Mar 2017, at 22:09, Paul Bernal  wrote:
> 
> Dear friends,
> 
> I have one dataframe which contains 378 observations, and another one,
> containing 362 observations.
> 
> Both dataframes have two columns, one date column and another one with the
> number of transits.
> 
> I wanted to come up with a code so that I could fill in the dates that are
> missing in one of the dataframes and replace the column of transits with
> the value NA.
> 
> I have tried several things but R obviously complains that the length of
> the dataframes are different.
> 
> How can I solve this?
> 
> Any guidance will be greatly appreciated,
> 
> Best regards,
> 
> Paul
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Aggregate matrix in a 2 by 2 manor

2016-07-31 Thread Anthoni, Peter (IMK)
 v, ncol = nnc )
> }
> test1 <- aggregate.nx.ny.expand.grid( tst )
> test2 <- aggregate.nx.ny.array.apply( tst )
> test3 <- aggregate.nx.ny.array.aperm.mean( tst )
> test4 <- aggregate.nx.ny.array.aperm.apply( tst )
> library(microbenchmark)
> microbenchmark(
>  aggregate.nx.ny.expand.grid( tst, 2, 2, mean, na.rm = TRUE )
> , aggregate.nx.ny.array.apply( tst, 2, 2, mean, na.rm = TRUE )
> , aggregate.nx.ny.array.aperm.mean( tst, 2, 2, na.rm = TRUE )
> , aggregate.nx.ny.array.aperm.apply( tst, 2, 2, mean, na.rm = TRUE )
> )
> #Unit: milliseconds
> # exprmin
> #   aggregate.nx.ny.expand.grid(tst, 2, 2, mean, na.rm = TRUE) 628.528322
> #   aggregate.nx.ny.array.apply(tst, 2, 2, mean, na.rm = TRUE) 846.883314
> #aggregate.nx.ny.array.aperm.mean(tst, 2, 2, na.rm = TRUE)   8.904369
> # aggregate.nx.ny.array.aperm.apply(tst, 2, 2, mean, na.rm = TRUE) 619.691851
> # lq   mean medianuq  max neval cld
> # 675.470967  916.39630  778.54090  873.9754 2452.695   100  b
> # 920.831966 1126.94691 1000.33830 1094.9233 3412.639   100   c
> #   9.191747   21.98528   10.30099   15.9169  158.687   100 a
> # 733.246331  936.73359  757.58383  844.2016 2824.557   100  b
> 
> 
> On Sat, 30 Jul 2016, Jeff Newmiller wrote:
> 
>> For the record, the array.apply code can be fixed as below, but then it is 
>> slower than the expand.grid version.
>> 
>> aggregate.nx.ny.array.apply <- function(dta,nx=2,ny=2, FUN=mean,...)
>> {
>>  a <- array(dta, dim = c(ny, nrow( dta ) %/% ny, nx, ncol( dta ) %/% nx))
>> apply( a, c(2, 4), FUN, ... )
>> }
>> 
>> -- 
>> Sent from my phone. Please excuse my brevity.
>> 
>> On July 30, 2016 11:06:16 AM PDT, "Anthoni, Peter (IMK)" 
>>  wrote:
>>> Hi all,
>>> 
>>> thanks for the suggestions, I did some timing tests, see below.
>>> Unfortunately the aggregate.nx.ny.array.apply, does not produce the
>>> expected result.
>>> So the fastest seems to be the aggregate.nx.ny.expand.grid, though the
>>> double for loop is not that much slower.
>>> 
>>> many thanks
>>> Peter
>>> 
>>>> tst=matrix(1:(1440*360),ncol=1440,nrow=360)
>>>> system.time( {for(i in 1:10)
>>> tst_2x2=aggregate.nx.ny.forloop(tst,2,2,mean,na.rm=T)})
>>>  user  system elapsed
>>> 11.227   0.073  11.371
>>>> system.time( {for(i in 1:10)
>>> tst_2x2=aggregate.nx.ny.interaction(tst,2,2,mean,na.rm=T)})
>>>  user  system elapsed
>>> 26.354   0.475  26.880
>>>> system.time( {for(i in 1:10)
>>> tst_2x2=aggregate.nx.ny.expand.grid(tst,2,2,mean,na.rm=T)})
>>>  user  system elapsed
>>> 9.683   0.055   9.763
>>>> system.time( {for(i in 1:10)
>>> tst_2x2=aggregate.nx.ny.array.apply(tst,2,2,mean,na.rm=T)})
>>>  user  system elapsed
>>> 7.693   0.055   7.800
>>> 
>>>> tst.small=matrix(1:(8*4),ncol=8,nrow=4)
>>>> aggregate.nx.ny.forloop = function(data,nx=2,ny=2, FUN=mean,...)
>>> + {
>>> +   nlon=nrow(data)
>>> +   nlat=ncol(data)
>>> +   newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
>>> +   dim(newdata)
>>> +   for(ilon in seq(1,nlon,nx)) {
>>> + for(ilat in seq(1,nlat,ny)) {
>>> +   ilon_new=1+(ilon-1)/nx
>>> +   ilat_new=1+(ilat-1)/ny
>>> +   newdata[ilon_new,ilat_new] = FUN(data[ilon+0:1,ilat+0:1],...)
>>> + }
>>> +   }
>>> +   newdata
>>> + }
>>>> aggregate.nx.ny.forloop(tst.small)
>>>[,1] [,2] [,3] [,4]
>>> [1,]  3.5 11.5 19.5 27.5
>>> [2,]  5.5 13.5 21.5 29.5
>>>> 
>>>> aggregate.nx.ny.interaction = function(data,nx=2,ny=2, FUN=mean,...)
>>> + {
>>> +
>>> +   nlon=nrow(data)
>>> +   nlat=ncol(data)
>>> +   newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
>>> +   newdata[] <- tapply( data, interaction( (row(data)+1) %/% 2,
>>> (col(data)+1) %/% 2 ), FUN, ...)
>>> +   newdata
>>> + }
>>>> aggregate.nx.ny.interaction(tst.small)
>>>[,1] [,2] [,3] [,4]
>>> [1,]  3.5 11.5 19.5 27.5
>>> [2,]  5.5 13.5 21.5 29.5
>>>> 
>>>> aggregate.nx.ny.expand.grid = function(data,nx=2,ny=2, FUN=mean,...)
>>> + {
>>> +   ilon <- seq(1,ncol(data),nx)
>>> +   ilat <- seq(1,nrow(data),ny)
>>> +   cells <- as.matrix(expand.grid(ilat, ilon))
>>> +   blocks <- apply(cells, 1, function(x)
>>> data[x[1]:(x[1]

Re: [R] Aggregate matrix in a 2 by 2 manor

2016-07-30 Thread Anthoni, Peter (IMK)
Hi all,

thanks for the suggestions, I did some timing tests, see below.
Unfortunately the aggregate.nx.ny.array.apply, does not produce the expected 
result.
So the fastest seems to be the aggregate.nx.ny.expand.grid, though the double 
for loop is not that much slower.

many thanks
Peter

> tst=matrix(1:(1440*360),ncol=1440,nrow=360)
> system.time( {for(i in 1:10) 
> tst_2x2=aggregate.nx.ny.forloop(tst,2,2,mean,na.rm=T)})
   user  system elapsed 
 11.227   0.073  11.371 
> system.time( {for(i in 1:10) 
> tst_2x2=aggregate.nx.ny.interaction(tst,2,2,mean,na.rm=T)})
   user  system elapsed 
 26.354   0.475  26.880 
> system.time( {for(i in 1:10) 
> tst_2x2=aggregate.nx.ny.expand.grid(tst,2,2,mean,na.rm=T)})
   user  system elapsed 
  9.683   0.055   9.763 
> system.time( {for(i in 1:10) 
> tst_2x2=aggregate.nx.ny.array.apply(tst,2,2,mean,na.rm=T)})
   user  system elapsed 
  7.693   0.055   7.800 

> tst.small=matrix(1:(8*4),ncol=8,nrow=4)
> aggregate.nx.ny.forloop = function(data,nx=2,ny=2, FUN=mean,...) 
+ {
+   nlon=nrow(data)
+   nlat=ncol(data)
+   newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
+   dim(newdata)
+   for(ilon in seq(1,nlon,nx)) {
+ for(ilat in seq(1,nlat,ny)) {
+   ilon_new=1+(ilon-1)/nx
+   ilat_new=1+(ilat-1)/ny
+   newdata[ilon_new,ilat_new] = FUN(data[ilon+0:1,ilat+0:1],...)
+ }
+   }
+   newdata
+ }
> aggregate.nx.ny.forloop(tst.small)
 [,1] [,2] [,3] [,4]
[1,]  3.5 11.5 19.5 27.5
[2,]  5.5 13.5 21.5 29.5
> 
> aggregate.nx.ny.interaction = function(data,nx=2,ny=2, FUN=mean,...) 
+ {
+   
+   nlon=nrow(data)
+   nlat=ncol(data)
+   newdata=matrix(NA,nrow=nlon/nx,ncol=nlat/ny)
+   newdata[] <- tapply( data, interaction( (row(data)+1) %/% 2, (col(data)+1) 
%/% 2 ), FUN, ...)
+   newdata
+ }
> aggregate.nx.ny.interaction(tst.small)
 [,1] [,2] [,3] [,4]
[1,]  3.5 11.5 19.5 27.5
[2,]  5.5 13.5 21.5 29.5
> 
> aggregate.nx.ny.expand.grid = function(data,nx=2,ny=2, FUN=mean,...) 
+ {
+   ilon <- seq(1,ncol(data),nx)
+   ilat <- seq(1,nrow(data),ny)
+   cells <- as.matrix(expand.grid(ilat, ilon))
+   blocks <- apply(cells, 1, function(x) data[x[1]:(x[1]+1),x[2]:(x[2]+1)])
+   block.means <- colMeans(blocks)
+   matrix(block.means, nrow(data)/ny, ncol(data)/nx)
+ }
> aggregate.nx.ny.expand.grid(tst.small)
 [,1] [,2] [,3] [,4]
[1,]  3.5 11.5 19.5 27.5
[2,]  5.5 13.5 21.5 29.5
> 
> aggregate.nx.ny.array.apply = function(data,nx=2,ny=2, FUN=mean,...) {
+   a <- array(data, dim = c(ny, nrow( data ) %/% ny, ncol( data ) %/% nx))
+   apply( a, c(2, 3), FUN, ... )
+ }
> aggregate.nx.ny.array.apply(tst.small)
 [,1] [,2] [,3] [,4]
[1,]  1.5  5.5  9.5 13.5
[2,]  3.5  7.5 11.5 15.5



> On 28 Jul 2016, at 00:26, David Winsemius  wrote:
> 
> 
>> On Jul 27, 2016, at 12:02 PM, Jeff Newmiller  
>> wrote:
>> 
>> An alternative (more compact, not necessarily faster, because apply is still 
>> a for loop inside):
>> 
>> f <- function( m, nx, ny ) {
>> # redefine the dimensions of my
>> a <- array( m
>>, dim = c( ny
>>   , nrow( m ) %/% ny
>>   , ncol( m ) %/% nx )
>>   )
>> # apply mean over dim 1
>> apply( a, c( 2, 3 ), FUN=mean )
>> }
>> f( tst, nx, ny )
> 
> Here's an apparently loopless strategy, although I suspect the code for 
> interaction (and maybe tapply as well?) uses a loop.
> 
> 
> tst_2X2 <- matrix(NA, ,ncol=4,nrow=2)
> 
> tst_2x2[] <- tapply( tst, interaction( (row(tst)+1) %/% 2, (col(tst)+1) %/% 2 
> ), mean)
> 
> tst_2x2
> 
> [,1] [,2] [,3] [,4]
> [1,]  3.5 11.5 19.5 27.5
> [2,]  5.5 13.5 21.5 29.5
> 
> -- 
> David.
> 
> 
>> 
>> -- 
>> Sent from my phone. Please excuse my brevity.
>> 
>> On July 27, 2016 9:08:32 AM PDT, David L Carlson  wrote:
>>> This should be faster. It uses apply() across the blocks. 
>>> 
>>>> ilon <- seq(1,8,nx)
>>>> ilat <- seq(1,4,ny)
>>>> cells <- as.matrix(expand.grid(ilat, ilon))
>>>> blocks <- apply(cells, 1, function(x) tst[x[1]:(x[1]+1),
>>> x[2]:(x[2]+1)])
>>>> block.means <- colMeans(blocks)
>>>> tst_2x2 <- matrix(block.means, 2, 4)
>>>> tst_2x2
>>>   [,1] [,2] [,3] [,4]
>>> [1,]  3.5 11.5 19.5 27.5
>>> [2,]  5.5 13.5 21.5 29.5
>>> 
>>> -
>>> David L Carlson
>>> Department of Anthropology
>>> Texas A&M University
>>> College Station, TX 77840-4352
>>> 
>>> 
>>> 
>>> -Original Message-
>>> From: R-help [mailto:r-help-boun...@r-poject.org] On Behalf Of Anthoni,
>>> Peter 

[R] Aggregate matrix in a 2 by 2 manor

2016-07-27 Thread Anthoni, Peter (IMK)
Hi all,

I need to aggregate some matrix data (1440x720) to a lower dimension (720x360) 
for lots of years and variables

I can do double for loop, but that will be slow. Anybody know a quicker way?

here an example with a smaller matrix size:

tst=matrix(1:(8*4),ncol=8,nrow=4)
tst_2x2=matrix(NA,ncol=4,nrow=2)
nx=2
ny=2
for(ilon in seq(1,8,nx)) {
  for (ilat in seq(1,4,ny)) {
ilon_2x2=1+(ilon-1)/nx
ilat_2x2=1+(ilat-1)/ny
tst_2x2[ilat_2x2,ilon_2x2] = mean(tst[ilat+0:1,ilon+0:1])
  }
}

tst
tst_2x2

> tst
 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]159   13   17   21   25   29
[2,]26   10   14   18   22   26   30
[3,]37   11   15   19   23   27   31
[4,]48   12   16   20   24   28   32

> tst_2x2
 [,1] [,2] [,3] [,4]
[1,]  3.5 11.5 19.5 27.5
[2,]  5.5 13.5 21.5 29.5


I though a cast to 3d-array might do the trick and apply over the new 
dimension, but that does not work, since it casts the data along the row.
> matrix(apply(array(tst,dim=c(nx,ny,8)),3,mean),nrow=nrow(tst)/ny)
 [,1] [,2] [,3] [,4]
[1,]  2.5 10.5 18.5 26.5
[2,]  6.5 14.5 22.5 30.5


cheers
Peter

__
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Extracting point data using longitude and latitude from netcdf file using R

2016-01-09 Thread Anthoni, Peter (IMK)
Hi Peter,

the start in nc_varget requires a latitude and longitude index, not the 
latitude and longitude in double format.
So you need to figure out what index your latitude and longitude correspond to, 
which will depends on what data are in your netCDF.

it might have looked like that it worked for a positive latitude, but you got 
the data from the latitude index 6 or 7, depends on how the double was 
transformed into an integer.

best regards
Peter

> On 09 Jan 2016, at 12:28, Peter Tuju via R-help  wrote:
> 
> I have data file in netcdf with three dimensions (x, y, t) and I want to 
> extract a variable RAINC and RAINNC 
> using longitude and latitude for a single point location with all the time, 
> but no lucky. The syntax is as follows;;
> setwd( "/run/media/tuju/0767090047/extract_wrf_txt_file" )
> rm( list = ls() )  
> library( ncdf4 ) 
> inp_file <- nc_open( "wrfout_d01_2016-01-07.nc" )time <- ncvar_get( inp_file, 
> "Times" )  # Reading the time variabledar_lon <- 
> 39.2dar_lat <- -6.87
> RAINC <- ncvar_get( inp_file, varid = "RAINC", start  = c( dar_lon, dar_lat, 
> 1 ), count = c( 1, 1, -1 ) )
> RAINNC <- ncvar_get( inp_file, varid = "RAINNC", start  =  c( dar_lon, 
> dar_lat, 1 ), count = c( 1, 1, -1 ) )
> RAIN <- RAINC + RAINNCRAIN_TABLE <- cbind( time, RAIN )
> write.table( RAIN_TABLE, "Dar_es_Salaam.txt", row.names = FALSE, 
>  col.names = c( "Valid Forecast Time",  "Rain (mm)", sep = "\t " )
> 
> # But no lucky with the red bolded syntax as I end up with the following 
> error message> RAINC <- ncvar_get( inp_file, varid = "RAINC", start  = c( 
> Lon[2], Lat[2], 1 ), count = c( 1, 1, -1 ) )
> Error in Rsx_nc4_get_vara_double: NetCDF: Index exceeds dimension bound
> Var: RAINC  Ndims: 3   Start: 0,4294967289,38 Count: 17,1,1
> Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, 
> addOffset,  : 
>   C function R_nc4_get_vara_double returned error
> However when I cahnge the latitude to postive it works fine. Note latitudes 
> in the file data ranges from -16.71505 to 7.787529 as shown below;
> head(ncvar_get(inp_file, "XLAT"))
> [1] -16.71505 -16.71505 -16.71505 -16.71505 -16.71505 -16.71505
>> tail(ncvar_get(inp_file, "XLAT"))
> [1] 7.787529 7.787529 7.787529 7.787529 7.787529 7.787529
> ## So, how can I get the syntax correct? Please help _
> Peter  E. Tuju
> Dar es Salaam
> T A N Z A N I A
> --
> 
> 
> 
> 
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] assigning values to elements of matrixes

2015-12-23 Thread Anthoni, Peter (IMK)
Hi,

How about the following:
foo2 <- function(s,i,j,value)
{
  M = get(paste("M_",s,sep=""))
  M[i,j] = value
  assign(paste("M_",s,sep=""),M, envir = .GlobalEnv) 
}

foo2("a",1,2,15)

cheers
Peter

> On 23 Dec 2015, at 09:44, Matteo Richiardi  wrote:
> 
> I am following the example I find on ?assign:
> 
> a <- 1:4
> assign("a[1]", 2)
> 
> This appears to create a new variable named "a[1]" rather than
> changing the value of the vector.
> 
> Am I missing something here? How can I assign a value to a specified
> element of a vector/matrix?
> 
> Of course, my problem is slightly more involved, but I guess the above
> is its core. For those interested, I have the two matrixes
> M_a <- matrix(0,10,10)
> M_b <- matrix(0,10,10)
> 
> I want to have a function that assigns a value to a specific element
> of one of the matrix, as in
> foo <- function(s,i,j,value){
>  assign(paste("M_",s)[i,j],value)
> }
> 
> This however does not work:
> 
> foo('a',1,1,1)
> 
> Error in paste("M_", s)[1, j] : incorrect number of dimensions
> 
> Following the ?assign help, I tried
> 
> foo2 <- function(s,i,j,value){
>  assign(paste("M_",s,"[i,j]"),value, envir = .GlobalEnv)
> }
> 
> but this produces the problem I described above (namely, a new
> variable is created, rather than replacing the specified element of
> the matrix).
> 
> I know that in this example I could simply pass the matrix as an
> argument to the function, but I am interested in understanding how
> 'assign' work.
> 
> Many thanks for your help.
> 
> Matteo
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Multiple if function

2015-09-15 Thread Anthoni, Peter (IMK)
Hi,

I guess this might work too and might be quite speedy:

ASBclass = factor(c(1,2,2,3,2,1))
Flow = c(1,1,1,1,1,1)

mult = ((ASBclass==1) * 0.1 + (ASBclass==2) * 0.15 + (ASBclass==3) * 0.2)
deviation = mult * Flow

or with the more complex arithmetic:

deviation = ((ASBclass==1) * (Flow*2) + (ASBclass==2) * (Flow+3) + 
(ASBclass==3) * sqrt(Flow))

cheers
Peter



> On 16 Sep 2015, at 04:20, Charles C. Berry  wrote:
> 
> On Tue, 15 Sep 2015, Bert Gunter wrote:
> 
>> Thanks to both Davids.
>> 
>> I realize that these things are often a matter of aesthetics -- and
>> hence have little rational justification -- but I agree with The Other
>> David: eval(parse) seems to me to violate R's soul( it makes R a macro
>> language instead of a functional one).
>> 
>> However, mapply(... switch) effectively loops through the frame row by
>> row. Aesthetically, I like it; but it seems inefficient. If there are
>> e.g. 1e6 rows in say 10 categories, I think Jeff's approach should do
>> much better.  I'll try to generate some actual data to see unless
>> someone else beats me to it.
> 
> Use mapply like this on large problems:
> 
> unsplit(
>mapply(
>function(x,z) eval( x, list( y=z )),
>expression( A=y*2, B=y+3, C=sqrt(y) ),
>split( dat$Flow, dat$ASB ),
>SIMPLIFY=FALSE),
>dat$ASB)
> 
> Chuck
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Multiple if function

2015-09-15 Thread Anthoni, Peter (IMK)
Hi Maria, 

Why not exploit some simple boolean facts (FALSE=0, TRUE=1) in your 
calculation, so

ASBclass = c(1,2,2,3,2,1)
Flow = c(1,1,1,1,1,1)

factor = ((ASBclass==1) * 0.1 + (ASBclass==2) * 0.15 + (ASBclass==3) * 0.2)
deviation = factor * Flow

cheers
Peter



> On 15 Sep 2015, at 12:56, Maria Lathouri  wrote:
> 
> Dear all, 
> 
> I am writing as I would like your help. I have a dataframe with two columns, 
> ASB and Flow, where the the first one has values 1, 2 or 3 and the second 
> flow data. Something like that:
> ASBclassFlow1  11.51   9.2
> 2  10.5
> 3   6.7  ...  ...
> I would like to produce a third column named eg. deviation where it would get 
> me values based on if ASBclass is 1, multiply Flow by 0.1; if ASBclass is 2 
> then multiply Flow by 0.15 and if ASBclass is 3 then multiply by 0.2.
> 
> If (ASBclass=1) { deviation<-Flow*0.1}
> If (ASBclass=2) { deviation<-Flow*0.15}If (ASBclass=1) { deviation<-Flow*0.2}
> I am not sure whether I should add the else function and how can I combine 
> these separate functions.
> 
> Can anyone help me on that?
> Thank you very much. 
> 
> Kind regardsMaria
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Error -> cannot open file 'specdata/001.csv': No such file or directory; Windows 8, R Version 3.2.1

2015-08-16 Thread Anthoni, Peter (IMK)
Hi Nikita,

To check whether the files are really there, run the following at the R prompt:

getwd()
list.files()

> *C:/Users/acer/My Documents/specdata/rprog-data-specdata/specdata*

Your working path looks like you either need to set it to:
  C:/Users/acer/My Documents/specdata/rprog-data-specdata
or change your read.csv to:
   df <- read.csv("001.csv")

cheers
Peter



> On 15 Aug 2015, at 19:06, Nikita Dinger  wrote:
> 
> I am having a problem in opening the excel files in specdata folder.
> 
> I have completed coding the R program for the assignment but when I run the
> following commands in the R console,
> 
> *source("pollutantmean.R")*
> *> pollutantmean("specdata", "nitrate", 23)*
> 
> I get an error message stating
> 
> *Error in file(file, "rt") : cannot open the connection*
> *In addition: Warning message:*
> *In file(file, "rt") :*
> *  cannot open file 'specdata/023.csv': No such file or directory*
> 
> I tried everything and reset my Working Directory to
> 
> 
> *C:/Users/acer/My Documents/specdata/rprog-data-specdata/specdata*
> 
> After the last specdata folder are all the excel sheets numbered 001 to 332.
> 
> I searched the internet and all other options available, and got a solution
> to open it using the following command:
> 
> *df <- read.csv("specdata/001.csv")*
> 
> 
> This generated the following error message
> 
> *Error in file(file, "rt") : cannot open the connection*
> *In addition: Warning message:*
> *In file(file, "rt") :*
> *  cannot open file 'specdata/001.csv': No such file or directory*
> 
> I have tried various other commands also such as
> 
> 
> *path <- c(paste("./",directory, "/",formatC(id[i], width=3,
> flag=0),".csv",sep=""))*
> 
> However, all the commands show some error.
> 
> What shall I do?
> I am using the 3.2.1 version of R on a Windows 8 laptop.
> 
> Regards,
> Nikita Dinger
> 
>   [[alternative HTML version deleted]]
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Varying name of output tables from looped process of list of spdf objects

2015-07-27 Thread Anthoni, Peter (IMK)
Hi Cecilia,

print(fnp) and print(a_10) are not exactly the same

>> print(fnp)
> $a_10
> class   : SpatialPolygonsDataFrame 
vs.
>> print(a_10)
> class   : SpatialPolygonsDataFrame 


Looks like mget returns a list.

Can you try: 
fnp <- get(afiles[ifile])[1]
print(fnp)
#the $a_10 should be gone.

cheers
Peter




> On 27 Jul 2015, at 18:33, SisoL  wrote:
> 
> Hi Peter,
> 
> Thank you for your reply. The method for looping seems to work, but
> gDistance will not recognise the input. I am puzzled because when I
> print(fnp), and print (a_10) they look exactly the same, but when I try to
> run gDistance{rgeos} with a_10 it works, but with fnp it throws an error.
> Please see below. Any other suggestions?
> 
>> afiles <- ls(pattern= "a_") #OK
>> print(afiles)
> [1] "a_1"  "a_10"
> 
>> fnp<- mget(afiles[ifile]) #OK
> 
>> print(fnp)
> $a_10
> class   : SpatialPolygonsDataFrame 
> features: 5 
> extent  : 825796.9, 831270.1, 815666.9, 816562.5  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0
> +y_0=0 +ellps=aust_SA +units=m +no_defs 
> variables   : 5
> names   : ID, GRIDCODE,Shape_Leng,Shape_Area, Count 
> min values  :  1,1, 94.2722341164, 6.47354525991, 1 
> max values  :  5,1, 2305.45647624, 92111.8528756, 1 
> 
>> print(a_10)
> class   : SpatialPolygonsDataFrame 
> features: 5 
> extent  : 825796.9, 831270.1, 815666.9, 816562.5  (xmin, xmax, ymin,
> ymax)
> coord. ref. : +proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0
> +y_0=0 +ellps=aust_SA +units=m +no_defs 
> variables   : 5
> names   : ID, GRIDCODE,Shape_Leng,Shape_Area, Count 
> min values  :  1,1, 94.2722341164, 6.47354525991, 1 
> max values  :  5,1, 2305.45647624, 92111.8528756, 1   
> 
>> distance.matrix<- gDistance(fnp, spgeom2= NULL, byid=T)
> Error in (function (classes, fdef, mtable)  : 
>  unable to find an inherited method for function ‘is.projected’ for
> signature ‘"list"’
> 
> 
> 
> 
> --
> View this message in context: 
> http://r.789695.n4.nabble.com/Varying-name-of-output-tables-from-looped-process-of-list-of-spdf-objects-tp4710369p4710408.html
> Sent from the R help mailing list archive at Nabble.com.
> 
> __
> R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
> 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@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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] Varying name of output tables from looped process of list of spdf objects

2015-07-27 Thread Anthoni, Peter (IMK)
Hi Cecilia,

Alternative solution

#...
afiles <- ls(pattern= "a_")

for (ifile in 1:length(afiles))
{
  fnp = mget(afiles[ifile])
#... do something with fnp

  outfile <- 
file.path("/Users/sisolarrosa/Documents/PhD/R_work/AF/IIC/conefor_inputs/", 
paste0("distances_", afiles[ifile], ".txt"));
#...
}

cheers
Peter



> On 27 Jul 2015, at 12:50, Jim Lemon  wrote:
> 
> Hi Cecilia,
> I _think_ that the error is occurring in the call to paste0. You may
> be able to get what you want like this:
> 
> paste0("distances_",deparse(substitute(fnp)),".txt")
> 
> Jim
> 
> 
> On Mon, Jul 27, 2015 at 5:06 AM, Larrosa, Cecilia
>  wrote:
>> Hi,
>> 
>> This is a repost from here 
>> (http://r.789695.n4.nabble.com/Writing-output-of-a-looped-process-with-pdfs-tt4710348.html),
>>  due to the post not being complete originally. I am running R studio on OS 
>> X Yosemite 10.10.4 (Mac). I appreciate you help very much!
>> 
>> The objective: I have 100 shapefiles that need to undergo the same process.
>> 
>> The process: I use gDistance{rgdal} to calculate the distance between all 
>> features (polygons) within each layer, and output a txt file.
>> 
>> The problem: I need the name of the output txt file to contain the name of 
>> the shapefile, but the shapefiles are read into R as 
>> SpatialPolygonsDataFrames (spdf) and I cannot find a way to use the name of 
>> the spdf objects as character in order to make it vary with each iteration.
>> 
>> My questions to you: Do you know a way to solve the problem or an 
>> alternative way to fulfil the objective? I have come to determine the 
>> problem after searching about the error message, have I interpreted 
>> correctly?
>> 
>> 
>> Here is a minimal dataset for replicability:
>> 
>>> dput(a_1)
>> new("SpatialPolygonsDataFrame"
>>, data = structure(list(ID = 1:3, GRIDCODE = c(1L, 1L, 1L), Shape_Leng = 
>> c(3349.48347556,
>> 1618.93904903, 893.268790786), Shape_Area = c(309430.38861, 90015.8325676,
>> 47507.0325775), Count = c(1L, 1L, 1L)), .Names = c("ID", "GRIDCODE",
>> "Shape_Leng", "Shape_Area", "Count"), row.names = 0:2, class = "data.frame")
>>, polygons = list(> "sp")>,
>>,
>>)
>>, plotOrder = 1:3
>>, bbox = structure(c(476685.625393809, 311791.86152084, 508519.585393809,
>> 312935.41622084), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"
>> ), c("min", "max")))
>>, proj4string = new("CRS"
>>, projargs = "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 
>> +y_0=0 +ellps=aust_SA +units=m +no_defs"
>> )
>> )
>>> dput(a_10)
>> new("SpatialPolygonsDataFrame"
>>, data = structure(list(ID = 1:5, GRIDCODE = c(1L, 1L, 1L, 1L, 1L), 
>> Shape_Leng = c(1691.7247095,
>> 2305.45647624, 1022.64650591, 1172.27848042, 94.2722341164),
>>Shape_Area = c(6.47354525991, 92111.8528756, 65.7173995386,
>>19042.7776647, 415.253663691), Count = c(1L, 1L, 1L, 1L,
>>1L)), .Names = c("ID", "GRIDCODE", "Shape_Leng", "Shape_Area",
>> "Count"), row.names = 0:4, class = "data.frame")
>>, polygons = list(> "sp")>,
>>,
>>,
>>,
>>)
>>, plotOrder = c(2L, 4L, 5L, 3L, 1L)
>>, bbox = structure(c(825796.904693809, 815666.86152084, 831270.106493809,
>> 816562.46752084), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"
>> ), c("min", "max")))
>>, proj4string = new("CRS"
>>, projargs = "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 
>> +y_0=0 +ellps=aust_SA +units=m +no_defs"
>> )
>> )
>> 
>> 
>> Here is the code that I have been using:
>> 
>> ###Load packages
>> library(rgdal)
>> library(gdistance)
>> 
>> ###Read forest shape files
>> setwd("/Users/sisolarrosa/Documents/PhD/R_work/AF/IIC/R_Quest/")
>> shps<- dir(getwd(), "*.shp")
>> shps <- gsub('.{4}$', '', shps)
>> for (shp in shps) assign(shp, readOGR(".",layer=shp))
>> 
>> ###Create list of spdf objects
>> fnps<- mget(ls(pattern= "a_"))
>> 
>> ###For each spatial layer (object in the list), calculate distance between 
>> all polygons within layer
>> for (fnp in fnps)
>> {
>>  distance.matrix<- gDistance(fnp, spgeom2= NULL, byid=T);
>>  row.names(distance.matrix) <- paste(1:nrow(distance.matrix), sep="”);   
>>  # did this because gDistance changed the IDs of the features from [1 to 
>> ...] to [0 to ...], not sure why
>>  colnames(distance.matrix)<- paste(1:ncol(distance.matrix), sep="”); 
>># same as above
>>  dists.melt <- 
>> melt(distance.matrix)[melt(upper.tri(distance.matrix))$value,];  #use only 
>> lower triangle of the distances matrix
>>  outfile <- 
>> file.path("/Users/sisolarrosa/Documents/PhD/R_work/AF/IIC/conefor_inputs/", 
>> paste0("distances_", fnp, ".txt"));
>>  write.table(dists.melt, outfile,row.names=FALSE, col.names=FALSE)
>> }
>> 
>> And this is the error message:
>> 
>> Error in as.character.default(> "SpatialPolygonsDataFrame">) :
>>  no method for coercing this S4 class to a vector
>> 
>> Thank you very much!!
>> Cecilia
>> 
>> 
>>[[alternative HTML version deleted]]
>> 
>>