BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
HyunJee
Fluorite | Level 6

I am trying to develop a map of the state of Oklahoma zip codes that includes color coding to identify how many physician types are located in particular zip codes. I have attached a sample dataset of count of physician types and the corresponding zip codes for the state of Oklahoma. There are some zip codes that are not in Oklahoma, but I planned to exclude them by only keep zip codes where the state = Oklahoma.- I have been having trouble creating the appropriate files to create the exact map I am needing. The below code is what I have written so far, but I know that it is incorrect. I appreciate any help that can be given. Thank you.

proc sort data=zipdata2;

by zip;

run;

data temp (keep = zip count x y city state);

merge zipdata2 sashelp.zipcode;

by zip;

run;

data temp2;

set temp;

if (count ne . and x NE . and y NE .) then do;

      x=atan(1)/45 * x *-1;

      y=atan(1)/45 * y;

      output;

    end;

/* put out a message for bad/missing zipcodes */

else if (count ne . and x=.) then do; 

    put "WARNING: Zipcode " zip " wasn't located.";

    end;

run;

data all;

set maps.states (where=(40)) temp2;

run;

proc gmap data=temp2 map=temp2;

id state zip;

choro count /anno=temp2;

run;

quit;

1 ACCEPTED SOLUTION

Accepted Solutions
MikeZdeb
Rhodochrosite | Level 12

Hi ... here's some SAS code that uses an OK ZCTA shapefile and your data to produce some maps.  One is just a ZCTA outline map, the other uses your MD data.  If you've never done a map and you'd like to read a bit (for free, though that's my book that Robert cites), try this ...

The Basics of Map Creation with SAS/GRAPH

There's also a paper about importing shapefiles (though it focuses on census tracts and county subdivisions) ...

From Shapefiles to Maps Via PROC MAPIMPORT and PROC GMAP

SAS code (maps attached) ...


* shapefile from ...
http://legacy.jefferson.kctcs.edu/techcenter/GIS%20data/Country/USA/State/OK/Polygons/ ;

proc mapimport datafile="z:\zt40_d00.shp" out=mymap;
select zcta zt40_d00_i;
rename zcta=zip;
run;

* segment fix;

data mymap;
set mymap;
segment=(100*zt40_d00_i)+segment;
drop zt40_d00_i;
run;

* mymap data set contains unprojected X/Y coordinates ... project them;

proc gproject data=mymap out=zproj degrees eastlong project=robinson;
id zip;
run;

* graphics options;

goptions reset=all ftext='calibri' htext=2 gunit=pct;

* create an outline map of OK ZCTAS;

proc gmap data=zproj map=zproj;
id zip ;
choro segment / nolegend levels=1 coutline=white;
run;
quit;

* change ZIP in your data to a character variable (to match imported data);

data okmds;
set z.zipdata2 (rename=(zip=temp));
zip = put(temp,z5.);
keep zip count;
run;

* title plus white space on sides of map;
title1 h=4 'OKLAHOMA ... PHYSICIANS BY ZIP' ls=2;
title2 a=90 ls=2;
title3 a=-90 ls=2;

* some colors from ... http://colorbrewer2.org/ ;

pattern1 v=ms c=cxD0D1E6;
pattern2 v=ms c=cx74A9CF;
pattern3 v=ms c=cx3690C0;
pattern4 v=ms c=cx0570B0;

* relocate LEGEND ... into Texas panhandle;

legend1 mode=share origin=(10,35)pct across=1 label=(position=top 'PHYSICIANS') shape=bar(3,4)pct;

* some ranges for MD counts;

proc format;
value count
1 = '1'
2-3 = '2-3'
4-5 = '4-5'
6-high = '6+'
;
run;

* map your data;
proc gmap data=okmds map=zproj all;
id zip ;
choro count / discrete coutline=blue legend=legend1;
format count count.;
run;
quit;


ok_zip_outline.pngok_zip_md_counts.png

View solution in original post

9 REPLIES 9
GraphGuy
Meteorite | Level 14

sashelp.zipcode just has the centroid of the zip codes, not the border/boundary of them.  Therefore you cannot use that data as the map=.  If you want the zip code boundaries, you'll need to download them and import them into a SAS data set that can be used as a map.


http://robslink.com/SAS/democd17/zcta.htm

http://robslink.com/SAS/democd17/zcta_info.htm

If you're going to use the sashelp.zipcodes to annotate points/markers on a map, then you can do that on the SAS county maps, such as ...

http://robslink.com/SAS/democd36/ski_unblinged.htm

http://robslink.com/SAS/democd36/ski_unblinged_info.htm

HyunJee
Fluorite | Level 6

Thank you Robert.

I was unable to open the two last links you provided. The window that pops up says they cannot be found.

The first two links are helpful that you provide. Though I am unsure the code to write to produce the map I need.

I have included the new code below I am using. I should state that I have never created maps in SAS before and this is my first time trying this. I am a bit confused on how to create them and appreciate any help that can be given. Thanks.

*creating a dataset with the correct zip code variables for merging;

data zipdata2;

set zipdata;

length zcta $5.;

zcta = zip_1;

zip = zip_1;

run;

*creating a dataset of the zip code boundaries;

proc mapimport

datafile="./zt40_d00.shp"

out=mymap;

run;

*x/y is really long/lat degrees, so converted them to radians,

so I can gproject them...;

data mymap; set mymap;

  /* Convert long/lat to radians */

  x=atan(1)/45 * x;

  y=atan(1)/45 * y;

run;

*Assigned unique segment values to avoid stray lines between disjoint

3-digit zcta's;

data mymap; set mymap;

segment=(100*ZT40_D00_I)+segment;

run;

*Project the map (note the use of 'eastlong'!);

proc gproject data=mymap out=mymap dupok eastlong project=robinson;

id zcta;

run;

*merging my data with the mymap data;

proc sort data=zipdata2;

by zcta;

run;

proc sort data=mymap;

by zcta;

run;

data combined;

merge zipdata2( in=a) mymap (in=b);

by zcta;

if a and b;

run;

*combining the mymap and my data with te sashelp.zipcode data;

data sas_zipcode;

set sashelp.zipcode;

where state = 40;

length zcta $5.;

zcta = zip;

run;

data combined2;

merge combined sas_zipcode;

by zip;

run;

data mymap2;

merge mymap sas_zipcode;

by zcta;

run;

*creating a graph of the data;

proc gmap data=combined2 map=mymap2;

id state zcta;

choro count;

run;

quit;

ballardw
Super User

If you are trying to show each type of physician in each ZIP you may want to use the physician type as a BY variable producing one map per physician type.

I have done something similar and unless your final map display is going to be huge a statewide map is unlikely to display ZIP code areas within dense population areas such as Oklahoma City very well. I usually end up at a city level for maps displayed on web sites or in publications. It may also be relatively easy to find a city location SHP file in the same projection as other state features and avoid the projection completely.

HyunJee
Fluorite | Level 6

I am actually trying to show the count of physicians, the type does not matter in this instance. This is a map that is not going to be published, but used a visual tool in presentations given.

I am wondering if anyone can look over my code to see what I am missing or adding that is causing the map to not be developed. I do not have any errors in my log, but no map is produced. Thank you for your help.

GraphGuy
Meteorite | Level 14

I don't think you need to sort your map and then merge it with your zip code data.

A map is a map, and data is data ... you do not need to merge them.

You are probably at the stage where Mike Zdeb's book would be very useful to you:

http://www.amazon.com/Maps-Made-Using-Carpenters-Software/dp/1590470931

GraphGuy
Meteorite | Level 14

The link should be fixed now.

MikeZdeb
Rhodochrosite | Level 12

Hi ... here's some SAS code that uses an OK ZCTA shapefile and your data to produce some maps.  One is just a ZCTA outline map, the other uses your MD data.  If you've never done a map and you'd like to read a bit (for free, though that's my book that Robert cites), try this ...

The Basics of Map Creation with SAS/GRAPH

There's also a paper about importing shapefiles (though it focuses on census tracts and county subdivisions) ...

From Shapefiles to Maps Via PROC MAPIMPORT and PROC GMAP

SAS code (maps attached) ...


* shapefile from ...
http://legacy.jefferson.kctcs.edu/techcenter/GIS%20data/Country/USA/State/OK/Polygons/ ;

proc mapimport datafile="z:\zt40_d00.shp" out=mymap;
select zcta zt40_d00_i;
rename zcta=zip;
run;

* segment fix;

data mymap;
set mymap;
segment=(100*zt40_d00_i)+segment;
drop zt40_d00_i;
run;

* mymap data set contains unprojected X/Y coordinates ... project them;

proc gproject data=mymap out=zproj degrees eastlong project=robinson;
id zip;
run;

* graphics options;

goptions reset=all ftext='calibri' htext=2 gunit=pct;

* create an outline map of OK ZCTAS;

proc gmap data=zproj map=zproj;
id zip ;
choro segment / nolegend levels=1 coutline=white;
run;
quit;

* change ZIP in your data to a character variable (to match imported data);

data okmds;
set z.zipdata2 (rename=(zip=temp));
zip = put(temp,z5.);
keep zip count;
run;

* title plus white space on sides of map;
title1 h=4 'OKLAHOMA ... PHYSICIANS BY ZIP' ls=2;
title2 a=90 ls=2;
title3 a=-90 ls=2;

* some colors from ... http://colorbrewer2.org/ ;

pattern1 v=ms c=cxD0D1E6;
pattern2 v=ms c=cx74A9CF;
pattern3 v=ms c=cx3690C0;
pattern4 v=ms c=cx0570B0;

* relocate LEGEND ... into Texas panhandle;

legend1 mode=share origin=(10,35)pct across=1 label=(position=top 'PHYSICIANS') shape=bar(3,4)pct;

* some ranges for MD counts;

proc format;
value count
1 = '1'
2-3 = '2-3'
4-5 = '4-5'
6-high = '6+'
;
run;

* map your data;
proc gmap data=okmds map=zproj all;
id zip ;
choro count / discrete coutline=blue legend=legend1;
format count count.;
run;
quit;


ok_zip_outline.pngok_zip_md_counts.png
MikeZdeb
Rhodochrosite | Level 12

Hi ... one more thing, if you'd like to add some city locations to orient folk using the map ...

* add three cities ... use PROC GEOCODE to find lat/long;
data okcities;
retain state 'OK';
input city & $20. @@;
datalines;
Oklahoma City    Tulsa    Norman
;

proc geocode data=okcities out=okcities_xy (keep=city x y);
run;

* combine unprojected map and city locations;

data map_plus;
set mymap okcities_xy;
run;

* project combined data set;

proc gproject data=map_plus out=map_plus degrees eastlong asis project=robinson;
id zip;
run;

* separate project map from projected city locations;

data
map (where=(city is missing))
cities (where=(city is not missing) keep=city x y)
;
set map_plus;
run;

* create an annotate data set to add cities to map (transparent color ... blue, might not be best choice);

data cities;
retain xsys ysys '2' hsys '3' function 'label' color 'a0000ffaa'  style '"calibri/bo"' size 3 when 'a';
set cities ;
text = upcase(city);
run;

goptions dev=png gsfname=gout xpixels=1024 ypixels=768;

filename gout 'z:\ok.png' ;

* map your data with cities;
proc gmap data=okmds map=map all anno=cities;
id zip;
choro count / discrete coutline=blue legend=legend1;
format count count.;
run;
quit;


ok.png
poaform
Calcite | Level 5

hey mike, thanks a lot for the post and the info!

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

Click image to register for webinarClick image to register for webinar

Classroom Training Available!

Select SAS Training centers are offering in-person courses. View upcoming courses for:

View all other training opportunities.

Discussion stats
  • 9 replies
  • 1772 views
  • 0 likes
  • 5 in conversation