Facundo Muñoz | 19 Dec 16:53 2014

kriging with non-euclidean distances

This is a recurring question from time to time.

I finally packaged [1] some old code (2009) that modifies certain
functions in geoR in order to use custom distance matrices in the
construction of empirical variograms (variog), the fit of variogram
models (likfit) and the kriging (krige.conv).

I used it at the time with “cost-based” distances, which need to be
computed elsewhere.
I used |GRASS GIS|, but nowadays there is |gdistance| which can complete
a full-R solution.
I also included some sample data and examples of usage.

I hope it is more easily accesible now.

Best regards

[1] https://github.com/famuvie/geoRcb


	[[alternative HTML version deleted]]

R-sig-Geo mailing list
R-sig-Geo <at> r-project.org

Raster projection alignment issues

Hi all,

I posted this over at gis.stackexchange.com, but thought this user group may be interested.

I downloaded  raster data from  http://gdo-dcp.ucllnl.org/downscaled_cmip_projections/dcpInterface.html#Welcomein

The BCSD-CMIP5 data is found under the tab "Projections: Complete Archives".
The files are in netCDF format, but am running into some issues. The original raster data (let's call it
rast1) when imported has the extent

xmin 235.25
xmax 293
ymin 25.125
ymax 52.875

and a crs of +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

The other objects (both sp and raster) I am working with contain negative values for xmin and xmax

xmin        : -122.8021

xmax        : -75.38249

ymin        : 31.78881

ymax        : 47.15952

When I attempt to reproject the raster data into Albers Equal Area (+proj=aea +lat_1=20 +lat_2=60
+lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83

(Continue reading)

Eko Susilo | 19 Dec 09:02 2014

How do I read and plot .nc data in R?

Could somebody help me to figure out my data

I have .nc data around Indonesia water (attached below)

I use R to open and plot it. This is my syntax

nc <- open.ncdf(choose.files('D:/R/yfin_pelikan/yfin_input', caption = 
'Select netcdf Files .....'))
print (nc)

# Read the variable

adl <- nc$var[[1]]
varsize <- adl$varsize
ndims <- adl$ndims
nt <- varsize[ndims]

start <- rep(1,ndims) # begin with start=(1,1,1,...,1)
start[ndims] <- j # change to start=(1,1,1,...,i) to read timestep i
count <- varsize # begin w/count=(nx,ny,nz,...,nt), reads entire var
count[ndims] <- 1 # change to count=(nx,ny,nz,...,1) to read 1 tstep

data <- get.var.ncdf(nc, adl, start=start, count=count )

# plot the map

jet.colors <- colorRampPalette(c('#00007F', 'blue', '#007FFF', 
'cyan','#7FFF7F', 'yellow', '#FF7F00', 'red', '#7F0000'))

plot.adl <- levelplot(data,
(Continue reading)

Rutherford Platt | 18 Dec 16:44 2014

gdal19.dll is missing from your computer


I'm a relatively new user, so bear with me.

I'm using a Windows 7 computer and R 3.1.   I wrote some scripts back in August that worked perfectly.  Now, when
the script reaches 'require(rgdal)' I get several error messages in a row that say: "The program can't
start because gdal19.dll is missing from your computer.  Try reinstalling the program to fix this problem"

I then get the following output:
Loading required package: rgdal
Loading required package: sp
rgdal: version: 0.8-16, (SVN revision 498)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
Path to GDAL shared files: C:/Users/rplatt.GETTYSBURG/Documents/R/win-library/3.1/rgdal/gdal
GDAL does not use iconv for recoding strings.
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: C:/Users/rplatt.GETTYSBURG/Documents/R/win-library/3.1/rgdal/proj

It seems like rgdal may have successfully installed, so why do I get all these error messages?  And how can I
get rid of them?


	[[alternative HTML version deleted]]
Michael Treglia | 17 Dec 00:02 2014

Extract with Large Rasters

Hi All,

I'm working with some pretty huge rasters and trying to extract data to
SpatialPointsDataFrames, and running into some issues. I'm using the raster
package for this. (sessionInfo at the end of this e-mail)

Here are the specs of my raster:

class       : RasterLayer
dimensions  : 59901, 49494, 2964740094  (nrow, ncol, ncell)
resolution  : 9.332575, 9.332575  (x, y)
extent      : 222855.6, 684762, 3762092, 4321122  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0
data source :
names       : n35w092_n39_w096_UTM_cubic
values      : 0, 837.6777  (min, max)

When I run the following line, with a vector points layer in the same CRS
as the raster, I get the error below.
>points <at> data <- data.frame(points <at> data, extract(elev, points)

Error in .readCellsGDAL(x, uniquecells, layers) :
  NAs are not allowed in subscripted assignments
In addition: Warning message:
In .readCells(x, cells, 1) : NAs introduced by coercion

> traceback()
8: .readCellsGDAL(x, uniquecells, layers)
(Continue reading)

Philippe Liege | 16 Dec 16:33 2014

gIntersection fails


I'm looking to clip the Europe map after reading this post:


The procedure in running fine when using the shape file NUTS_RG_60M_2010.shp


s <- readOGR('NUTS_2010_60M_SH/data', layer="NUTS_RG_60M_2010")

s1 <- s[s$STAT_LEVL_==2 & substring(s$NUTS_ID,1,2) %in%

# plot(s1)

clip.extent <- as(extent(-10, 19, 36, 51), "SpatialPolygons")

# clip.extent <- as(extent(-10, 40, 60, 72), "SpatialPolygons")

proj4string(clip.extent) <- CRS(proj4string(s))

gI <- gIntersects( s1 , clip.extent , byid = TRUE )

out <- lapply( which(gI) , function(x){ gIntersection( s1[x,] , clip.extent)
} )

# But let's look at what is returned
(Continue reading)

Brooks E.G.E. | 16 Dec 16:07 2014

Extract from raster produces NA when polygon has holes

Dear all

I am using extract to 'extract' the values from raster to polygons, which is working fine until I noticed a
couple of values are coming up as NA. On closer inspection these polygons have holes; it is a lake layer so
these are representing islands. I don't know why expressly the extract function doesn't like holes - is
there any way round this, or do I need to manually remove these? (I don't want to remove them as ideally I
don't want to extract the data which belongs to the islands rather than the lakes).

I am thinking this may be obvious for people in the know, but if it needs clarifying, or dummy data, apologies
and please let me know.

Many thanks

	[[alternative HTML version deleted]]
Zilefac Elvis via R-sig-Geo | 14 Dec 00:02 2014

Position of colorbar and color vs label breaks R raster

Dear All,
I have a raster layer (lr) with the following attributes:


class       : RasterLayer 

dimensions  : 137, 229, 31373  (nrow, ncol, ncell) 
resolution  : 0.136, 0.0808  (x, y) 
extent      : -120.0072, -88.8632, 48.95328, 60.02288  (xmin, xmax, ymin, ymax) 
coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
data source : in memory 
names       : myfile100 
values      : 106.6673, 172.4463  (min, max) 

I write the raster to a directory like this:
writeRaster(r, file="myfile100.tif", format="GTiff", overwrite=TRUE) 

lr = raster("myfile100.tif")
Long=c(-120.0072, -88.8632) 
Lat=c(48.95328, 60.02288)
brks <- seq(minValue(lr), maxValue(lr), by=10) 
nb <- length(brks)-1 
colors<-c("darkred", "red3","orange","gold1", "yellow", "lightskyblue","steelblue3",

(Continue reading)

Agustin Lobo | 10 Dec 13:39 2014

Problem with plotting overlaid maps

I have a problem plotting a raster layer with its legend overlaid on
a GM map. I can circunvent the problem (see below), but then
I get another problem: the resulting map has a geometry that is dependent
upon the geometry of the window. Is there any solution?

This is what I do:


colores <- rev(heat_hcl(9,c=c(80,30), l=c(30,90),
power=c(1/5,1.3),alpha = 0.75))
system("wget https://dl.dropboxusercontent.com/u/3180464/hcL1G2.tif",intern=TRUE)
g2 <- raster("hcL1G2.tif")

migmap <- gmap(x = g1, type = "hybrid", zoom = 5)
g2GM <- projectRaster(from = g2, crs = CRS("+init=epsg:3857"))
g2GM[g2GM<0] <- 0

#' Display without legend: no problems but no legend...
plot(log10(g2GM),zlim=c(1,5),col=colores, breaks=seq(1,5,by=0.5),
     legend=FALSE, add=TRUE)
#' The map is not distorted (shape independent of window geometry)

#' with legend we must do it this way
plot(log10(g2GM),zlim=c(1,5),col=colores, breaks=seq(1,5,by=0.5),
(Continue reading)

Simone Mordue (PGR | 10 Dec 12:43 2014

FW: subscript out of bounds, getvertices adehabitathr

Hope someone can help, I have 10 animals and have generated random points based on recorded polygons for
each animal. I have run this 100 times per animal. My aim is to generate a mean utilized distribution for
each animal based on all 100 runs. so far I have used this code:

xybat <- subset(bat.master, bat.master$id =="Y2889a",select=x:loopno )

# change to spatial points
xy <- xybat[1:2 ] # first two rows save as coords
df <- xybat[-1:-3 ]# remove unneded columns for ud

SPDF <- SpatialPointsDataFrame(coords=xy, data=df ) # combine df and

udHR <- kernelUD(SPDF, h = 200, grid=habitat)

## I would proceed using the raster packages
# stack the UD's as a raster
ud <- stack(lapply(udHR, raster))

## You can now check the first one

## take the mean
plot(udm <- mean(ud))

## now convert it back to a estUD, this is a bit of a hack and not the nicest way to do it

Y2889a<- udHR[[1]]
Y2889a <at> grid <- as(udm, "GridTopology")

I do this first bit of code for each animal, here is a subset of just two of them. I can visulise them using
(Continue reading)

Vasya Pupkin | 10 Dec 07:03 2014

bootstrapping linear model and residuals kriging by RSAGA

Dear ladies and gentlemen, I don’t really know it this is a suitable platform for asking such questions,
but stackoverflow could not help me, so I am trying my luck here:).
What I am doing: I am trying to make a soil erodibility map. I have soil erodibility measurements at points
and I do a linear regression analysis to predict soil erodibility from predictors (morphological and
remotely sensed). Then I apply the regression equation to the predictor rasters and get a raster of soil
erodibility. Then I calculate residuals and extra/inter-polate with ordinary kriging (global) using
my own formula discribing the variogramm. Then I summ up the map from regression equation and the
interpolated rsiduals to create the final map, thus mimiking regression kriging as you have suggested.
The regression was done in R.
Now I want to validate the final map with bootstrap as it is implimented in R "boot" package. In the madual it
says you can wrap any complicated model in the "boot" sequence. So I wanted to write this whole
regression+residual_kriging model in the boostrap in R. After googling I fould this RSAGA package that
provides a direct access to SAGA's functionality. And there is this "rsaga.geoprocessor" module which
works like CMD for SAGA, and will not return the results back to R environment for further calculations but
will save them to a file, like SAGA's CMD normally does. And there is this "pick.from.points" module which
has sort of deeper integration with SAGA's modules (not all though) and is what I think I need but I can't
figure out how I can use my own variogram formula and how I can return the results back to R environment for
further calculations. So in other words I need to do a regression and then predict the values and
extrapolate residuals for excluded points, so the RSAGA module should return the values of extrapolated
residuals back to “boot” calculation.
What I want to do is something like this (R script):
       lm( resp~predict, ...)+
       pick.from.points( lm(resp~predict)$residuals, method = "krige", model = "a+b^c",...))
Do you have any idea how it could be done? If not, could you please advise on any other ways or anybody who could
know? The approach – regression + residuals kriging is critical, but not the software, so if there are
any other packages for doing this, please advice. 

Best wishes,
(Continue reading)