Bookmark and Share

Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.


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

Re: st: Aggregating spatial units for SPMAP


From   Robert Picard <[email protected]>
To   "[email protected]" <[email protected]>
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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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, [email protected]
>>> 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
>>> <[email protected]> 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/


© Copyright 1996–2018 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   Site index