Re: st: Aggregating spatial units for SPMAP (sorry, no simple solution)

[Scott Merryman <scott.merryman@gmail.com>]

Though this does not aggregate regions, one can use the line
suboption (along with the polyfirst option) to draw a regional
border:

Scott


cd "C:/data"
use "Italy-RegionsCoordinates.dta",clear
gen aggregated_region = 1 if inrange(_ID,1,10)
replace aggregated_region = 2 if inrange(_ID,11,20)

save foo,replace
qui {
foreach l in 8 9 10 11 12 {
use foo,clear
       keep if _ID == `l'
       duplicates drop
       save foo`l',replace
       }
//area 12 in region 2 borders areas 9 and 10 in region 1
foreach l in 9 10 {
use foo12,clear
gen order = _n
merge 1:1 _X _Y using foo`l'
keep if _m == 3
gen gr = `l'
sort order
save foo12_`l',replace
}

//Area 11 in region 2 borders areas 8, 9, and 10 in region 1
foreach l in 10 9 8   {
use foo11,clear
gen order = _n
merge m:1 _X _Y using foo`l'
keep if _m == 3
gen gr = `l'
sort order
save foo11_`l',replace
}

use foo11_10
drop if order > 102
save,replace

use foo12_9
append  using foo12_10 foo11_10 foo11_9 foo11_8
save foo2,replace
}

use "Italy-RegionsData.dta", clear

spmap pop98 using "foo", id(id)  ///
  fcolor(Greys)  os(none ..)    ///
       polygon(data("foo") by(aggregated_region) fcolor(none none) ///
       oc(red blue) osize(vthick vthick)) ///
       polyfirst line(data(foo2) color(blue) size(thick)) ///
       label(x(x) y(y) label(id)) name(region,replace)

fs foo*
   foreach f in `r(files)' {
            rm `f'


On Thu, Jun 7, 2012 at 5:22 AM, Lars E. Kroll <mail@lkroll.de> wrote:
> Dear Statalist-Members,
>
> I have follow up to the question regarding aggregating spatial units
> (Statalist thread st: Aggregating spatial units), unfortunately a
> previous simple answer on the list proved to be wrong (see code
> example below) and I'm searching for a more general solution to the
> problem.
>
> I appreciated any hints and I promise to release the final working code in SSC.
>
> Yours
> Lars
>
> PROBLEM:
> The aggregation of spatial units (polygons) is somehow working. WHAT
> IS CURRENTLY NOT WORKING is getting rid of the BORDERS WITHIN
> aggregated areas (polygons).
>
> WHY:
> Often, regional data (shape files) are freely available for some level
> of area differentiation but not for another higher area level, which
> is needed by the researcher.
>
> GOAL:
> I would like to provide/use an ADO, that can aggregate spatial units
> and drop inner borders of the aggregated (multi part) higher level
> areas. Currently we can draw high level multi part areas but have to
> keep the inner boarders of the constituting sub areas. (Higher level
> Areas are marked by BLUE and RED boarders in the example below).
>
> EXAMPLE CODE:
> version 12.1
> set more off
> capture clear
> // Load relevant ADOs
> capture net install spmap.pkg
> capture net get spmap.pkg
>
> // Define Tempfiles
> tempfile test_data
> tempfile test_coordinates_regions
> tempfile test_coordinates_aggr_regions
> set seed 123456789
>
> // Generate Regions in Master Datafile
> sysuse "Italy-RegionsData.dta", clear
> gen aggregated_region = 1 if inrange(id,1,10)
> replace aggregated_region = 2 if inrange(id,11,20)
> gen test = uniform()*10
> format test %1.0f
> save "`test_data'.dta"  , replace
>
> // Generate Regions in Coordinates Datafile
> sysuse "Italy-RegionsCoordinates.dta", clear
> * Simplify map data
> keep if (mod(_n,5) == 0) | (_X>=.)
> save "`test_coordinates_regions'.dta" , replace
> * Generate Region File
> gen aggregated_region = 1 if inrange(_ID,1,10)
> order aggregated_region
> replace aggregated_region = 2 if inrange(_ID,11,20)
> list if inlist(_ID, 1,2) , sepby(_ID ) noobs
> di as result "PLEASE FILL IN YOUR MAGIC STATA-CODE THAT ERASES THE
> INNER BOARDERS HERE!!"
> * Save new Coordinates
> save "`test_coordinates_aggr_regions'.dta" , replace
>
> // Draw the MAP
> use "`test_data'.dta" , clear
> spmap test using "`test_coordinates_regions'" , id(id)          ///
>     clnumber(5) ocolor(black ..) osize(none ..)    ///
>     title("Problem of aggregation with SPMAP-data files",
> size(*0.8))           ///
>         subtitle("RED= Region 1, BLUE = Region 2") ///
>     legstyle(3) legend(ring(1) position(3))                            ///
>     plotregion(icolor(white)) graphregion(icolor(white))               ///
>     polygon(data("`test_coordinates_aggr_regions'")
> by(aggregated_region) osize(thick ..) ocolor(red blue)   )
>
> // CLEAN UP THE MESS
> capture erase "`test_data'.dta"
> capture erase "`test_coordinates_regions'.dta"
> capture erase "`test_coordinates_aggr_regions'.dta"
> exit
>
> IDEAS FOR A SOLUTION:
> Steps:
> 1. Group polygons of regions with shared boarders (to allow for multi
> part regions)
> 2. Sort polygon within these sub regions by their centroids
> 3. The tricky part: Erase any borders that are within the sub regions
> and keep the outer boarder
>
> Helps:
> - Somehow, one must "walk" along the boarder of a polygon UNTIL there
> is another boarder of another Polygon within the same region, with a
> more "outer" boarder.
>  If you can code these criteria, you solved the main puzzle.
> - But.. maybe there are some path finding algorithms out there, that
> can do all the tricks?
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/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/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/