Sometimes it is difficult to know what parameters to use when projecting map data onto Esri or OpenStreetMaps in PROC SGMAP. The shapefile .PRJ file contains everything you need to set these parameters in PROC GPROJECT.

In a previous blog, you saw how to lookup projections at a web site, so they can be used with Esri or OpenStreetMaps. But, what if you can’t find the projection you need?

In this example, you are mapping boundaries for British Columbia using the PROC SGMAP CHOROMAP statement using projections you define yourself.

Here is the complete code for this example.

Import the Shapefile

Download the 2011 shapefile for British Columbia under “GIS Downloads” on this web site, then import the data using PROC MAPIMPORT and use PROC SGMAP.

proc mapimport out=BC1
  datafile="u:\canada\CD_2011.shp";
  id cdname;
run;
 
title1 "British Columbia Choromap with Non-Projected Data";
proc sgmap mapdata=BC1 noautolegend;
  openstreetmap;
  choromap / mapid=cdname; 
run;

Here is the output:

That looks good.

Next, try PROC SGMAP with an OpenStreetMap background:

title1 "British Columbia Choromap with Non-Projected Data";
proc sgmap mapdata=BC1 noautolegend;
  openstreetmap;
  choromap / mapid=cdname; 
run;

Where did your map areas go?

Use PROC GPROJECT

There are many different map projections in use across the globe. The British Columbia shapefile is using a different map coordinate system and cannot be projected onto an OpenStreetMap.

To fix this, you need to use PROC GPROJECT like this:

proc gproject data=BC1 out=BC
  from="parameters from the .PRJ file"
  to="EPSG:4326";
  id cdname;
run;

You need to supply the “from” information above, but what if you cannot find the correct projection parameters from web sites or from the SASHELP.PROJ4DEF dataset?

Build your own using data from the .PRJ file included with the shapefile package. Here is a web page that outlines the details about these parameters.

For our uses, the parameters can be summarized as follows:

+proj = PROJCS
+dataum = DATUM
+ellps = Spheroid
+x_0 = False_Easting
+y_0 = False_Northing
+lon_0 = Central_Merdian
+lat_0 = Latitude_Of_Origin
+lat_1 = Standard_Parallel_1
+lat_2 = Standard_Parallel_2
+units = UNIT
+no_defs

Here is the .PRJ file for your data:

Using the .PRJ file and definitions on the web page, fill in the parameters and include them in PROC GPROJECT:

proc gproject data=BC1 out=BC
  from="+proj=aea +datum=NAD83 +ellps=GRS80
        +x_0=1000000.0 +y_0=0.0
        +lon_0=-126.0 
        +lat_0=45.0 +lat_1=50.0 +lat_2=58.5      
        +units=m +no_defs"
  to="EPSG:4326";  /* always use EPSG:4326 for OSM and Esri */
  id cdname;
run;

Run PROC SGMAP again using the projected data and the boundaries now match the OpenStreetMap:

That is what you expected. The boundaries match the background map.

Do you see the dark jagged edges of the barrier islands on the coast? In the next blog you will work on making those easier to visualize by reducing the number of points used to draw the lines.