Nick Matzke | 27 Apr 00:37 2015

Plot phylogeny with "increasing" or "decreasing" node ages?

Hi all,

In FigTree, there is an option (Left menu->Trees->Order nodes->increasing,
or decreasing) to plot trees and order them such that the higher nodes/tips
in the tree plot at the top or bottom of the plot.

Does anything like this exist for plotting trees in R?  Or do I have to
hunt-n-peck and manually rotate nodes 1-by-1 until it looks right?

(yes, I could save the tree out from FigTree and use that, but I want to
automate this plotting function without a FigTree step)


	[[alternative HTML version deleted]]
drflying | 26 Apr 07:36 2015

自动回复:Re: Count Pixel with Specific of Values Appearing Sequentially in Rasterstack

R-sig-Geo mailing list
R-sig-Geo <at>

Count Pixel with Specific of Values Appearing Sequentially in Rasterstack

Dear Colleagues,
Any idea on how to count the number of times a particular number (e.g. 1) appears consecutively within each
pixel of a raster stack/brick? I found a similar question online
but was not able to view/get any answer to solve such problem.
Thanks in advance.

	[[alternative HTML version deleted]]
drflying | 24 Apr 23:25 2015

自动回复: Can plotKML produces <region> kmls to plot raster stack for large areas?

R-sig-Geo mailing list
R-sig-Geo <at>
Amit Boshale via R-sig-Geo | 24 Apr 23:24 2015

Can plotKML produces <region> kmls to plot raster stack for large areas?

Hi All, 

I want to display MODIS rasterStack (300 images, area 1200 km x 1200 km)on Google Earth using plotKML, can I
do that?
The code provided in RasterBrickTimeSeries in the tutorial worked well for small areas, but I don't know
what to do for large areas
Tuturial code (working on my data)
> require(raster)
> require (plotKML)
> data(LST)
> gridded(LST) <- ~lon+lat
> proj4string(LST) <- CRS("+proj=longlat +datum=WGS84")
> dates <- sapply(strsplit(names(LST), "LST"), function(x){x[[2]]})
> datesf <- format(as.Date(dates, "%Y_%m_%d"), "%Y-%m-%dT%H:%M:%SZ")
> # begin / end dates +/- 4 days:
> TimeSpan.begin = as.POSIXct(unclass(as.POSIXct(datesf))-4*24*60*60, origin="1970-01-01")
> TimeSpan.end = as.POSIXct(unclass(as.POSIXct(datesf))+4*24*60*60, origin="1970-01-01")
> # pick few climatic stations:
> pnts <- HRtemp08[which(HRtemp08$NAME=="Pazin")[1],]
> pnts <- rbind(pnts, HRtemp08[which(HRtemp08$NAME=="Crni Lug - NP Risnjak")[1],])
> pnts <- rbind(pnts, HRtemp08[which(HRtemp08$NAME=="Cres")[1],])
> coordinates(pnts) <- ~Lon + Lat
> proj4string(pnts) <- CRS("+proj=longlat +datum=WGS84")
> # get the dates from the file names:
> LST_ll <- brick(LST[1:5])
> LST_ll <at> title = "Time series of MODIS Land Surface Temperature (8-day mosaics) images"
> LST.ts <- new("RasterBrickTimeSeries", variable = "LST", sampled = pnts,
+     rasters = LST_ll, TimeSpan.begin = TimeSpan.begin[1:5], TimeSpan.end = TimeSpan.end[1:5])
> data(SAGA_pal)
(Continue reading)

Juta Kawalerowicz | 23 Apr 19:04 2015

Spatial overlay vs. extract and aggregation methods


I have a spatial polygon data frame (spdf1) and try to aggregate
information from it over polygons in a larger spatial polygon data
frame(spdf2). I wanted to run some simple test to see that my code is sound
(on nc.sids from maptools), but run into a problem which I don't
understand. As far as I can see, there are 2 ways to do this:

1) rasterise the smaller spdf1, extract and calculate over larger spdf2

nc.sids <- readShapeSpatial(system.file("shapes/sids.shp",
                            IDvar="FIPSNO", proj4string=CRS("+proj=longlat

r<-raster(ncol=180, nrow=180)
rp<-rasterize(nc.sids, r, 'BIR74')
plot(nc.sids, add=TRUE)

#this will take time on larger files...
v <- extract(rp, nc.sids, weights=TRUE)
res<-sapply(v, function(x) if (!is.null(x)) {sum(apply(x, 1, prod),
na.rm=TRUE) / sum(x[,2])} else NA )
cor(res, nc.sids$BIR74)
#seems to be working ok.
(Continue reading)

Luciano La Sala | 23 Apr 16:36 2015

Problem clipping SpatialPolygonDataFrame

Hello everyone,

I have a SpatialPolygonsDataFrame for Argentina which I need to crop 
using a bounding box or any other method. I tried the following but the 
new map turns out to have the same exact extent as the original one:

# Info on my Data Frame.

 > summary(arg.provinces)
Object of class SpatialPolygonsDataFrame
         min       max
x -73.57778 -53.59183
y -55.06153 -21.77855
Is projected: FALSE
proj4string :
[+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]

 > extent(arg.provinces)
class       : Extent
xmin        : -73.57778
xmax        : -53.59183
ymin        : -55.06153
ymax        : -21.77855

# Trying to clip it

out <- crop(arg.provinces, extent(-73.57778, -53.59183, -55.06153, 
(Continue reading)

Gilles Benjamin Leduc | 23 Apr 02:45 2015

adding a scale

Hi all,

I am ploting map from shapefiles (curently read with readOGR from rgdal). 
My PhD supervisor asked me to add a scale… Before going in crazy computation and the arrow function I
wonder… Is there a function that can do it automatically (or easily) 

Thanks in advance 


R-sig-Geo mailing list
R-sig-Geo <at>
Beaulieu, Jake | 22 Apr 20:14 2015

plotKML throws status 37 warning

I am having trouble running the example code provided in the plotKML help file.  The example is as follows:

coordinates(eberg) <- ~X+Y
proj4string(eberg) <- CRS("+init=epsg:31467")
## subset to 20 percent:
eberg <- eberg[runif(nrow(eberg))<.1,]
## Not run: ## bubble type plot:

Here is what happens when I run the code:

> plotKML(eberg["CLYMHT_A"])

Plotting the first variable on the list

KML file opened for writing...

Reprojecting to +proj=longlat +datum=WGS84 ...

Writing to KML...

Closing  eberg__CLYMHT_A__.kml

access to 'C:\Users\jbeaulie\DOCUME~1\PLOTKM~1\EBERG_~1.KML' denied

Warning message:

(Continue reading)

Stephen Stewart | 22 Apr 01:28 2015

MODIS Package 0.10-23 (16/04/2015)

Hi everyone,

I was testing a fresh install of the MODIS package on my laptop yesterday in
preparation for a workshop I am running later on in the week.

I came across an error in that I have previously not encountered with MODIS


MODISoptions(gdalPath='c:/OSGeo4W64/bin', localArcPath =
"C:/GIS_in_R/MODIS/",outDirPath = "C:/GIS_in_R/MODIS/DONE/") # Set the GDAL
and output paths.

start_date <- as.Date("1-1-2014","%d-%m-%Y") 

end_date <- as.Date("1-1-2014","%d-%m-%Y") 

dates <- transDate(start_date,end_date) 

runGdal("MOD11A2",extent = extent(vic), begin = dates$beginDOY, end =
dates$endDOY, outProj="+proj=longlat +datum=WGS84 +no_defs")

Error in file (description = lockfile, open = "wt") :
cannot open the connection
In addition: warning message:
In file(description = lockfile, open = "wt") :
cannot open file
'C:/GIS_in_R/MODIS/DONE/.auxilliaries/MOD11A2.005.LPDAAC.lock': No such file
or directory.
(Continue reading)

Guido Schulz | 20 Apr 14:05 2015

How to do residual & regression deletion diagnostics on spautolm/sarlm objects?

In stats there are nice functions to carry out residual & regression
deletion diagnostics on lm or glm objects, such as:


lm.object <- lm(sr ~ pop75 + dpi + ddpi, data = LifeCycleSavings)

I would like to do these (or similar) diagnostics on a spdep::spautolm
(type "SAR") object, but neither plot.lm nor influence.measures accepts
spdep::spautolm/sarlm objects.

Can you recommend a workaround?
Or are these diagnostics simply not appropriate for SAR type regressions?