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

["Lars E. Kroll" <mail@lkroll.de>]

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/