Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Merging two provinces into one while working with maps (shapefile)

    Hi everyone

    I am trying to draw a map from 1928, but the shapefile that I have is from 2016. The new map contains more provinces (states) than the division that I have in 1928.

    For example: what today are two different provinces, back then was just one. I am trying to merge the two regions in order to have one single delimitation. But I have not succeed so far, I used the mergepoly but the only thing that it does is merging all regions into one.

    Does someone happen to know how could I solve this inconvenient?
    I attached the map, and the regions selected are some of the provinces that I would like to merge.

    Thanks









    Attached Files

  • #2
    Here's an example that shows how to merge midwest states within the map of the conterminous 48 US states. The example requires mergepoly, geo2xy, and spmap, all from SSC. You will also need to install the geo2xy ancillary datasets as well as a small dataset of US regions using the following commands:

    Code:
    net get geo2xy, from("http://fmwww.bc.edu/repec/bocode/g")
    copy "http://robertpicard.com/stata/us_regions.dta" "us_regions.dta"
    The important detail to understand is that a shapefile has two parts stored in separate datasets:
    1. the database (with one observation per feature)
    2. the coordinates for each feature (one observation per point; connecting consecutive points creates the feature's polygon).
    So to get there, you need to create a new shapefile with a unique identifier for each merged region in the database. The coordinates dataset must include the location of points that describes the merged region(s). This requires a fair amount of data management gymnastics and you must be careful not to mess up the order point points that describe each feature:

    Code:
    * reduce to conterminous 48 states plus DC, add obs counter to coordinates
    use "geo2xy_us_data.dta", clear
    list if inlist(_ID,14,39,42)
    drop if inlist(_ID,14,39,42)
    save "us48_db.dta", replace
    use "geo2xy_us_coor.dta", clear
    gen long obs = _n
    drop if inlist(_ID,14,39,42)
    save "us48_coor.dta", replace
    
    * Identify midwest states and use 9999 as the identifier to use
    use "us48_db.dta", clear
    rename STATEFP STATE
    merge 1:1 STATE using "us_regions.dta", nogen
    keep if region == "Midwest Region":region
    gen long _IDmidwest = 9999
    
    * we now have a database of one obs per state in the Midwest
    isid _ID, sort
    list _IDmidwest _ID NAME region
    save "midwest_states.dta", replace
    
    * merge all midwest states
    mergepoly using "us48_coor.dta", coor("midwest_coor.dta") replace
    list
    save "midwest_db.dta", replace
    
    * go back to the original US States shapefile and remove midwest states
    use "us48_db.dta", clear
    merge 1:1 _ID using "midwest_states.dta", keep(master) nogen
    
    * append the single database record for the midwest and save the new combo db
    append using "midwest_db.dta"
    replace _ID = _IDmidwest if !mi(_IDmidwest)
    replace NAME = "Midwest States" if !mi(_IDmidwest)
    isid _ID, sort
    save "us+mid_db.dta", replace
    
    * adjust the midwest coor identifier to match and append to US coor
    use "midwest_coor.dta", clear
    replace _ID = 9999
    gen long obs = _n
    save "midwest_coor.dta", replace
    use "us48_coor.dta", clear
    append using "midwest_coor.dta"
    sort _ID obs
    save "us+mid_coor.dta", replace
    
    * never make a map from geographich coordinates, you must use a map projection
    geo2xy _Y _X, replace projection(albers)
    save "us+mid_coor_XY.dta", replace
    
    * finally, make a map of the US with all midwest states represented by a single polygon
    use "us+mid_db.dta", clear
    spmap using "us+mid_coor_XY.dta", id(_ID)
    
    graph export "us_midwest.png", width(800) replace
    Here's the listing of midwest states that are grouped:
    Code:
    . list _IDmidwest _ID NAME region
    
         +------------------------------------------------+
         | _IDmid~t   _ID           NAME           region |
         |------------------------------------------------|
      1. |     9999     5      Wisconsin   Midwest Region |
      2. |     9999    11       Michigan   Midwest Region |
      3. |     9999    16      Minnesota   Midwest Region |
      4. |     9999    17       Nebraska   Midwest Region |
      5. |     9999    19           Ohio   Midwest Region |
         |------------------------------------------------|
      6. |     9999    20       Illinois   Midwest Region |
      7. |     9999    21       Missouri   Midwest Region |
      8. |     9999    22           Iowa   Midwest Region |
      9. |     9999    23   South Dakota   Midwest Region |
     10. |     9999    32        Indiana   Midwest Region |
         |------------------------------------------------|
     11. |     9999    48         Kansas   Midwest Region |
     12. |     9999    51   North Dakota   Midwest Region |
         +------------------------------------------------+
    and here's the map generated

    Click image for larger version

Name:	us_midwest.png
Views:	1
Size:	74.5 KB
ID:	1521677


    I plan to clean up this example and include it in the next update to mergepoly.

    Comment


    • #3
      Hi Robert, thank you for your answer!

      I did all the steps but came across the following error:

      * never make a map from geographic coordinates, you must use a map projection
      geo2xy _Y _X, replace project(albers)
      save "coor_XY.dta", replace

      *ERROR:
      geo2xy _Y _X, project(albers)
      latitude _Y must be between -90 and 90

      I continued running the commands but in the end, the result was not good, I guess because of that error.

      Comment


      • #4
        I assume that you are trying to adapt the code to your situation and so the error suggests that you have projected coordinates and not lat/lon. Just skip that step and call spmap using the coordinates dataset equivalent to "us+mid_coor.dta" in the example.

        Comment


        • #5
          Robert Picard This -mergepoly- command is So Great Job! Thank you!

          Comment


          • #6
            @Robert Picard, thanks a lot i had the same problem and i could solved it whit your commands!!!

            Comment

            Working...
            X