Qiuhua Ma | 28 Aug 19:05 2015

criterion for the best fitting model if using GMM?


I have two big datasets:10,000 observations and 40,000.

I tried to use ML so that I can use maximized LL or AIC to decide the best
fitting model! However it took so long and most of the time R was turned
off automatically.

GMM seems like a better option. I searched the literature but cannot find
the criterion for the best fitting model. Any idea?



	[[alternative HTML version deleted]]
Ben Tupper | 28 Aug 16:33 2015

sp and latticeExtra: add colorbar for layer drawn 'under'


I have drawn on the following posting to draw a raster 'under' a set of polygons as shown in the image link
below.  The data used is available at the third link (about 650KB).

Posting: http://stackoverflow.com/questions/27062768/how-to-plot-additional-raster-with-spplot

Image: https://dl.dropboxusercontent.com/u/8433654/sp-over-plot.png

Example data: https://dl.dropboxusercontent.com/u/8433654/example.RData

Here is a brief description of each variable plotted.

R is a single layer RasterStack of sea surface temperatures
qa is a SpatialPolygonsDataFrame of fisheries management polygons
spdfbb is also a SpatialPolygonsDataFrame that contains a single polygon (created from the extent of qa)

   colorkey = FALSE,
   sp.layout = list(b = list('sp.polygons', qa, col = 'grey', first = FALSE))) + 
   latticeExtra::as.layer(spplot(R), under = TRUE)

Is it possible add a colorbar for the raster, R, as it might appear if I called...


Session info is below my signature.

(Continue reading)

Andy Bunn | 27 Aug 20:28 2015

Warning on projection (.getCRSfromGridMap4) from netcdf to brick

Hello. I'm using data out of a netcdf object from here:

The metadata has the projection info as:
#  <coordTransform name="albers_conical_equal_area"
#  <parameter name="grid_mapping_name" value="albers_conical_equal_area"/>
#  <parameter name="latitude_of_projection_origin" value="0.0 "/>
#  <parameter name="longitude_of_central_meridian" value="-120.0 "/>
#  <parameter name="standard_parallel" value="34.0 40.5 "/>
#  <parameter name="false_easting" value="0.0 "/>
#  <parameter name="false_northing" value="-4000.0 "/>
#  <parameter name="units" value="km"/>
#  <parameter name="semi_major_axis" value="6378137.0 "/>
#  <parameter name="inverse_flattening" value="298.257222101 "/>
#  </coordTransform>

I can load the data fine using raster::brick but I'm having a warning
about the projection info that I don't understand:

R> fname <- "csl_cwd_pet_test.nc"
R> petBrick <- brick(fname,varname="HST_Monthly_pet")
Warning message:
In .getCRSfromGridMap4(atts) : cannot process these parts of the CRS:
(Continue reading)

Michael Koehler | 27 Aug 14:54 2015

passing levels to predict::raster

|Dear List, I want to predict a model a rasterstack that includes 
categorical data. The outcome is all NA. However, converting the same 
rasterstack/brick to a dataframe and pedict the model gives proper 
results. What am I doing wrong  with passing the levels to 
predict::raster? ||Here is an example: 
library(raster)library(rasterVis)library(randomForest)# make a 
<-round(r)# make faktorf <-as.factor(r)is.factor(f)# get some none-sense 
levelsx <-levels(f)[[1]]x$code <-paste("A",letters[10:20])x$code2 
<-paste("B",letters[10:20])x$code3 <-letters[10:20]levels(f)<-x 
f<-deratify(f)# make a bricklevels(f)set.seed(2)# get some none-sense

fit and predict a random forest with 
itranfor<-randomForest(y~.,data=xx,ntree=100)predict(ranfor)# try to 
predict with a

gives an empty raster# convert the raster to a dataframe and predict 
works just fine! |
|#R version 3.1.2(2014-10-31)Platform:x86_64-w64-mingw32/x64 (64-bit) |I was asking that also here:
But I did not get a response.


Michael Kohler

	[[alternative HTML version deleted]]
(Continue reading)

Tim Appelhans | 27 Aug 12:49 2015

Re: How to plot multiple semi-variogram from a single dataset efficiently in R?

I see that in your lapply() call you create an object called 'plot'. 
This is bad practice as 'plot' is reserved for a function call in base R 
and may get you in trouble in certain circumstances. That is why in the 
code I gave you (and in the code that follows) this object is called 
'plt' rather than plot.

The following will create a variogram for every hour of each day, label 
them according to YYYY-mm-dd HH:MM:SS (standard POSIX datetime format) 
and save each as a file following a naming convetion of 
'plot_YYYYmmddHH.png' to avoid spaces in the file names.


seoul311 <- read.csv("Downloads/seoul1to7.csv")
seoul311 <- na.omit(seoul311)

### first we split seoul311 by time into a list rather than subsetting 
seoul311_splt <- split(seoul311, seoul311$time)

a<-as.POSIXct(names(seoul311_splt), format="%Y%m%d%H")

vars <- lapply(seq(seoul311_splt), function(i) {

(Continue reading)

Tim Peterson | 27 Aug 02:48 2015

Factorial kriging in GStat - tips for implementation?

Hi all,

I'm looking to undertake factorial colocated cokriging with GStat but 
factorial kriging does not appear to be an available option in GStat (eg 
help.search("factorial kriging") ). There was a post in 2007 requesting 
this feature (http://marc.info/?l=r-sig-geo&m=119810046000747) and a 
proposed approach in 2000 
(https://www.mail-archive.com/gstat-info <at> geog.uu.nl/msg00045.html) but 
no updates since. As far as I can tell, the required changes would be to 
'gls.c' lines 412-423, specifically the call to GCV0 for calculating C0 
(right-hand side of the kriging matrix) should be edited so that only 
the covariance for the required factor is calculated.

Anyway, before I consider adding this feature to GStat, I'd be grateful 
if anyone could share tips or thoughts. Lastly, any changes I make will 
(of course) be freely shared.

Many thanks,


Dr. Tim Peterson

The Department of Infrastructure Engineering
The University of Melbourne, 3010 Australia
T: +61 3 8344 9950, M: +61 0438 385 937

Dept. profile : 
(Continue reading)

Tim Peterson | 27 Aug 02:45 2015

GStat 'predict.c' parallel estimation

Hi all,

I'm considering editing the GStat kriging file 'predict.c' (while loop 
at line 151-165) for parallel estimation. I require parallel estimation 
for both kriging and indicator simulation and will be applying it to a 
region of ~30 million grid cells, which will be repeated for 300 time 

I appreciate that snow can be used for parallel kriging

but, because of my computation burden, I require parallel calculation 
using Xeon Phi cards - hence I think the changes should be in 
'predict.c'. Also, in looking at the GStat C-code I suspect that the 
global variables in 'glvars.c' may complicate the parallel calculations. 
So anyway, before I embark upon this task (most likely using openMP), 
I'd be grateful if anyone could share tips or thoughts (and my edits 
will be freely shared.)

Also, some may be questioning the parallel estimation of each 
simulation. Considering that I'm estimating ~30 million cells, and will 
only be using ~70 cores per simulation, I am happy to accept the 
deficiency that recently estimated cells will not be used in the 
estimation of nearby cells within a concurrent parallel cycle. However, 
the estimates do need to be accessible to latter parallel cycles.

Lastly, my backup option is to use openMP to parallelise the GSLib 
Fortran subroutines 'cokb3d.for' and 'sisim.for' (ie outside of R). This 
is mainly because the subroutines are fairly standalone and don't use 
(Continue reading)

Beaulieu, Jake | 26 Aug 19:17 2015

Stackoverflow: how-to-adjust-the-size-of-the-points-when-using-plotkml-in-r


Please follow this link to a plotKML question posted to Stackoverlfow.  I am not the original poster, but
have the same question.  The suggested solution does not work.



Jake J Beaulieu, PhD
United States Environmental Protection Agency
Office of Research and Development
National Risk Management Research Laboratory
Water Supply and Water Resources Division
Urban Waters Management Branch
26 W MLK Dr.
MS 698
Cincinnati, OH 45268
Phone: 513-569-7842

	[[alternative HTML version deleted]]
Tim Appelhans | 26 Aug 10:22 2015

Re: How to plot multiple semi-variogram from a single dataset efficiently in R?

first of all, please always reply to the list as well. As I said, other 
people might have similar problems and will thus be able to find a 
soultion in the mail archives.

If you run the code I gave you, you will see that it works.

The relevant bit is

layer(panel.text(x = 1000, y = 900, adj = c(0, 1),
                      labels = dates[panel.number()]))

where the text gets plotted at x = 1000 and y = 900. In your case, y = 
900 is beyond the limits of the y-scale.

Replacing the above with this should fix it:

layer(grid.text(x = unit(0.1, "npc"),
                     y = unit(0.9, "npc"),
                     just = c("left", "top"),
                     label = dates[panel.number()]))

This will use normalised x and y scales ranging from 0 to 1 no matter 
what the absolute units are.


On 26.08.2015 10:12, Uzzal wrote:
(Continue reading)

Tim Appelhans | 26 Aug 09:41 2015

Re: How to plot multiple semi-variogram from a single dataset efficiently in R?

it should show exactly that in the top left corner of each variogram. It 
does on my computer.

Also, I include R-si-geo again as other people might have similar problems.


On 26.08.2015 08:06, Uzzal wrote:
> I am extreamly sorry for my previous mail. Now its working.I forgot to 
> loaded "latticeExtra" package.
> The plots I got ,it all have individual names . But every plot has 7 
> variogram.is it possible to naming each variogram in each plot also 
> like "20120301","20120302",..........."201203007"?
> Uzzal
>     --- Original Message ---
>     *From : * "Tim Appelhans"<tim.appelhans <at> gmail.com>
>     *To : * "Uzzal"<uzzal <at> gist.ac.kr>
>     *Date : * 2015/08/26 Wednesday AM 1:14:25
>     *Subject : * Re: [R-sig-Geo] How to plot multiple semi-variogram
>     from a single dataset efficiently in R?
>     How's this:
>     ### in order to use latticeCombineGrid() you need to use
>     ##################
(Continue reading)

Qiuhua Ma | 25 Aug 19:26 2015

Kelejian-robinson Test Command?


I run Moran's I test and LM test for spatial autocorrelation. However I
cannot find the command for KR test. Is this command availabel in R?



	[[alternative HTML version deleted]]