topic Re: Grouped State Boundaries Not Drawn Correctly Using Maps.Counties in Graphics Programming

Grouped State Boundaries Not Drawn Correctly Using Maps.Counties

cjinsf — Wed, 12 Mar 2014 00:04:59 GMT

Hi!

Any help would be greatly appreciated or an explanation of why maps.counties does this!

For some reason when I create "divisions" of states and map the boundaries there are random lines drawn where there shouldn't be any; in the attached 140311_DivBordProb_dens_le_6.png there is a line on the state border where Apache, NM and McKinley, NM touch- it looks like it runs the length of McKinley county.

There are also random lines drawn on county lines in Mississippi and Alabama (among other states; I understand why VA has a bunch of red dots).

In the gproject step I tried playing with density; when density <= 5 the line on the AZ/ NM border disappears, but a random line is drawn on the eastern county line of Parmer, TX (compare 140311_DivBordProb_dens_le_5.png)

I really need to use maps.counties since there's a ton of other data sets that get annotated (roads, groups of counties, etc) to the map and those all work fine; just not sure why I'm having a problem at state level- thank you for any pointers!

Simplified code for Division Boundaries only (this is done in EG 5.1):

/*select states to draw*/

%let stopt = not in;

%let stlst = 'AK','HI','PR';

goptions reset=all cback=white noborder htitle=12pt htext=10pt;

GOPTIONS DEVICE=HTML;

proc gproject data=maps.counties out=base_data_1 dupok;

id state county;

where fipstate(state) &stopt. (&stlst.) /*and density <= 6*/;

run;

/* Create a data set with only STATE borders*/

proc gremove data=base_data_1 out=stborder;

by state;

id county;

run;

/* create state border annotation data*/

data anno_state_border;

set stborder;

by state segment;

retain

size 2.5 /*line thickness*/

color 'black' /*yellow border*/

xsys '2' /*Use the map coordinate system*/

ysys '2'

when 'a' ; /* Annotate after the map is drawn*/

/*For the first point in each polygon or the

first point of the interior polygon, set

FUNCTION to 'POLY'*/

if first.segment or (lag(x)=. and lag(y)=.)

then function='POLY ';

/* For other points, set FUNCTION to 'POLYCONT'*/

else function='POLYCONT';

/*Don't output points with missing values*/

if x and y then output;

run;

/*create data sets for DIVISION*/

proc sql;

create table Divisions as select

case

when state in (17,18,26,39,55) then 'East North Central'

when state in (1,21,28,47) then 'East South Central'

when state in (34,36,42) then 'Mid-Atlantic'

when state in (4,8,16,30,32,35,49,56) then 'Mountain'

when state in (9,23,25,33,44,50) then 'New England'

when state in (2,6,15,41,53) then 'Pacific'

when state in (10,11,12,13,24,37,45,51,54,72) then 'South Atlantic'

when state in (19,20,27,29,31,38,46) then 'West North Central'

when state in (5,22,40,48) then 'West South Central'

end as Division

from base_data_1

;

quit;

/* Create a data set with only DIVISION borders*/

/*need to resort by division or you will get a 'by' error*/

proc sort data=divisions out=div_border_1;

by division;

run;

proc gremove data=div_border_1 out=div_border_2;

by division;

id state county;

run;

data anno_div_border;

set div_border_2;

by division segment;

retain

size 5 /*line thickness*/

/*color 'aFFFF000a' yellow border*/

color 'red'

xsys '2' /*Use the map coordinate system*/

ysys '2'

when 'a' ; /*Annotate after the map is drawn*/

/* For the first point in each polygon or the

first point of the interior polygon, set

FUNCTION to 'POLY'*/

if first.segment or (lag(x)=. and lag(y)=.)

then function='POLY ';

/*For other points, set FUNCTION to 'POLYCONT*/

else function='POLYCONT';

/*Don't output points with missing values*/

if x and y then output;

run;

/* Combine the Annotate data sets. */

data annoall;

set anno_state_border anno_div_border ;

run;

/*draw the map*/

pattern v = me c = white r = 52;

options maxmemquery = 6M;

ods listing;

filename gout "/var/sas/projects/winoutput/140311_DivBordProb.png";

goptions

device=png

gsfname= gout

xpixels=2500 ypixels=2500;

Title1 h=12pt "Divisions";

proc gmap map=base_data_1 data=base_data_1 anno=annoall;

id state county;

choro county / coutline=black nolegend;

legend1 shape=bar(1,1) pct;

run;

quit;

Re: Grouped State Boundaries Not Drawn Correctly Using Maps.Counties

GraphGuy — Wed, 12 Mar 2014 13:12:49 GMT

Hi CJ,

Looks like you're trying some pretty intricate stuff here, and have a pretty good understanding. I've got a few suggestions that will probably make life a bit simpler!

Rather than trying to gremove the county boundaries within states, I would recommend using maps.states which has already done that for you. The counties map just has too much detail, such as the city/counties in Virginia, and a few big lakes, and maybe county boundaries on the left & right of wide rivers, etc, etc - lots of things that might leave unwanted borders.

Another thing - I notice you use squ a bit ... I believe sql can sometimes change the order of obsns in the data, which can also make 'odd' things happen in map data sets (because the order is very important).

So I've created a new version of your code, using maps.states for the state and region borders, and I replaced an sql query with a data step, and made a few other little changes. I think this should get you on the right track!

/*select states to draw*/
%let stopt = not in;
%let stlst = 'AK','HI','PR';

goptions reset=all cback=white noborder htitle=12pt htext=10pt;
GOPTIONS DEVICE=HTML;

data counties; set maps.counties (where=(fipstate(state) &stopt. (&stlst.)));
flag=1;
run;

data states; set maps.states (where=(fipstate(state) &stopt. (&stlst.)));
original_order=_n_;
flag=2;
run;

data combined; set counties states;
run;
proc gproject data=combined out=combined dupok;
id state county;
run;
data counties states; set combined;
if flag=1 then output counties;
else if flag=2 then output states;
run;

/* create state border annotation data*/
data anno_state_border; set states;
   by state segment;
   retain
      size 2.5        /*line thickness*/
      color 'black' /*yellow border*/
      xsys '2'       /*Use the map coordinate system*/
      ysys '2'
      when 'a' ;   /* Annotate after the map is drawn*/
      /*For the first point in each polygon or the
      first point of the interior polygon, set
      FUNCTION to 'POLY'*/
   if first.segment or (lag(x)=. and lag(y)=.)
      then function='POLY    ';
      /* For other points, set FUNCTION to 'POLYCONT'*/
   else function='POLYCONT';
      /*Don't output points with missing values*/
   if x and y then output;
run;

/*create data sets for DIVISION*/
data divisions; set states;
length division $100;
if state in (17,18,26,39,55) then division='East North Central';
if state in (1,21,28,47) then division='East South Central';
if state in (34,36,42) then division='Mid-Atlantic';
if state in (4,8,16,30,32,35,49,56) then division='Mountain';
if state in (9,23,25,33,44,50) then division='New England';
if state in (2,6,15,41,53) then division='Pacific';
if state in (10,11,12,13,24,37,45,51,54,72) then division='South Atlantic';
if state in (19,20,27,29,31,38,46) then division='West North Central';
if state in (5,22,40,48) then division='West South Central';
run;

proc sort data=divisions out=divisions;
by division segment original_order;
run;

/* Create a data set with only DIVISION borders*/
proc gremove data=divisions out=divisions;
by division;
id state;
run;

data anno_div_border;
   set divisions;
   by division segment;
   retain
      size 5        /*line thickness*/
      /*color 'aFFFF000a' yellow border*/
        color 'red'
      xsys '2'       /*Use the map coordinate system*/
      ysys '2'
      when 'a' ;   /*Annotate after the map is drawn*/
      /* For the first point in each polygon or the
      first point of the interior polygon, set
      FUNCTION to 'POLY'*/
   if first.segment or (lag(x)=. and lag(y)=.)
      then function='POLY    ';
      /*For other points, set FUNCTION to 'POLYCONT*/
   else function='POLYCONT';
      /*Don't output points with missing values*/
   if x and y then output;
run;

/* Combine the Annotate data sets. */
data annoall;
set anno_state_border anno_div_border ;
run;

/*draw the map*/
pattern v = me c = white r = 52;

options maxmemquery = 6M;

ods listing;

filename gout "&name..png";
goptions
device=png
gsfname= gout
xpixels=2500 ypixels=2500;

Title1 h=12pt "Divisions";

proc gmap map=counties data=counties anno=annoall;
   id state county;
   choro county / coutline=gray99 nolegend;
   legend1 shape=bar(1,1) pct;
run;
quit;

Re: Grouped State Boundaries Not Drawn Correctly Using Maps.Counties

cjinsf — Wed, 12 Mar 2014 16:48:04 GMT

Thank you! Combining the maps.states and maps.counties is great; I do this when combining roads, select locations, city, etc and then blow them out into their own anno sets- didn't think of it at the very beginning of the project using the SAS datasets (probably b/c when first starting to learn mapping figuring out which data sets used which coordinate system & gproject was a little confusing)

Yes, sql does do reordering (the code in the post is condensed down to illustrate the problem); I actually do (monotonic()) as sortord to maps.counties to be able to re-sort correctly in later steps.

Your SAS/Graph Examples is awesome- I've been playing with animating the maps and doing graphs w/ .svg; fun stuff!

Thanks again!