US county-level map

[AncaTilea]

Hi.

I am trying to create a county level map for the US. I have certain prevalence that I need to show. The rotation of the map seems ever so slightly off and I am not sure what I am missing:

*-------county level map------------------------------------;

*without the GPROJECT I get a weird looking inverse map;

    proc gmap data = maps.county map =  maps.county all;

    where state not in (2, 15);

    id state county;

    choro county/coutline=black discrete nolegend statistic = mean;

    run;quit;

    *-------------------inverse the X,Y coordinates for counties;

    proc gproject data = maps.county out = counties;

id county;

    run;

    goptions reset = all;

    proc gmap data = counties map = counties all;

    where state not in (2, 15);

    id state county;

    choro county/coutline=black discrete nolegend statistic = mean;

    run;quit;

It is somewhat better...but not really.

Any ideas?

Thank you.

Sincerely,

Anca.

[Fugue]

The COUNTY and COUNTIES data sets are unprojected (see https://support.sas.com/documentation/cdl/en/graphref/63022/HTML/default/viewer.htm#gproject-concept...).

Default GPROJECT options should produce an adequate looking map (see https://support.sas.com/documentation/cdl/en/graphref/63022/HTML/default/viewer.htm#gr23n01-ex.htm)?

[AncaTilea]

I know, but I can't make it work.

What is the "default" of GPROJECT, and as you can see - if you ran the code, it doesn't give a straight map.

(I removed the ID county from PROC GPROJECT with the same result - crooked graph)

Thanks.

[Fugue]

Clearly, I need to study these PROCS a bit more to understand the options better. But I adapted the SAS example code and did a quick and dirty, which generated some warnings messages but produced adequate looking maps (final map outputs).

goptions reset=all border;

data us48; /* unprojected map values */

   set maps.counties;

   if state ne 2 and state ne 15 and state ne 72;

run;

title "United States Map";

  pattern value=mempty color=blue;

proc gmap map=us48 data=us48 all density=4;

   id state county;

   choro state county / nolegend levels=1;

run;

proc gproject data=us48

              out=us48proj; /* projected map data */

   id state county ;

run;

proc gmap map=us48proj

          data=us48proj all density=4;

   id state county ;

   choro state county / nolegend levels=1;

run;

quit;

[AncaTilea]

Thank you.

I still can't get my code to run.

If I used maps.counties instead of map.county I get a 90degrees angle map.

It must be the data I am feeding in that messes up my coordinates.

I appreciate your help.

Anca.

[Fugue · Accepted solution]

Perhaps your county/counties data has been previously projected?

Anyway, I played with the code a bit more and cleaned/simplified it. This works fine in SAS 9.3. I've using the MAPS.COUNTY dataset provided by SAS.

goptions reset=all border;

data us48;

   set maps.counties;

   if state ne 2 and state ne 15 and state ne 72;

run;

title "United States Map";

  pattern value=mempty color=blue;

proc gproject data=us48

              out=us48proj;

   id state ;

run;

proc gmap map=us48proj

          data=us48proj all density=4;

   id county ;

   choro state / nolegend levels=1;

run;

quit;

[Fugue]

Actually, I tested using both data sets (COUNTY, COUNTIES). Both worked fine.

[AncaTilea]

Right.

I think my problem - which I still haven't solved, is with the data set I am merging with the county level data from SAS.

I am not trying to just generically map the counties of the US.

I have county level data and I need to map that data. To do so, I am merging it with the SAS data (maps.county or maps.counties) and then PROC GMAP.

If I run any of the codes posted in the answers here I too get the right US map.

So, my problem is in my data.

Cheers.

[overmar]

So I ammended the example from SAS/GRAPH(R) 9.2: Reference, Second Edition to come up with this code, and it creates the map with the counties and the right colors, I don't really understand how this is projected as I am used to using ArcGIS to project maps and then using the ESRI Bridge from SAS to Arc, but the code does work.

goptions reset=all border;

data us48;

   set maps.county;

   if state ne 2 and state ne 15 and state ne 72;

run;

title "United States Map";

pattern value=mempty color=blue;

proc gmap map=us48 data=us48 all density=4;

   id county;

   choro state / nolegend levels=1;

run;

proc gproject data=us48

              out=us48proj;

   id county;

run;

proc gmap map=us48proj

          data=us48proj all density=4;

   id county;

   choro county/coutline=black discrete nolegend statistic = mean;

run;

quit;

[AncaTilea]

Thanks.

It is the same as 's answer. Having this cyber-discussion has made me realize that the problem is in my data.

Cheers.

Anca.

[GraphGuy]

Fugue is definitely on the right track ... here is a little more info about "why" your first gproject attempt looked crooked...

In your first attempt, you were projecting the entire US (including Alaska), and gproject was basically trying to make the map look good/centered (from Maine, to Alaska) - this is very tough with Alaska in there, trying to cover an area that spans many-many degrees of longitude.  And then when you were removing Alaska (after-the-fact) the resulting map would have looked skewed to the right ... except that gmap automatically re-centers the map... therefore instead of looking skewed to the right, it looked "crooked".  Try plotting the map without removing Alaska, and you'll see that the continental US looks "crooked" there too, and that crookedness becomes more visually obvious when Alaska isn't there to visually balance it out.

So, as Fugue recommended, subset the map first, before gprojecting. I might clean up/simplify the code even further like this:

proc gproject data=maps.county (where=(state not in (2,15))) out=counties;
id county;
run;

proc gmap data=counties map=counties all;
id state county;
choro county/coutline=black discrete nolegend statistic=mean;
run;

[AncaTilea]

Oh, how wonderful.

I completely overlooked 's subsetting part. I kept reading it as the standard PROC GPROJECT.

I am so glad the mystery of my crooked map is solved.

This is going to be a good week!

Thank you !!!

Anca.