Jianyun Wu | 28 Nov 06:33 2014

How to delete some polygons from 'SpatialPolygonsDataFrame'

Dear list,

I usually use 'readOGR()' from rgdal package in R to input the ESRI
So the inputted data are SpatialPolygonsDataFrame.

I am processing the New South Wales, Australia map for a spatial analysis.
There is an island far away from the mainland, which I want to delete it
from the map when I plot NSW.  Because it makes the whole plot small due to
its long distant.

Below is the area contains the polygons of this island and surrounding
islands. There are 14 " <at> Polygons" contributing to the entire area. From
some experiments, I know that it is Polygon 3-14 making the distant

Is there a simple way that I can delete those small polygons from this area
polygon and their attributes from the data

Thank you


str(nsw_lga[153, ])
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  .. <at>  data       :'data.frame':    1 obs. of  3 variables:
  .. ..$ STATE_CODE: Factor w/ 9 levels "1","2","3","4",..: 1
  .. ..$ LGA_CODE11: Factor w/ 565 levels "10050","10110",..: 153
  .. ..$ LGA_NAME11: Factor w/ 564 levels "Adelaide (C)",..: 483
  .. <at>  polygons   :List of 1
(Continue reading)

Tomislav Hengl | 27 Nov 14:35 2014

Question about derivative work - what is the license for map derived using e.g. spatial "predict" function?

Dear list,

I have a question about licensing the data that is produced by spatial 
prediction from point data. My ideas is that a map produced by using 
e.g. geostatistics from point data is a new data product and as such 
does not falls under the regulations of the original license used for 
the point data (so if the license for the point data is restrictive, the 
license for the output maps does not have to respect this). Consider for 

R> library(gstat)
R> library(sp)
R> demo(meuse, echo=FALSE)
R> m <- vgm(.59, "Sph", 874, .04)
R> x <- krige(log(zinc)~1, meuse, meuse.grid, model = m)

The produced map "x" can be considered a new data product. There is 
absolutely no way that one could reproduce the original input points 
from this map, hence it should be considered "a non-derivative work". 
Only if we would derive a map using interpolation technique that allows 
re-constructions of points (e.g. Thiessen polygons) the license would 
need to be respected.

Or am I mistaken? (I know this is a type of a question for lawyers in 
fact, but any experience / opinion you have is welcome)


"To qualify as a derivative work, the derivative must use a substantial 
(Continue reading)

Ottow, Bouke Pieter | 27 Nov 13:27 2014

RasterBrickTimeSeries in plotKML without starting GE


I am using the function of the package plotKML to plot a RasterBrickTimeSeries in GoogleEarth. I have it
working with the plotKML() function which also immediately opens GE. I am using the function in a model I am
writing and therefor I just want it to write the kml without opening it. Is there such a function as well?
Like the kml_layer.xxx() or kml() functions? It seems that it is missing for RasterBrickTimeSeries?

Kind regards,
Bouke Pieter Ottow
Manuel Ribeiro | 27 Nov 13:24 2014

Cross-variogram and respective envelopes

Hello, i am trying to calculate (and plot) cross-variograms and their
respective envelopes (bootstrap confidence intervals of empirical

I used gstat and geoR packages, but up to now I did not succeed. From what I
read (and tried) it seems these are not suited for this combination of

I have no experience in programming functions, so if anyone knows any
available function(s) to do these tasks, please let me know.

Thanks for your help, Manuel

Manuel C Ribeiro

Cerena - Centro de Recursos Naturais e Ambiente
Instituto Superior Técnico
Av. Rovisco Pais  1049-001 Lisboa Portugal
Tel: (+351) 21 841 74 41
Fax: (+351) 21 841 73 89

	[[alternative HTML version deleted]]

R-sig-Geo mailing list
R-sig-Geo <at> r-project.org
(Continue reading)

Tomislav Hengl | 25 Nov 23:15 2014

Re: questions on spatial interpolation exercise(NL)

On 25-11-2014 21:41, Mengxi Yang wrote:
> attached file is monv_201301.Please have a look.

Janurary has 31 days hence you need to adjust the number of rows in that 

 > tmp <- readLines(paste(getwd(), "/", "monv_201301.txt", sep=""))
 > tmp <- gsub(pattern=' . ', replacement=' 0.0 ', tmp, fixed = TRUE)
 > tmp <- gsub(pattern='[[:space:]]{2,7}', replacement=',', tmp)
 > tmp[8]
[1] "# 
 > write.table(t(matrix(unlist(strsplit(tmp, ",")[-c(1:8)]), nrow=42)), 
"tmp.txt", col.names=FALSE, row.names=FALSE)
 > monv_201301 <- read.table("tmp.txt", col.names=c("empty1", "NR", 
paste("D", 1:31, sep=""), "I", "NORM", "II", "NORM", "III", "NORM", 
"MND", "NORM", "Locatie"))
 > monv_201301$empty1 <- NULL
 > monv_201301$Locatie <- as.factor(gsub(pattern='[[:space:]]+$', 
replacement='', monv_201301$Locatie))
 > str(monv_201301)
'data.frame':   325 obs. of  41 variables:
  $ NR     : int  760 458 680 572 89 664 678 560 910 441 ...
  $ D1     : num  7.2 17.2 10 9.6 12 11.9 10.2 14.3 13.6 15.5 ...
  $ D2     : num  0.4 2.8 4.7 2.7 3.2 4.5 3 2.8 0.8 3.2 ...
  $ D3     : num  0.9 0.9 1.7 0.5 0.5 1.4 1 0.5 0.9 0.2 ...

(Continue reading)

Masha Pautrel | 25 Nov 06:09 2014

problem with NY8 <- readOGR(".", "NY8_utm18")

Dear all,

I try to use NY8 data set, but I can't open the file.

Here the errors message

 > library(sp)
 > library(rgdal)
 > NY8 <- readOGR(".", "NY8_utm18")
Erreur dans ogrInfo(dsn = dsn, layer = layer, encoding = encoding, 
use_iconv = use_iconv,  :
   Cannot open file
 > cities <- readOGR(".", "NY8cities")
Erreur dans ogrInfo(dsn = dsn, layer = layer, encoding = encoding, 
use_iconv = use_iconv,  :
   Cannot open file
 > TCE <- readOGR(".", "TCE")
Erreur dans ogrInfo(dsn = dsn, layer = layer, encoding = encoding, 
use_iconv = use_iconv,  :
   Cannot open file

Thank you for yours help

Masha Pautrel
Moshood Agba Bakare | 24 Nov 18:34 2014

Variable search radius for interpolation

Hi All,
Could anyone guide me on how to define variable search radius in GSTAT
while interpolating into a grid cell? The GSTAT documentation used MAXDIST
to define fixed search radius. (2) Is there a function in GSTAT to
determine the point density while interpolating into different grids size?

Thank you

	[[alternative HTML version deleted]]
Michael Krabbe Borregaard | 24 Nov 14:56 2014

Re: shifting points in plot window

This was very helpful. I have a package that depends on raster, and have included a note on this in the documentation.

From: Michael Sumner:

This is a known problem with no known fix. I've been trying to isolate it but without success.

If anyone can see how to fix it we can help get it into raster - I believe that it is somewhere in the plot
method/s for the raster package. It occurs for any vector data overlay, and for either useRaster (for
image()) set TRUE or FALSE.

Cheers, Mike

On Mon, 17 Nov 2014 21:46 mkborregaard <mkborregaard <at> snm.ku.dk<mailto:mkborregaard <at> snm.ku.dk>> wrote:
When I plot a raster image (package raster) and then plot points on top (type
SpatialPoints, package sp) everything aligns well. But if then resize the
window, the alignment disappears, and the points move around. Can I avoid
this behavior? The points and the raster are both defined as lat/long data.
proj4string is "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs"
for both.

I am sorry for posting a very simple question that I should be able to find
the resolution of, but I have not been able to find this.

View this message in context: http://r-sig-geo.2731867.n2.nabble.com/shifting-points-in-plot-window-tp7587439.html
Sent from the R-sig-geo mailing list archive at Nabble.com.

R-sig-Geo mailing list
(Continue reading)

Mark Wynter | 24 Nov 05:19 2014

R crashes performing gstat IDW with omax argument

I'm trying to use gstat's IDW function with an omax argument
Addition of the omax argument causes R to crash - works fine without it.

*** Error in `/usr/lib/R/bin/exec/R': corrupted double-linked list:
0x00000000038c8540 ***
I'm running R v 3.1 on Ubuntu 12.04

Can anyone else replicate this error on R 3.1 with the following script?


#make some random xyz point data
x = round(runif(100), 2)
y = round(runif(100), 2)
z = data.frame(rainfall=round(runif(100), 2)*1000)

#convert to spatial points data frame
xy = cbind(x, y)
xy.sp = SpatialPoints(xy)
xy.spdf = SpatialPointsDataFrame(xy.sp, z)

#prepare grid surface
x.range <- range(xy.spdf <at> coords[,1])
y.range <- range(xy.spdf <at> coords[,2])
grd <- expand.grid(x=seq(from=x.range[1]-0.5, to=x.range[2]+0.5, by=0.01),
(Continue reading)

Mengxi Yang | 23 Nov 23:51 2014

questions on spatial interpolation exercise(NL)

Dear list,
I have a question on write.table. I found an nice exercise on spatial
interpolation from

In this exercise, the date is 200906. but when I change to 201310. there is
a wrong message like:

> tmp <- readLines(paste(getwd(), "/", "monv_201310.txt", sep=""))
> tmp <- gsub(pattern=' . ', replacement=' 0.0 ', tmp, fixed = TRUE)
> tmp <- gsub(pattern='[[:space:]]{2,7}', replacement=',', tmp)
> write.table(t(matrix(unlist(strsplit(tmp, ",")[-c(1:8)]), nrow=41)),
"tmp.txt", col.names=FALSE, row.names=FALSE)
Warning message:
In matrix(unlist(strsplit(tmp, ",")[-c(1:8)]), nrow = 41) :
  data length [13625] is not a sub-multiple or multiple of the number of
rows [41]

I checked several month, some success but some not.
Does anyone know what is the wrong with it? should I change something?

Thank you very much.

Mengxi Yang

	[[alternative HTML version deleted]]
Zilefac Elvis | 23 Nov 18:13 2014

subplots of maps in R using hydroTSM package

Please I need help with maps in R using the hydroTSM package.
I am doing some map work using this wonderful package but encountered a problem using "par".
My intention is to generate 4 plots and show them on same window as follows:

if (dev.cur() == 1) x11(width=14,height=8)
par(mfrow=c(4,4),mai = c(0.8, 0.6, 0.1, 0),oma=c(0,1.4,6,0.1))


## Loading the monthly time series of precipitation within the Ebro River basin.

## Loading the gis data

## Loading the shapefile with the subcatchments

## Projection for the Subcatchments file
# European Datum 50, Zone 30N
p4s <- CRS("+proj=utm +zone=30 +ellps=intl +units=m +no_defs")

## Selecting the first day of 'EbroPPtsMonthly' for all the stations.
# The first column of 'EbroPPtsMonthly' has the dates
x.ts <- as.numeric(EbroPPtsMonthly[1, 2:ncol(EbroPPtsMonthly)])

## Setting the name of the gauging stations
(Continue reading)