Bookmark and Share

Notice: On March 31, it was announced that Statalist is moving from an email list to a forum. The old list will shut down at the end of May, and its replacement, statalist.org is already up and running.


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

st: Aggregating spatial units for SPMAP


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/


© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   Site index