Imported shapefile boundaries not aligning with US outline

[bgladd]

Hello,

I created a choropleth map using an imported shapefile and it's missing some portions of state boundaries. I have tried to add an outline of the US to use as an anno layer in my proc gmap code, but the US outline is not lining up with the existing boundaries. I saw recommendations to combine the data sets, project the combined file and then separate them again to use in gmap proc, but that caused errors that prevented the US outline from showing up at all. I have pasted the first set of code that produces the map and misaligned overlay below: 

/*Creating new US data set with added var region so a US outline can be created*/
data newus;
 set maps.us (where=(state not in (2, 15)));
 original_order=_n_;
 flag=2;
 region=1;
run;
 
proc gremove data=newus out=usoutline;
 by region;
 id id;
run;
 
/* Defining our .shp file as a map in SAS */
proc mapimport datafile="/projects/HRR_Bdry.SHP" out=hrrmap create_id_;
 select HRR_BDRY_I HRRNUM HRRCITY; 
run;
 
data hrrmap;
 set hrrmap;
 flag=1;
run;
 
proc gproject data=hrrmap degrees eastlong out=hrrproj;
 id hrrnum;
run;
 
 
/*Create us border annotation data*/
data anno_us_border;
 set usoutline (where=(state not in (2, 15)));
 by region segment;
 flag=2;
 retain size 1.5 color 'black' xsys '2' ysys '2' when 'a';
 if first.segment or (lag(x)=. and lag(y)=.)  then function='POLY    ';
   else function='POLYCONT';
   if x and y then output;
run;
 
proc gmap data = datout.prate map=hrrproj anno=anno_us_border;;
  id hrrnum;
  choro prate / levels=4 nolegend;
run;
quit;

The code above does produce a map, but the US outline boundaries are just slightly off from the HRR boundary. I am using SAS 9.4 (TS1M5) and the MAPSGFK US map. The HRR shapefile I imported is an unprojected NAD83 file. I did try combining the unprojected HRR shapefile with the US outline file and then projecting both, but that produced errors that prevented the US outline from showing up at all. I've pasted that code below: 

/*Creating new US data set with added var region so a US outline can be created*/
data newus;
 set maps.us (where=(state not in (2, 15)));
 original_order=_n_;
 flag=2;
 region=1;
run;
 
proc gremove data=newus out=usoutline;
 by region;
 id id;
run;
 
/* Defining our .shp file as a map in SAS */
proc mapimport datafile="/projects/HRR_Bdry.SHP" out=hrrmap create_id_;
 select HRR_BDRY_I HRRNUM HRRCITY; 
run;
 
data hrrmap;
 set hrrmap;
 flag=1;
run;
 
data combined;
 set usoutline hrrmap;
 run;
 
proc gproject data=combined degrees eastlong                   
                   PARALLEL1 = 30.753722
                 PARALLEL2 = 43.173294
                 MERIDIAN = -95.841165  out=combined dupok;
 id region hrrnum;
run;
 
data us hrr;
 set combined;
 if flag=1 then output hrr;
  else if flag=2 then output us;
 run;
 
data anno_state_border;
 set us;
 by region ;
 function='POLY';
run;
 
proc gmap data = datout.racerate map=hrr anno=anno_state_border;
  id hrrnum;
  choro hisprate / levels=4 nolegend ;
run;
quit;

If anyone has any suggestions, I'd greatly appreciate it.

Thanks.

[GraphGuy]

The sas-supplied maps.us is a projected map, containing only the projected X & Y values. This is an old-old legacy map that has been shipped with SAS a long time, and the original projection parameters are 'unknown', and therefore trying to project another map the exact same way (to combine them) is not really possible. (You might come close, using trial-and-error ... but it's just not a "good practice".)

I would recommend using the newer mapsgfk.us instead, which contains unprojected lat/long values.

data newus; set mapsgfk.us_states (where=(statecode not in ('AK' 'HI')) drop = x y);
run;

Then if the imported shapefile has unprojected lat/long values, you can project that and your newus any way you want 🙂

Or, if your shapefile has projected values, and you know what proj.4 projection was used, you can use that same proj.4 projection to project newus the same, or you can use that proj.4 to 'unproject' your shapefile coordinates to get the unprojected lat/long (for example, using gproject parameters such as ... from="EPSG:2264" to="EPSG:4326") . And then once you have both maps in unprojected lat/long, you can project both maps with whatever projection you want.

My rule-of-thumb is start by getting everything in unprojected lat/long 🙂