From f616d2d3a9a9012d62bdbda64dd37c8b2cec9b9f Mon Sep 17 00:00:00 2001
From: Philippe Verley <philippe.verley@ird.fr>
Date: Wed, 7 Sep 2016 23:24:34 +0000
Subject: [PATCH] This commit adds a new option to the ads::mimetic function:
the user may provide a raster in the arguments for defining a mask.
mimetic(..., mask=NULL) mask must either inherit from RasterStack or
RasterLayer at the moment (we might change this in future commit). Every non
NA values are considered as allowed values for positioning new points. It
only works for rectangular windows without triangle at the moment (we will
extend the option to other window type in future commits).
In the C sources, I updated the functions mimetic_rect(), exchange_point_rect() and created a new function s_alea_rect_mask based on s_alea_rect. The main differences is that instead of
generating new points anywhere within the sampling window, it only generates points within the raster limits inside the sampling window.
The precision parameter might not be implemented properly for the masked point process, it should be double checked.
The updated functions take the additional arguments nmask, xmask, ymask, dx and dy:
nmask: number of non-NA points of the raster
xmask: x-coordinates of the non-NA points in the raster
ymask: y-ccordinates of the non-NA points in the raster
dx: dimension of one raster cell along the x-axis
dy: dimension of one raster cell along the y-axis
The ads::mimetic function does not check whether the spatial point pattern and the raster are consistent in terms of CRS. It assumes at the moment that they are expressed in the same CRS. The function gets rid of any point from the raster that are outside the rectangular sampling window.
---
R/mimetic.R | 24 +++++++++++-
src/Zlibs.c | 103 +++++++++++++++++++++++++++++++++++++++++++++++-----
src/Zlibs.h | 5 ++-
3 files changed, 120 insertions(+), 12 deletions(-)
diff --git a/R/mimetic.R b/R/mimetic.R
index c428df8..ab0c004 100755
--- a/R/mimetic.R
+++ b/R/mimetic.R
@@ -2,7 +2,7 @@
#RP 11/06/2013
###################################################
-mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50) {
+mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50,mask=NULL) {
# checking for input parameters
stopifnot(inherits(x,"fads")||inherits(x,"spp"))
if(inherits(x,"fads")) {
@@ -50,13 +50,35 @@ mimetic<-function(x,upto=NULL,by=NULL,prec=NULL,nsimax=3000,conv=50) {
PACKAGE="ads")
}
else {
+ if (!missing(mask) && !is.null(mask)) {
+ stopifnot(inherits(mask, "RasterStack")||inherits(mask, "RasterLayer"))
+ if (inherits(mask, "RasterStack"))
+ mask <- mask[[1]]
+ mask.coord <- coordinates(mask)[!is.na(values(mask)),]
+ mask.n <- nrow(mask.coord)
+ mask.x <- mask.coord[,1]
+ mask.y <- mask.coord[,2]
+ mask.dx <- abs(mask@extent@xmax-mask@extent@xmin)/mask@ncols
+ mask.dy <- abs(mask@extent@ymax-mask@extent@ymin)/mask@nrows
+ mask.spp <- spp(x=mask.coord[,1],y=mask.coord[,2],window=p$window)
+ res<-.C("mimetic_rect",
+ as.integer(p$n),as.double(p$x),as.double(p$y),as.double(surface),
+ as.double(xmin),as.double(xmax),as.double(ymin),as.double(ymax),
+ as.double(mask.spp$n),as.double(mask.spp$x), as.double(mask.spp$y),as.double(mask.dx),as.double(mask.dy),
+ as.double(prec),as.integer(tmax),as.double(by),
+ as.double(lobs),as.integer(nsimax),as.integer(conv),cost=double(nsimax),
+ g=double(tmax),k=double(tmax),xx=double(p$n),yy=double(p$n),mess=as.integer(1),
+ PACKAGE="ads")
+ } else {
res<-.C("mimetic_rect",
as.integer(p$n),as.double(p$x),as.double(p$y),as.double(surface),
as.double(xmin),as.double(xmax),as.double(ymin),as.double(ymax),
+ as.integer(0),as.double(0), as.double(0),as.double(0),as.double(0),
as.double(prec),as.integer(tmax),as.double(by),
as.double(lobs),as.integer(nsimax),as.integer(conv),cost=double(nsimax),
g=double(tmax),k=double(tmax),xx=double(p$n),yy=double(p$n),mess=as.integer(1),
PACKAGE="ads")
+ }
}
}
else if("circle"%in%p$window$type) {
diff --git a/src/Zlibs.c b/src/Zlibs.c
index 362f14b..e089f91 100755
--- a/src/Zlibs.c
+++ b/src/Zlibs.c
@@ -2572,7 +2572,7 @@ int intertype_rect_ic(int *point_nb1, double *x1, double *y1, int *point_nb2, do
r=0;
while (erreur!=0)
{
- erreur=mimetic_rect(point_nb1, x1, y1, surface, xmi, xma, ymi, yma, prec, t2, dt, lt, nsimax, conv, cost, gt, kt, x, y, 0);
+ erreur=mimetic_rect(point_nb1, x1, y1, surface, xmi, xma, ymi, yma, 0, NULL, NULL, 0, 0, prec, t2, dt, lt, nsimax, conv, cost, gt, kt, x, y, 0);
r=r+erreur;
if (r==*rep)
{
@@ -3688,6 +3688,30 @@ void s_alea_rect(int point_nb, double x[], double y[],
PutRNGstate();
}
+void s_alea_rect_mask(int point_nb, double x[], double y[],
+ double nmask, double *xmask, double *ymask, double dx, double dy,
+ double xmi, double xma, double ymi, double yma,
+ double p)
+ {
+ int i, out;
+ size_t ii;
+ double xr, yr;
+ GetRNGstate();
+ for (i=0; i<point_nb; i++)
+ {
+ ii=nmask*unif_rand();
+ out=1;
+ while (out)
+ {
+ x[i]=*(xmask+ii)+((2*unif_rand()-1)*(dx/p))*p;
+ y[i]=*(ymask+ii)+((2*unif_rand()-1)*(dy/p))*p;
+ out=(x[i]<xmi)||(x[i]>xma)||(y[i]<ymi)||(y[i]>yma);
+ }
+ //Rprintf("s_alea_rect_mask, ii %d, xmask %f, ymask %f, x %f, y %f\n", ii, *(xmask+ii), *(ymask+ii), x[i], y[i]);
+ }
+ PutRNGstate();
+ }
+
/*pour un zone circulaire*/
void s_alea_disq(int point_nb, double *x, double *y, double x0, double y0, double r0, double p)
{
@@ -7741,11 +7765,14 @@ int randomdist(int *vec, int nb_type, double *mat, double *matp)
/*Mimetic point process as in Goreaud et al. 2004 */
/******************************************************/
-int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *xmi, double *xma, double *ymi, double *yma,
+int mimetic_rect(int *point_nb, double *x, double *y, double *surface,
+ double *xmi, double *xma, double *ymi, double *yma,
+ double *nmask, double *xmask, double *ymask, double *dx, double *dy,
double *prec, int *t2, double *dt, double *lobs, int *nsimax, int *conv, double *cost,
double *g, double *k, double *xx, double *yy, int *mess)
{
int i, compteur_c=0, r=0, erreur=0;
+ size_t ii;
int compteur=0;
double *l;
double cout, cout_c;
@@ -7756,9 +7783,30 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
double offsetx=((*xmi)<0)?(*xmi):0.;
double offsety=((*ymi)<0)?(*ymi):0.;
decalRect(*point_nb, x, y, xmi, xma, ymi, yma);
+ if ((*nmask)>0)
+ {
+ if (offsetx<0.)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ double avant= *(xmask+ii);
+ *(xmask+ii)-=offsetx;
+ }
+ }
+ if (offsety<0.)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ *(ymask+ii)-=offsety;
+ }
+ }
+ }
//creation of a initial point pattern and cost
- s_alea_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *prec);
+ if ((*nmask)>0)
+ s_alea_rect_mask(*point_nb, x, y, *nmask, xmask, ymask, *dx, *dy, *xmi, *xma, *ymi, *yma, *prec);
+ else
+ s_alea_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *prec);
erreur=ripley_rect(point_nb, x, y, xmi, xma, ymi, yma, t2, dt, g, k);
if (erreur!=0) return -1;
cout=0;
@@ -7774,7 +7822,7 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
cout_c=0;
while (compteur<*nsimax)
{
- cout_c=echange_point_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, intensity, *prec, cout, lobs, t2, dt, g, k);
+ cout_c=echange_point_rect(*point_nb, x, y, *xmi, *xma, *ymi, *yma, *nmask, xmask, ymask, *dx, *dy, intensity, *prec, cout, lobs, t2, dt, g, k);
if (cout==cout_c)
compteur_c++;
else
@@ -7810,17 +7858,40 @@ int mimetic_rect(int *point_nb, double *x, double *y, double *surface, double *x
*xma+=offsetx;
*ymi+=offsety;
*yma+=offsety;
+ // Shift mask into original location
+ if ((*nmask)>0)
+ {
+ if (offsetx<0)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ xmask[ii]=xmask[ii]+offsetx;
+ }
+ }
+ if (offsety<0)
+ {
+ for (ii=0; ii<(*nmask); ii++)
+ {
+ ymask[ii]=ymask[ii]+offsety;
+ }
+ }
+ }
free(l);
return r;
}
-double echange_point_rect(int point_nb, double *x, double *y, double xmi, double xma, double ymi, double yma, double intensity, double p, double cout, double * lobs, int *t2, double *dt, double *g, double *k)
+double echange_point_rect(int point_nb, double *x, double *y,
+ double xmi, double xma, double ymi, double yma,
+ double nmask, double *xmask, double *ymask, double dx, double dy,
+ double intensity, double p, double cout, double * lobs,
+ int *t2, double *dt, double *g, double *k)
{
double xr, yr, xcent[4], ycent[4], n_cout[4], *l, xprec, yprec;
- int erreur, max, i, j, num;
+ int erreur, max, i, j, num, out;
+ size_t ii;
vecalloc(&l, *t2);
GetRNGstate();
- num=unif_rand()* (point_nb); // numero de l'arbre que l'on retire
+ num=unif_rand()*(point_nb); // numero de l'arbre que l'on retire
xprec=x[num];
yprec=y[num];
xr=xma-xmi;
@@ -7835,8 +7906,22 @@ double echange_point_rect(int point_nb, double *x, double *y, double xmi, double
for (j=0; j<4; j++)
{
- xcent[j]=xmi+(unif_rand()*(xr/p))*p; // coordonn�es (x,y) du point tir�
- ycent[j]=ymi+(unif_rand()*(yr/p))*p;
+ if (nmask>0)
+ {
+ ii=nmask*unif_rand();
+ out=1;
+ while (out)
+ {
+ xcent[j]=xmask[ii]+((2*unif_rand()-1)*(dx/p))*p;
+ ycent[j]=ymask[ii]+((2*unif_rand()-1)*(dy/p))*p;
+ out=(xcent[j]<xmi)||(xcent[j]>xma)||(ycent[j]<ymi)||(ycent[j]>yma);
+ }
+ }
+ else
+ {
+ xcent[j]=xmi+(unif_rand()*(xr/p))*p; // coordonn�es (x,y) du point tir�
+ ycent[j]=ymi+(unif_rand()*(yr/p))*p;
+ }
x[num]=xcent[j];
y[num]=ycent[j];
erreur=ripley_rect(&point_nb, x, y, &xmi, &xma, &ymi, &yma, t2, dt, g, k);
diff --git a/src/Zlibs.h b/src/Zlibs.h
index 9130b52..7b16155 100755
--- a/src/Zlibs.h
+++ b/src/Zlibs.h
@@ -78,6 +78,7 @@ int intertypelocal_tr_disq(int *,double *,double *,int *,double *,double *,doubl
double *,double *,double *,double *,double *,double *,int *,double *,double *,double *);
void s_alea_rect(int,double[],double[],double,double,double,double,double);
+void s_alea_rect_mask(int,double[],double[],double,double *,double *,double,double,double,double,double,double,double);
void s_alea_disq(int,double *,double *,double,double,double,double);
void s_alea_tr_rect(int,double *,double *,double,double,double,double,int,double *,double *,double *,double *,
double *,double *,double);
@@ -155,13 +156,13 @@ int rao_tr_disq(int *,double *,double *,double *,double *,double *,int *,double
int rao_tr_disq_ic(int *,double *,double *,double *,double *,double *,int *,double *,double *,double *,double *,double *,double *,int *,double *,int *,int *,double *,
int *,int *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,int *);
-int mimetic_rect(int *,double *,double *, double *,double *,double *,double *,double *,double *, int *,double *,double *,int *,int *,double *,double *,double *,double *,double *,int *);
+int mimetic_rect(int *,double *,double *, double *,double *,double *,double *,double *,double *,double *,double *,double *,double *,double *, int *,double *,double *,int *,int *,double *,double *,double *,double *,double *,int *);
int mimetic_disq(int *,double *,double *,double *,double *,double *,double *,double *, int *, double *, double *, int *, int *, double *,double *,double *,double *,double *,int *);
int mimetic_tr_rect(int *,double *,double *, double *,double *,double *,double *,double *,int *, double *, double *, double *, double *, double *, double *,
double *, int *, double *, double *, int *, int *, double *,double *, double *,double *,double *,int *);
int mimetic_tr_disq(int *,double *,double *, double *,double *,double *,double *,int *, double *,double *,double *,double *,double *,double *cy,
double *, int *, double *, double *, int *, int *, double *,double *,double *,double *,double *,int *);
-double echange_point_rect(int,double *,double *,double,double,double,double,double,double,double,double *,int *,double *,double *,double *);
+double echange_point_rect(int,double *,double *,double,double,double,double,double,double *,double *,double,double,double,double,double,double *,int *,double *,double *,double *);
double echange_point_disq(int,double *,double *,double,double,double,double,double,double,double *,int *,double *,double *,double *);
double echange_point_tr_rect(int,double *,double *,double,double,double,double,int *,double *,double *,double *,double *,double *,double *,
double,double,double,double *,int *,double *,double *,double *);
--
GitLab