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 |
Alexandru Cojocaru <scojocaru@gmail.com> |

To |
statalist@hsphsun2.harvard.edu |

Subject |
st: Aggregating spatial units for SPMAP |

Date |
Wed, 26 Dec 2012 10:34:30 -0500 |

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/

**Follow-Ups**:**Re: st: Aggregating spatial units for SPMAP***From:*Robert Picard <picard@netbox.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>

- Prev by Date:
**Re: st: RE: replace missing values in different waves** - Next by Date:
**RE: st: Re: margins after xtlogit,fe** - Previous by thread:
**Re: st: Aggregating spatial units for SPMAP** - Next by thread:
**Re: st: Aggregating spatial units for SPMAP** - Index(es):