Multiple maps in spmap

[John MacDonald]

How do you use the code for spmap to overlay two shapefiles on top of each other?

[Friedrich Huebler]

Could you be more specific? Would merging the shapefiles be an option?

[John MacDonald]

I should have noted. I have two separate shapefiles that are in different polygons. I'd like to overlay the bigger polygons on top of the smaller ones. It looks like the spmap program has this functionality but there isn't an example in the help menu.

I can't merge the files without changing the polygons. I could merge the database file though and call up different coordinate files if that works.

[Friedrich Huebler]

Maybe you could do this as a layer, as described in this presentation by Maurizio Pisati, the author of spmap.

[wbuchanan]

John MacDonald just to make sure, are the contents in both files polygons? The Shapefile specification only allows one type of geometry per file, so if you had a file with Multipolygons and a file with Polygons it might require a bit more work to place all the data in a single file. The other issue that could cause problems is the way graphs are rendered in Stata. For example, if the larger polygons contain any color they would effectively mask the underlying smaller polygons.

[Robert Picard]

It's pretty easy to combine two shapefiles. The easiest is to negate the identifiers (e.g. replace _ID = -_ID) in one of the two shapefiles and then append the coordinates (and databases).

However spmap appears to have an option to overlay polygons. Here's a quick example. I first use mergepoly to create a shapefile of the boundaries of the whole country from the regions shapefile. To install mergepoly, type in Stata's Command window:

ssc install mergepoly

Then I overlay the outline over a regions map

* Make a shapefile for the outline of Italy; type -ssc install mergepoly-
*
use "Italy-RegionsData.dta", clear
mergepoly id using "Italy-RegionsCoordinates.dta", coord("italy_coor.dta") replace

use "Italy-RegionsData.dta", clear
spmap relig1 using "Italy-RegionsCoordinates.dta", id(id) ///
    polygon(data("italy_coor.dta") fcolor(none) ocolor(blue) osize(thick))

The resulting maps looks like

[John MacDonald]

This was very helpful and worked perfectly. Thank you Robert.
Is there a way to link polygons and join on the primary overlapping polygons? I was able to use the geoinpoly program you wrote. I simply changed the _ID variable name and then used the coordinate file like it was a point file to join to another coordinate file. This seems to have worked fine.

[Robert Picard]

Indeed, you can use geoinpoly to find which polygons from shapefile A overlap polygons from shapefile B. As you found out, you just consider the coordinates of the polygon vertices from shapefile A as points and geoinpoly will match these points to polygons from shapefile B. You can't form new polygons with the union or intersection of two shapefiles however.

[Mauricio Ramirez]

Hi, Robert. I've trying to use the instructions you gave above but didn't make it. After using mergepoly and trying to generate the map with spmap Stata retuns "file __POL.dta could not be opened" error r(603). Do you have any clue about it?

[Max Julian]

Hi Robert, I would like to create a dataset that indicates if one polygon is (partially) nested in another polygon. I tried geoinpoly, but it gives me an error message, because I need the polygons in Lat/Long format but I am not sure if the XY coordinates for a shape can easily be converted into Lat/Long format. Do you have any solution to my problem?
Thanks a lot!

[Robert Picard]

You do not have to convert your projected coordinates to lat/lon to use geoinpoly, there's a noprojection option.

[Kyle Smith]

I have a question regarding post #6 above. I am trying to do something similar to:

use "Italy-RegionsData.dta", clear
spmap relig1 using "Italy-RegionsCoordinates.dta", id(id) polygon(data("italy_coor.dta") ///
fcolor(none) ocolor(blue) osize(thick))

Specifically, I am using spmap to plot county-level data for different states. I would like the county plots to use colors, but I want to draw state boundaries in a diff color. The only change to the above (post #6) is that instead of the interior regions being different shades of gray, I would like them to be diff colors. But I would still like the exterior border (the one in blue), to be a color (not grey). I want the state borders all to be the same color. For example, when I change the last/third line of code to:

 fcolor(Blues) ocolor(green) osize(thick))

I get the following error:

 Problem with option polygon(): when no group variable is specified in suboption by(), suboption fcolor() does not accept palette names  r(198);

Any idea how to fix this?

Thanks in advance.