Re: st: Randomize points to map polygons using -spmap- ?

[Maurizio Pisati <maurizio.pisati@unimib.it>]

Dear Andrew,
good point. Below is a modified version of -csrpp- that generates and saves the polygon's ID of each point:

. use "ADMdb.dta", clear
. generate cases = 50
. csrpp cases using "ADMcoord.dta", id(_ID) saving("casecoord.dta") replace
. spmap using "ADMcoord.dta", id(_ID) point(data("casecoord.dta") x(_X) y(_Y) shape(i)) label(data("casecoord.dta") x(_X) y(_Y) label(_ID))

Best wishes,
Maurizio


===== CODE STARTS HERE =============================================================

*! -csrpp-: Generate completely spatially random points in a polygon           
*! Version 1.1 - 17 July 2013 (based on suggestion by Andrew Lover)            
*! Version 1.0 - 19 August 2006                                                
*! Author: Maurizio Pisati                                                     
*! Department of Sociology and Social Research                                 
*! University of Milano Bicocca (Italy)                                        
*! maurizio.pisati@unimib.it                                                   




*  ----------------------------------------------------------------------------
*  1. Define program                                                           
*  ----------------------------------------------------------------------------

program csrpp
version 9.2




*  ----------------------------------------------------------------------------
*  2. Define syntax                                                            
*  ----------------------------------------------------------------------------

syntax varname(numeric) [if] [in] [iweight/] using/,   ///
       ID(varname numeric)                             ///
      [Scale(numlist max=1 >0)]                        ///
       SAVing(string)                                  ///
      [REPLACE]                                        ///
      [noDOTS]




*  ----------------------------------------------------------------------------
*  3. Check syntax                                                             
*  ----------------------------------------------------------------------------

/* Check weight */
if ("`weight'" != "") {
   tempvar W
   qui gen `W' = `exp'
   qui su `W'
   if r(min) < 0 {
      di as err "{p}The weighting variable can take only nonnegative "   ///
                "values{p_end}"
      exit 198
   }
}
else local exp = 1

/* Check option id() */
capture isid `id'
if _rc {
   di as err "{p}Variable `id' specified in option {bf:{ul:id}()} "   ///
             "does not uniquely identify the polygons{p_end}"
   exit 459
}

/* Check option saving() */
if (substr(reverse("`saving'"),1,4) != "atd.") {
   local region "`saving'.dta"
}

/* Marksample */
marksample TOUSE
markout `TOUSE' `id'
qui count if `TOUSE'
if (r(N) == 0) error 2000




*  ----------------------------------------------------------------------------
*  4. Define basic objects                                                     
*  ----------------------------------------------------------------------------

/* Preserve data */
preserve

/* Select relevant observations */
qui keep if `TOUSE'

/* Set default scale factor */
if ("`scale'" == "") local scale = 1

/* Set option dots */
if ("`dots'" != "") local dots "no"
else local dots "yes"




*  ----------------------------------------------------------------------------
*  5. Data setup                                                               
*  ----------------------------------------------------------------------------

/* Generate id variable */
cap rename `id' _ID
sort _ID

/* Generate variables stating number of points to be generated */
unab varlist : `varlist'
cap rename `varlist' NP
qui replace NP = round(NP * `scale')

/* Generate weight variable */
qui gen W = `exp' * (1 / `scale')

/* Save temporary file */
keep _ID NP W
tempfile DATA
qui save `DATA', replace




*  ----------------------------------------------------------------------------
*  6. Polygon setup                                                            
*  ----------------------------------------------------------------------------

/* Open file */
use "`using'", clear

/* Generate subpolygon id */
quietly {
   gen POLY = (_X == .)
   replace POLY = sum(POLY)
   drop if _X == .
   sort POLY, stable
   tempfile POLY
   save `POLY', replace
}

/* Compute polygon area */
quietly {
   by POLY : gen double AREA = (_X*_Y[_n+1]) - (_X[_n+1]*_Y) if _n>1 & _n<_N
   by POLY : replace AREA = sum(AREA)
   by POLY : replace AREA = abs(AREA[_N] / 2)
}

/* Compute subpolygon weigths */
collapse AREA (count) NV = AREA, by(_ID POLY)
quietly {
   by _ID : gen double CUM = sum(AREA)
   by _ID : replace AREA = AREA / CUM[_N]
   drop CUM
}




*  ----------------------------------------------------------------------------
*  7. Prepare working dataset                                                  
*  ----------------------------------------------------------------------------

/* Merge files and adjust NP variable */
quietly {
   merge _ID using `DATA'
   keep if _merge == 3
   replace NP = round(AREA * NP)
   gen ID = _n
}

/* Select relevant variables and save dataset */
qui generate _ID_ = _ID
keep POLY ID NP NV W _ID_
order POLY ID NP NV W _ID_
sort POLY
qui save `DATA', replace

/* Select relevant records from polygon file */
quietly {
   use `POLY', clear
   merge POLY using `DATA'
   keep if _merge == 3
   keep ID _X _Y
   ren ID _ID
   order _ID
   sort _ID, stable
   save `POLY', replace   
}

/* Open data  */
use `DATA', clear
keep ID NP NV W _ID_
qui gen END = sum(NV)
qui gen START = END[_n-1] + 1
qui replace START = 1 in 1
drop NV
order ID NP START END W _ID_

/* Count polygons */
qui count
local NPOLY = r(N)

/* Count points to be generated */
su NP, mean
local NPOINTS = r(sum)

/* Merge polygon file */
qui merge using `POLY'
drop _merge ID _ID

/* Set number of records */
qui count
if (`NPOINTS' > r(N)) qui set obs `NPOINTS'

/* Generate variables */
qui gen ID = .
qui gen X  = .
qui gen Y  = .
qui gen _W = .




*  ----------------------------------------------------------------------------
*  8. Generate points and save dataset                                         
*  ----------------------------------------------------------------------------

/* Generate points */
mata: csrpp(`NPOLY', "`dots'")

/* Create and select relevant variables */
keep ID X Y _W
ren ID `id'
ren X _X
ren Y _Y
lab var `id' "Polygon ID"
lab var _X "x-coordinate"
lab var _Y "y-coordinate"
lab var _W "Weight"

/* Save dataset */
qui drop if `id' == .
qui compress
save "`saving'", `replace'




*  ----------------------------------------------------------------------------
*  9. End program                                                              
*  ----------------------------------------------------------------------------

restore
end








/*  -----------------------------------------------------------------------  */
/*  Call mata                                                                */
/*  -----------------------------------------------------------------------  */

version 9.2
mata:
mata clear
mata set matastrict on




/*  -----------------------------------------------------------------------  */
/*  csrpp()                                                                  */
/*  -----------------------------------------------------------------------  */
void csrpp(real scalar npoly, string scalar dots)
{


real matrix    POLY, COOR
real scalar    i, line
real scalar    np, id, we, fr, lr, fi, li

if (dots == "yes") {
   printf("\n")
   printf("{txt}Scanning polygons ({res}%g{txt})\n", npoly)
   printf("{txt}{hline 4}{c +}{hline 3} 1 {hline 3}{c +}{hline 3} 2 ") 
   printf("{txt}{hline 3}{c +}{hline 3} 3 {hline 3}{c +}{hline 3} 4 ")
   printf("{txt}{hline 3}{c +}{hline 3} 5\n")
   displayflush()
}

fi = 1
li = 0
for (i=1; i<=npoly; i++) {
   if (dots == "yes") {
      line = mod(i,50)
      if (line != 0  &  i < npoly) {
         printf("{txt}.")
         displayflush()
      }
      if (line == 0  &  i < npoly) {
         printf("{txt}. %5.0f\n", i)
         displayflush()
      }
      if (i == npoly) {
         printf("{txt}.\n")
         printf("\n")
         displayflush()
      }
   }
   np = _st_data(i, 1)
   if (np > 0) {
      id = _st_data(i, 5)
      we = _st_data(i, 4)
      li = li + np
      fr = _st_data(i, 2)
      lr = _st_data(i, 3)
      POLY = st_data((fr,lr), (6..7))
      COOR = sp_csrpp(np, POLY)
      st_store((fi,li), 8, J(np, 1, id))
      st_store((fi,li), (9..10), COOR)
      st_store((fi,li), 11, J(np, 1, we))
      fi = li + 1
   }
}


}




/*  -----------------------------------------------------------------------  */
/*  sp_csrpp()                                                               */
/*  -> sp_pip()                                                              */
/*  -> sp_polyarea()                                                         */
/*                                                                           */
/*  Returns matrix npoints-by-2 containing the coordinates of npoints        */
/*  completely spatially random points lying within polygon POLY             */
/*  -----------------------------------------------------------------------  */
real matrix sp_csrpp(real scalar npoints, real matrix POLY)
{


/* Setup */
real matrix    COOR
real matrix    XC, YC, XY, PIP, DONE
real scalar    xmin, xmax, ymin, ymax, areap, areab, aratio
real scalar    ndone, nmiss, ngen

/* Generate working objects */
xmin = min(POLY[., 1])
xmax = max(POLY[., 1])
ymin = min(POLY[., 2])
ymax = max(POLY[., 2])
areap = sp_polyarea(POLY)
areab = abs(xmax - xmin) * abs(ymax - ymin)
aratio = areab / areap

/* Generate points */
ndone = 0
while (ndone < npoints) {
   nmiss = npoints - ndone
   ngen = round(nmiss * aratio)
   if (ngen == 1) ngen = 10
   XC = xmin :+ uniform(ngen, 1) :* (xmax - xmin)
   YC = ymin :+ uniform(ngen, 1) :* (ymax - ymin)
   XY = XC , YC
   PIP = sp_pip(XY, POLY)
   DONE = select(XY, PIP)
   if (ndone == 0) COOR = DONE
   else COOR = COOR \ DONE
   ndone = rows(COOR)
}

/* Return results */
return(COOR[(1..npoints), .])


}




/*  -----------------------------------------------------------------------  */
/*  sp_pip()                                                                 */
/*  -> sp_spip()                                                             */
/*                                                                           */
/*  Returns matrix R-by-1 indicating whether each point in R-by-2 matrix     */
/*  POINTS lies within (value 1) or without (value 0) polygon POLY           */
/*  -----------------------------------------------------------------------  */

real matrix sp_pip(real matrix POINTS, real matrix POLY)
{


/* Setup */
real matrix    PIP
real scalar    np
real scalar    i, x, y

/* Generate working objects */
np = rows(POINTS)
PIP = J(np, 1, 0)

/* Check and report points status */
for (i=1; i<=np; i++) {
   x = POINTS[i, 1]
   y = POINTS[i, 2]
   PIP[i] = sp_spip(x, y, POLY)
}

/* Return results */
return(PIP)


}




/*  -----------------------------------------------------------------------  */
/*  sp_polyarea()                                                            */
/*                                                                           */
/*  Returns scalar containing the area of polygon POLY                       */
/*  -----------------------------------------------------------------------  */

real scalar sp_polyarea(real matrix POLY)
{


/* Setup */
real scalar    area
real scalar    nv, i, x1, y1, x2, y2

/* Generate working objects */
area = 0
nv = rows(POLY)
if (POLY[1,1] == POLY[nv,1] & POLY[1,2] == POLY[nv,2]) nv = nv - 1

/* Compute polygon area */
for (i=1; i<=nv; i++) {
   x1 = POLY[i, 1]
   y1 = POLY[i, 2]
   if (i == nv) {
      x2 = POLY[1, 1]
      y2 = POLY[1, 2]
   }
   else {
      x2 = POLY[i+1, 1]
      y2 = POLY[i+1, 2]
   }
   area = area + ((x2 - x1) * (y2 + y1))
}

/* Return results */
return(abs(area) * 0.5)


}




/*  -----------------------------------------------------------------------  */
/*  sp_spip()                                                                */
/*                                                                           */
/*  Returns scalar indicating whether single point (x,y) lies within         */
/*  (value 1) or without (value 0) polygon POLY                              */
/*  -----------------------------------------------------------------------  */

real scalar sp_spip(real scalar x, real scalar y, real matrix POLY)
{


/* Setup */
real scalar    pip
real scalar    nv, iwind, xlastp, ylastp, ioldq, inewq
real scalar    xthisp, ythisp
real scalar    i, a, b

/* Generate working objects */
nv = rows(POLY)
if (POLY[1,1] == POLY[nv,1] & POLY[1,2] == POLY[nv,2]) nv = nv - 1
iwind = 0
xlastp = POLY[nv, 1]
ylastp = POLY[nv, 2]
ioldq = sp_spip_aux(xlastp, ylastp, x, y)

/* Check and report point status */
for (i=1; i<=nv; i++) {
   xthisp = POLY[i, 1]
   ythisp = POLY[i, 2]
   inewq = sp_spip_aux(xthisp, ythisp, x, y)
   if (ioldq != inewq) {
      if (mod(ioldq+1, 4) == inewq) {
         iwind = iwind + 1
      }
      else if (mod(inewq+1, 4) == ioldq) {
         iwind = iwind - 1
      }
      else {
         a = (ylastp - ythisp) * (x - xlastp)
         b = xlastp - xthisp
         a = a + ylastp * b
         b = b * y
         if (a > b) {
            iwind = iwind + 2
         }
         else {
            iwind = iwind - 2
         }
      }
   }
   xlastp = xthisp
   ylastp = ythisp
   ioldq = inewq
}
pip = abs(iwind / 4)

/* Returns results */
return(pip)


}

/* Auxiliary function */
real scalar sp_spip_aux(real scalar xp, real scalar yp, real scalar xo,
                        real scalar yo)
{
real scalar iq
if(xp < xo) {
   if(yp < yo) iq = 2
   else iq = 1
}
else {
   if(yp < yo) iq = 3
   else iq = 0
}
return(iq)
}




/*  -----------------------------------------------------------------------  */              
/*  Exit Mata                                                                */
/*  -----------------------------------------------------------------------  */              

end

===== CODE ENDS HERE ===============================================================




Il giorno 17/lug/2013, alle ore 04.42, graupel75@gmail.com ha scritto:

> Dear Maurizio,
> 
> One small question, any suggestions on adding/keeping the respective
> polygon name for each coordinate pair in the "casecoord.dta" output
> from -csrpp- ? I tried to work it in the .ado you posted, but all
> attempts error out, and it appears that the coordinates are not
> sequential down the polygon list.
> 
> thanks again!
> 
> cheers-
> Andrew Lover
> ___________________________________________________________
> Epidemiologist
> Centre for Infectious Disease Epidemiology Research (CIDER)
> Saw Swee Hock School of Public Health
> National University of Singapore
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/


*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/