Notice: On March 31, it was **announced** that Statalist is moving from an email list to a **forum**. The old list will shut down on April 23, and its replacement, **statalist.org** is already up and running.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From |
Robert Picard <picard@netbox.com> |

To |
"statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |

Subject |
Re: st: Aggregating spatial units for SPMAP |

Date |
Wed, 26 Dec 2012 19:22:36 -0500 |

Looks like -mergepoly- created an enclave within the borders of Romania. If that's the case, then the only explanation I can think of is that the shapefile you are using is missing one or more county polygons. If all you want is the outer borders of Romania, then I would suggest you just drop the polygon that forms the enclave. On Dec 26, 2012, at 10:39 AM, Alexandru Cojocaru <scojocaru@gmail.com> wrote: > Dear Robert, > > [not changing the subject line as my problem is the same as in the > original post] > > I tried -mergepoly- (available from SSC) on a shapefile containing the > 41 regions of Romania. It produces the outer borders of Romania, but > also the inner (to Romania) border for 1 of the 41 counties. > > Any thoughts? > Sandu > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > Dear Robert, > > Just tried mergepoly (available from SSC) on a shapefile containing > the 41 counties of Romania. It produces the outer border of Romania, > but the graph also traces the internal (to Romania) border of 1 of the > regions. > > Any thoughts? > Sandu > > On Sat, Dec 22, 2012 at 12:22 PM, Robert Picard <picard@netbox.com> wrote: >> >> Thanks to Kit Baum, a much improved version of -mergepoly- with help >> file is now available through SSC. To install it, type -ssc install >> mergepoly- in Stata. Here's the program description: >> >> mergepoly is used to merge adjacent polygons from a shape >> boundary file. Shape boundary files are usually obtained from >> outside sources and converted to Stata datasets using -shp2dta- >> (SSC). The shape dataset contains one record per geographic feature and >> the >> coordinate dataset describes the shape of each feature in terms >> of points. The points, in order, form polygons. mergepoly creates >> new polygons that follow the outer border of the adjacent polygons >> by removing all shared line segments. >> >> Enjoy, Robert >> >> >> On Sun, Dec 16, 2012 at 3:33 PM, Robert Picard <picard@netbox.com> wrote: >>> I think I may have a solution. The following program takes a >>> shapefile and creates new polygons that have the shape of the >>> outer borders of the adjacent polygons. It builds on the fact >>> that adjacent polygons have common points along line segments >>> that are shared. The outer borders are composed of polylines >>> (connected line segments) that start and end with a shared point >>> with all in-between points unique to a single polygon. It's a bit >>> tricky to connect all the outer border polylines, but minimal >>> testing suggests that the program works. >>> >>> Let me know if the program works or not in your case. If no >>> problem is reported, I'll add bells and whistles and write up a >>> help file and post it to SSC. >>> >>> Hope this helps, Robert >>> >>> * -------------------------------------------------------------- >>> version 12 >>> clear all >>> >>> *! 1.0.0 16dec2012 Robert Picard, picard@netbox.com >>> program define mergepoly >>> >>> args id0 lat0 lon0 >>> >>> rename `lat0' lat >>> rename `lon0' lon >>> >>> gen obs = _n >>> >>> // create new id in case there are multipart polygons >>> gen id = sum(mi(lat,lon)) >>> sort id obs >>> qui by id: drop if _n == 1 >>> >>> // each polygon should start and end at the same lat/lon >>> by id: assert lat[1] == lat[_N] >>> by id: assert lon[1] == lon[_N] >>> >>> // adjacent polygons have common points >>> sort lat lon obs >>> by lat lon: gen N = _N >>> >>> // identify polylines that form the outside borders; >>> // they start and end at common points and have a >>> // run of non-adjacent points >>> sort id obs >>> by id: gen pstart = N > 1 & N[_n+1] == 1 >>> by id: gen pend = N > 1 & N[_n-1] == 1 >>> by id: gen inpoly = sum(pstart - pend) >>> qui replace inpoly = 1 if pend >>> gen polyline = sum(pstart) * inpoly >>> >>> // link polylines by grouping coordinates and >>> // identifying the next polyline >>> sort lat lon pend pstart obs >>> by lat lon: egen pp = total(pstart + pend) >>> assert pp <= 2 >>> by lat lon: assert pend if pp == 2 & _n == _N >>> by lat lon: assert pstart if pp == 2 & _n == _N - 1 >>> qui by lat lon: gen next = polyline[_n-1] if _n == _N & pp == 2 >>> >>> // handle special case of line segments on the outer border; >>> // these are not polylines, they just have 2 points and they >>> // create a bridge between two polylines. >>> sort lat lon pend obs >>> by lat lon: gen sameid = id == id[_N] >>> sort lat lon pend sameid obs >>> qui by lat lon: gen link = obs[1]+1 if pend & pp == 1 >>> sort lat lon pstart obs >>> qui by lat lon: replace link = obs[1] if pstart & pp == 1 >>> sort pp link pstart pend >>> qui by pp link: replace next = polyline[2] if pp == 1 & pend >>> qui by pp link: assert _N == 2 if pp == 1 & !mi(link) >>> >>> qui keep if polyline >>> sort polyline obs >>> by polyline: gen one = _n == 1 >>> >>> // some polylines may be polygons already >>> by polyline: gen ispolyline = lat[1] != lat[_N] | lon[1] != >>> lon[_N] >>> sort ispolyline polyline obs >>> qui gen newpolygon = sum(one) if !ispolyline >>> >>> // loop over each polyline in order and build new polygons that >>> // follow the outer borders of adjacent polygons >>> sum newpolygon, meanonly >>> local curpolygon = cond(mi(r(max)), 1, r(max) + 1) >>> sum polyline if mi(newpolygon), meanonly >>> local nextpl = r(min) >>> local more = !mi(r(N)) >>> local startpline `nextpl' >>> >>> qui gen pline_order = . >>> >>> local i 0 >>> while `more' { >>> >>> local i = `i' + 1 >>> qui replace pline_order = `i' if polyline == `nextpl' >>> qui replace newpolygon = `curpolygon' if polyline == >>> `nextpl' >>> >>> sum next if polyline == `nextpl', meanonly >>> if `nextpl' == r(max) local startpline = r(max) >>> local nextpl = r(max) >>> >>> if `nextpl' == `startpline' { >>> local curpolygon = `curpolygon' + 1 >>> sum polyline if mi(newpolygon), meanonly >>> local nextpl = r(min) >>> local startpline `nextpl' >>> } >>> >>> qui count if mi(newpolygon) >>> local more = r(N) >>> } >>> >>> // each new polygon must start with missing lat/lon >>> sort newpolygon pline_order obs >>> qui expand 2 if one >>> sort newpolygon pline_order obs >>> qui by newpolygon: replace lat = . if _n == 1 >>> qui by newpolygon: replace lon = . if _n == 1 >>> >>> rename newpolygon ext`id0' >>> keep `id0' lat lon ext`id0' >>> rename lat `lat0' >>> rename lon `lon0' >>> >>> end >>> >>> use "http://fmwww.bc.edu/repec/bocode/i/Italy-RegionsCoordinates.dta";, >>> clear >>> >>> mergepoly _ID _Y _X >>> >>> line _Y _X, lwidth(vvthin) lcolor(blue) cmissing(n) /// >>> aspectratio(1.2) xscale(off) yscale(off) /// >>> ylabel(minmax, nogrid) xlabel(minmax) >>> >>> * -------------------------------------------------------------- >>> >>> >>> >>> >>> On Fri, Dec 14, 2012 at 6:11 AM, Karsten Pfaff >>> <atze1985@zedat.fu-berlin.de> wrote: >>>> Dear Statalist members, >>>> I would like to refer first to the following thread that emerged some >>>> month earlier on statalist: >>>> >>>> http://www.stata.com/statalist/archive/2012-06/msg00384.html >>>> >>>> My problem is similar. I have a shapefile (coordinate file) containing >>>> polygon coordinates about a certain type of regions (calles "kreise") >>>> in >>>> Germany, which allows me to draw a grid of Germany that divides Germany >>>> in >>>> just these regions. However, I am not interested in the regions >>>> themselves >>>> but in some kind of higher order regions called "ror". Of course I also >>>> have information about which "kreise" are part of which "ror" (just >>>> imagine you have coordinates about the counties of the USA and based on >>>> these coordinates you would like to get the coordinates of the federal >>>> states, knowing which countries belong to which federal state). The >>>> data >>>> file looks as follows (illustrative example): >>>> >>>> ror kreise (=_ID) _Y _X >>>> 101 1 9.0468785 54.872552 >>>> 101 1 9.0646343 54.871183 >>>> 101 1 9.0649463 54.871159 >>>> 101 1 9.090514 54.70260 >>>> 101 2 >>>> 101 2 >>>> 101 2 some coordinates for kreise 2 & 3 >>>> 101 2 >>>> 101 2 >>>> 101 2 >>>> 102 3 >>>> 102 3 >>>> >>>> So obviously, ror 101 contains kreise 1 and 2. >>>> Based on these information, I would like to create a shapefile that >>>> draws >>>> polygones of each ror, and ignores the "within" borders. Is there any >>>> possibility to do so? I think there must be some way since all the >>>> necessary information to create a map containing but the rors are at >>>> hand. >>>> Any help or suggestion is highly appreciated. >>>> >>>> Best regards Karsten Pfaff >>>> >>>> >>>> >>>> >>>> * >>>> * 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/ > * > * 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/

**Follow-Ups**:**Re: st: Aggregating spatial units for SPMAP***From:*Alexandru Cojocaru <scojocaru@gmail.com>

**References**:**st: Aggregating spatial units for SPMAP***From:*"Karsten Pfaff" <atze1985@zedat.fu-berlin.de>

**Re: st: Aggregating spatial units for SPMAP***From:*Robert Picard <picard@netbox.com>

**Re: st: Aggregating spatial units for SPMAP***From:*Robert Picard <picard@netbox.com>

**st: Aggregating spatial units for SPMAP***From:*Alexandru Cojocaru <scojocaru@gmail.com>

- Prev by Date:
**st: Tobit graphics** - Next by Date:
**Re: st: Aggregating spatial units for SPMAP** - Previous by thread:
**st: Aggregating spatial units for SPMAP** - Next by thread:
**Re: st: Aggregating spatial units for SPMAP** - Index(es):