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 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]

Re: st: Aggregating spatial units for SPMAP


From   Alexandru Cojocaru <scojocaru@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Aggregating spatial units for SPMAP
Date   Thu, 27 Dec 2012 09:39:29 -0500

Dear Robert,

The issue does indeed appear to be with the shapefiles I had. I
downloaded shapefiles for Romania and Lithuania from a different
location, and for both countries -mergepoly- produces the intended
result, i.e. the outer borders of the country.

Thanks again,
Sandu

On Wed, Dec 26, 2012 at 8:59 PM, Robert Picard <picard@netbox.com> wrote:
> Unfortunately, I do not have access to the shapefiles that you are
> working with. If you are interested in macro-regions of Romania, you
> should use -mergepoly- on each of these macro-regions separately. In
> other words you should keep all _ID that are within a specific
> macro-region and use -mergepoly- to remove line segments that are
> shared between adjacent counties within that macro-region.
>
> The only explanation I can offer for why you have line segments in
> Lithuania after running -mergepoly- is that these segments are not
> shared between adjacent regions within the shapefile. If the input
> shapefile is not constructed such that each region is composed of
> shared line segments (for example if each region overlaps neighbor
> regions), except on the outer borders, then -mergepoly- will not be
> able to merge such polygons.
>
> On Wed, Dec 26, 2012 at 7:49 PM, Alexandru Cojocaru <scojocaru@gmail.com> wrote:
>> Dear Robert,
>>
>> Thanks much for your reply. I looked into the shapefile, and it
>> doesn't appear like it is missing a polygon (the _ID variable has 41
>> distinct values, equal to the number of counties). I have since also
>> tried -mergepoly- on Lithuania (I have a shapefile with 10 regions),
>> and there  after -mergepoly- in addition to the outer border of
>> Lithuania, also drawn are partial inner borders for several of the ten
>> regions. Again, all ten regions appear to be present in the shapefile.
>>
>> In terms of the desired outcome, it is not the outer border of
>> Romania, but rather a map of 8 macro-regions of Romania instead of 41
>> counties that I am after. (Lithuania is not of interest, per se, just
>> wanted to see if I can replicate the problem I'm getting with Romanian
>> data).
>>
>> As I am new to this, any suggestions you may have would be most useful.
>>
>> Cheers,
>> Sandu
>>
>>
>> On Wed, Dec 26, 2012 at 7:22 PM, Robert Picard <picard@netbox.com> wrote:
>>> 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/
>> *
>> *   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