Combining shapefiles of two adjacent countries

[Roberto Liebscher]

Dear Statalisters,

I have shapefiles of two adjacent countries (Germany and Switzerland) which I would like to combine into one large shapefile. Is there an easy way to achieve this?

This is how I started to convert the two single shapefiles into Stata-files but I have no idea how to proceed from here. Any help is highly appreciated.

clear
cd C:\Users\wwa594\Documents\Stata\spmap

//Germany
copy "http://biogeo.ucdavis.edu/data/gadm2.8/shp/DEU_adm_shp.zip" "DEU_adm_shp.zip", replace
unzipfile "DEU_adm_shp.zip"
shp2dta using DEU_adm0, database(gerdb) coordinates(gercoord) genid(id)
use gerdb, clear
spmap using gercoord, id(id)

//Switzerland
copy "http://biogeo.ucdavis.edu/data/gadm2.8/shp/CHE_adm_shp.zip" "CHE_adm_shp.zip", replace
unzipfile "CHE_adm_shp.zip"
shp2dta using CHE_adm0, database(chfdb) coordinates(chfcoord) genid(id2)
use chfdb, clear
spmap using chfcoord, id(id2)

[William Lisowski]

Let me make a hopelessly naive suggestion, because I have no experience with shapefiles, and only the faintest understanding.

Is it possible that the coordinates of the Swiss-German border are the precisely identical in both files? In that case, if something like the following pseudo code will identify every x-y coordinate pair that appears in both files and drop them from the combined file.

in gerdb, create source=1 and rename the x and y coordinates to x and y
in chfdb, create source=2 and rename the x and y coordinates to x and y
use gerdb
append using chfdb
bysort x y (source): drop if source[1]!=source[_N]

I hope this might prove helpful if an authoritative answer from a knowledgeable user does not appear.

[Roberto Liebscher]

Thanks William Lisowski, this put me on the right track.

The code:

clear
cd C:\Users\wwa594\Documents\Stata\spmap

//Germany
copy "http://biogeo.ucdavis.edu/data/gadm2.8/shp/DEU_adm_shp.zip" "DEU_adm_shp.zip", replace
unzipfile "DEU_adm_shp.zip", replace
shp2dta using DEU_adm1, database(gerdb) coordinates(gercoord) genid(id) replace
use gerdb, clear
spmap using gercoord, id(id)

//Switzerland
copy "http://biogeo.ucdavis.edu/data/gadm2.8/shp/CHE_adm_shp.zip" "CHE_adm_shp.zip", replace
unzipfile "CHE_adm_shp.zip", replace
shp2dta using CHE_adm1, database(chfdb) coordinates(chfcoord) genid(id) replace
use chfdb, clear
spmap using chfcoord, id(id)

//combine the two
use gercoord, clear
sum _ID
local a = `r(max)'
save gercoord, replace

use chfcoord, clear
replace _ID = _ID+`a'
save chfcoord2, replace

// Databases
use gerdb, clear
local a = _N
use chfdb, clear
replace id = id+`a'
append using gerdb
save alldb, replace

// Combine Coordinates
use gercoord, clear
append using chfcoord2
sort _ID
save allcoord, replace

use alldb, clear
sort id
spmap using allcoord, id(id)

The result:

[William Lisowski]

Ah, I see I misunderstood - my solution assumed the shapefiles only had external boundaries and your goal was a map of the combined perimeter of German and Switzerland with no border between the two. So now seeing your objective I realize the problem was the duplication of values of ID in the two files, and the solution was pre-processing the IDs before appending the pairs of files. And it's clear my solution suffered from having forgotten that "shape files" include two components, a database and a coordinates file.

Thank you for sharing your solution, and congratulations on your ability to discern the outlines of a solution in my naive suggestion.

[Roberto Liebscher]

Now I am facing a different problem. I wanted to add some points (cities) to the map and used Robert Picard's geo2xy to convert latitude and longitude into the (x, y) space used for the map. However, for Munich I get the following:

city	lat	lon	xcoord	ycoord
München	48.13512	11.58198	136.23608	-88.8455

The point (136.236, -88.8455) is way off the map. In fact, for my "basemap" all y-coordinates are positive. Can someone give me a hint what went wrong here?

[Roberto Liebscher]

I solved this issue by inserting the WGS84 values by hand using this website. However, I would still be interested in understanding why geo2xy did not give me the correct results. Anyways, I wish you guys a good start into the new year.

[William Lisowski]

Looking at help geo2xy we see that a number of projections are available, where each projection defines a different mapping of (lon,lat) to (x,y). And clicking on any of them brings us to an entry in the output of help geo2xy_proj that describes the projection. Within each projection, the section on "Spheroid and (x,y) coordinates units" describes how the projection is done.

Reading this briefly, you (and I) were expecting that (lon,lat)=(0,0) would project onto (x,y) (0,0). This is not the case for the web-mercator projection that is the default.

I'm in way over my head here, but it seems to me that one of the other projections, perhaps in combination with scaling, would give you what you need.

[Robert Picard]

Roberto, you do not show code of how you tried to use geo2xy so it's hard to see what's wrong and the same goes for how you fixed the problem.

Since you appear to use spmap to generate your maps, here's an example of how to generate a map that shows the location of two cities, all using unprojected coordinates:

* Example generated by -dataex-. To install: ssc install dataex
clear
input str8 city float(lat lon)
"München" 48.13512 11.58198
"Bern"     46.94792 7.444608
end
save "cities.dta", replace

use "alldb.dta", clear
spmap using "allcoord.dta", id(id) ///
    point(data("cities.dta") xcoord(lon) ycoord(lat) fcolor(blue))

and the resulting map

Here's the same example, this time with all coordinates projected using geo2xy. By default, the projection is the same as the one used in Google Maps.

* Example generated by -dataex-. To install: ssc install dataex
clear
input str8 city float(lat lon)
"München" 48.13512 11.58198
"Bern"     46.94792 7.444608
end
geo2xy lat lon , gen(ycoord xcoord)
save "cities.dta", replace

use "allcoord.dta"
geo2xy _Y _X, replace
save "allcoord_xy.dta", replace

use "alldb.dta", clear
spmap using "allcoord_xy.dta", id(id) ///
    point(data("cities.dta") xcoord(xcoord) ycoord(ycoord) fcolor(red))

graph export map2.png, width(800)

and the resulting map

As far as I can tell, both cities are correctly located on the map in both cases.

[Roberto Liebscher]

Thanks Robert Picard. The problem was that I did not "project" the original coordinates in "allcoord.dta". I thought that these points were already "projected" and tried to add the values I obtained from geo2xy. Thanks for clearing this up, Robert.