Workshop on “Creating maps with Stata and thematic applications” by Ben Jann, Hermosillo, October 24–25, 2023
In this exercise we want to draw maps of the world and see how different projections affect the result. Great shape data of the world can be found at https://www.naturalearthdata.com/.
1. Download the data on countries, rivers, and lakes at the 1:50m scale, convert to Stata format, and load the data into frames.
At https://www.naturalearthdata.com/, country data is available under category "Cultural", rivers and lakes are available under category "Physical".
. geoframe convert Countries data/ne_50m_admin_0_countries, replace (translating ne_50m_admin_0_countries from directory data/ne_50m_admin_0_countries) (importing .shp file) (importing .dbf file) (creating _ID spatial-unit id) (creating _CX coordinate) (creating _CY coordinate) file Countries_shp.dta created file Countries.dta created (type geoframe create Countries to load the data) . geoframe convert Rivers data/ne_50m_rivers_lake_centerlines, replace (translating ne_50m_rivers_lake_centerlines from directory data/ne_50m_rivers_lake_centerlines) (importing .shp file) (importing .dbf file) (creating _ID spatial-unit id) (creating _CX coordinate) (creating _CY coordinate) file Rivers_shp.dta created file Rivers.dta created (type geoframe create Rivers to load the data) . geoframe convert Lakes data/ne_50m_lakes, replace (translating ne_50m_lakes from directory data/ne_50m_lakes) (importing .shp file) (importing .dbf file) (creating _ID spatial-unit id) (creating _CX coordinate) (creating _CY coordinate) file Lakes_shp.dta created file Lakes.dta created (type geoframe create Lakes to load the data) . . geoframe create Countries, replace (reading shapes from Countries_shp.dta) (all observations in frame Countries_shp matched) (link to frame Countries_shp added) (current frame now Countries) Frame name: Countries Frame type: unit Feature type: <none> Number of obs: 242 Unit ID: _ID Coordinates: _CX _CY Area: <none> Linked shape frame: Countries_shp . geoframe create Rivers, replace feature(water) (reading shapes from Rivers_shp.dta) (all observations in frame Rivers_shp matched) (link to frame Rivers_shp added) (current frame now Rivers) Frame name: Rivers Frame type: unit Feature type: water Number of obs: 478 Unit ID: _ID Coordinates: _CX _CY Area: <none> Linked shape frame: Rivers_shp . geoframe create Lakes, replace feature(water) (reading shapes from Lakes_shp.dta) (all observations in frame Lakes_shp matched) (link to frame Lakes_shp added) (current frame now Lakes) Frame name: Lakes Frame type: unit Feature type: water Number of obs: 412 Unit ID: _ID Coordinates: _CX _CY Area: <none> Linked shape frame: Lakes_shp
2. Draw a world map with countries colored by type of economy and also include rivers and lakes. Exclude Antarctica from the map. Furthermore, color the background of the map (i.e. the sea) in the same color as lakes.
Variable ECONOMY contains the type of economy. Generate a numeric
variable by applying command encode to the variable
(and maybe clean the labels).
You can used global option plotregion(color()) to set the background color. The
default color for lakes used by geoplot is color("135 206 235*.5").
To exclude Antarctica, apply geoframe select if SOVEREIGNT!="Antarctica" to
frame Countries.
. // fix ECONOMY variable . frame change Countries . fre ECONOMY ECONOMY -- ECONOMY -------------------------------------------------------------------------------- | Freq. Percent Valid Cum. -----------------------------------+-------------------------------------------- Valid 1. Developed region: G7 | 7 2.89 2.89 2.89 2. Developed region: nonG7 | 49 20.25 20.25 23.14 3. Emerging region: BRIC | 4 1.65 1.65 24.79 4. Emerging region: MIKT | 4 1.65 1.65 26.45 5. Emerging region: G20 | 20 8.26 8.26 34.71 6. Developing region | 104 42.98 42.98 77.69 7. Least developed region | 54 22.31 22.31 100.00 Total | 242 100.00 100.00 -------------------------------------------------------------------------------- . encode ECONOMY, gen(economy) . forv i = 1/7 { 2. local lbl: label economy `i' 3. local lbl = substr(`"`lbl'"',4,.) 4. lab def economy `i' `"`lbl'"', modify 5. } . fre economy economy -- ECONOMY ------------------------------------------------------------------------------- | Freq. Percent Valid Cum. ----------------------------------+-------------------------------------------- Valid 1 Developed region: G7 | 7 2.89 2.89 2.89 2 Developed region: nonG7 | 49 20.25 20.25 23.14 3 Emerging region: BRIC | 4 1.65 1.65 24.79 4 Emerging region: MIKT | 4 1.65 1.65 26.45 5 Emerging region: G20 | 20 8.26 8.26 34.71 6 Developing region | 104 42.98 42.98 77.69 7 Least developed region | 54 22.31 22.31 100.00 Total | 242 100.00 100.00 ------------------------------------------------------------------------------- . // get rid of Antarctica . geoframe select if SOVEREIGNT!="Antarctica" (dropped 5,031 observations in frame Countries_shp) (dropped 1 observation in frame Countries) . // draw map . geoplot /// > (area Countries i.economy, color(viridis)) /// > (area Countries) /// > (line Rivers, lw(0.1)) /// > (area Lakes) /// > , legend(pos(sw) reverse rowgap(1) region(color(gs15))) /// > plotregion(color("135 206 235*.5")) tight
3. Generate a frame containing a grid using geoframe grid and add the
grid to the map.
Use a grid from -180 to 180 in steps of 15 degrees for the X axis and from -65 to 85 in steps of 15 degrees for the Y axis
. geoframe grid Grid, x(-180(15)180) y(-65(15)85) (new frame Grid created) (new frame Grid_shp created) . geoplot /// > (area Countries i.economy, color(viridis)) /// > (area Countries) /// > (line Rivers, lw(0.1)) /// > (area Lakes) /// > (line Grid) /// > , legend(pos(sw) reverse rowgap(1) region(color(gs15))) /// > plotregion(color("135 206 235*.5")) tight
4. Transform the data using the Web Mercator projection using geoframe project
and draw the transformed map.
You need to apply geoframe project to all frames
(Countries, Rivers, Lakes, Grid). Also make sure to store the transformed data
in new frames using option into(), so that the original data is preserved.
. frame Countries: geoframe project, into(pCountries) (new frame pCountries created) (new shape frame pCountries_shp created) (projecting _CX and _CY in frame pCountries using web_mercator) (projecting _X and _Y in frame pCountries_shp using web_mercator) . frame Rivers: geoframe project, into(pRivers) (new frame pRivers created) (new shape frame pRivers_shp created) (projecting _CX and _CY in frame pRivers using web_mercator) (_CX has values outside [-180,180]; using clipped values) (_CY has values outside [-90,90]; using clipped values) (projecting _X and _Y in frame pRivers_shp using web_mercator) . frame Lakes: geoframe project, into(pLakes) (new frame pLakes created) (new shape frame pLakes_shp created) (projecting _CX and _CY in frame pLakes using web_mercator) (projecting _X and _Y in frame pLakes_shp using web_mercator) . frame Grid: geoframe project, into(pGrid) (new frame pGrid created) (new shape frame pGrid_shp created) (projecting _CX and _CY in frame pGrid using web_mercator) (projecting _X and _Y in frame pGrid_shp using web_mercator) . geoplot /// > (area pCountries i.economy, color(viridis)) /// > (area pCountries) /// > (line pRivers, lw(0.1)) /// > (area pLakes) /// > (line pGrid) /// > , legend(pos(sw) reverse rowgap(1) region(color(gs15))) /// > plotregion(color("135 206 235*.5")) tight
5. Now draw a map using the Robinson projection.
. frame Countries: geoframe project robinson, into(pCountries) replace (new frame pCountries created) (new shape frame pCountries_shp created) (projecting _CX and _CY in frame pCountries using robinson) (projecting _X and _Y in frame pCountries_shp using robinson) . frame Rivers: geoframe project robinson, into(pRivers) replace (new frame pRivers created) (new shape frame pRivers_shp created) (projecting _CX and _CY in frame pRivers using robinson) (_CX has values outside [-180,180]; using clipped values) (_CY has values outside [-90,90]; using clipped values) (projecting _X and _Y in frame pRivers_shp using robinson) . frame Lakes: geoframe project robinson, into(pLakes) replace (new frame pLakes created) (new shape frame pLakes_shp created) (projecting _CX and _CY in frame pLakes using robinson) (projecting _X and _Y in frame pLakes_shp using robinson) . frame Grid: geoframe project robinson, into(pGrid) replace (new frame pGrid created) (new shape frame pGrid_shp created) (projecting _CX and _CY in frame pGrid using robinson) (projecting _X and _Y in frame pGrid_shp using robinson) . . geoplot /// > (area pCountries i.economy, color(viridis)) /// > (area pCountries) /// > (line pRivers, lw(0.1)) /// > (area pLakes) /// > (line pGrid) /// > , legend(pos(sw) reverse rowgap(1) region(color(gs15))) /// > plotregion(color("135 206 235*.5")) tight
A problem with this map is that the blue background goes beyond the edge of the globe. Here is how you can fix this. Trick is to generate a convex hull around the projected grid and then add the hull as a layer to the graph instead of setting a background color for the plotregion.
. frame pGrid: geoframe bbox Globe, hull (new frame Globe created) . frame Globe: geoframe set feature water . . geoplot /// > (area Globe) /// > (area pCountries i.economy, color(viridis)) /// > (area pCountries) /// > (line pRivers, lw(0.1)) /// > (area pLakes) /// > (line pGrid) /// > , legend(pos(sw) reverse rowgap(1) region(color(gs15))) tight
