Need help with some help with Geocoding and GMAP

[HarfordKaren]

I am trying to create a Harford County, Maryland map and map out out students who take the bus. I have the following code to get their geocode, but I am at a loss on how to do the GMAP go just get our county map.  Can someone help me please?

options nomlogic nomprint symbolgen nodate nonumber cleanup;

libname ora oracle dbprompt=yes path="irdb";

libname enroll "P:\Research\RESEARCH2\Internal\Ad_Hoc\data requests\Bus_Project";

libname pathout "C:\Geocode\Data";

proc import out=student_list

datafile="P:\Research\RESEARCH2\Internal\Ad_Hoc\data requests\Bus_Project\bus_students_2015_2.xlsx"

dbms=excel replace;

sheet="bus_students_2015";

run;

proc sort data=student_list;

by pidm;

run;

data student_list (keep=pidm address city state zip2);

set student_list;

by pidm;

if street_address_1 ~=' ';

address = street_address_1;

zip2=input(substr(zip,1,5),5.0);

run;

data student_list;

set student_list;

by pidm;

zip=zip2;

run;

proc geocode

method=street

data=student_list

out=work.geocoded

lookupstreet=pathout.usm

attributevar=(TRACT);

run;

[MikeZdeb]

Hi.  Here's code for a map of Harford County, MD (FIPS codes for state and county used to select that county).  I added some stars on the map at the center of each zip located in that county according to the data set SASHELP.ZIPCODE.  I figured that if you have geocoded your data, you might want to display the locations on the map.

One star (zip center) on the upper-left must be in two counties but is assigned to Harford.  If you look at ...

http://www.melissadata.com/lookups/CountyZip.asp?fips=24025

you'll see that several zips in that county cross county boundaries.  The code is partially based on this example ...

http://support.sas.com/kb/56/332.html

* find lat/long of zips in Harford County;
data centroids;
retain xsys ysys '2' hsys '3' function 'label' color 'yellow' size 5
       style 'wingdings' text 'ab'x when 'a';
set sashelp.zipcode (keep=state county x y);
where state eq 24 and county eq 25;
run;

* project data set centroids;
proc gproject data=centroids out=anno dupok
  parmin=mapsgfk.projparm parmentry=us_counties;
  id state;
run;

* draw a map;
proc gmap data=mapsgfk.us_counties map=mapsgfk.us_counties;
where state eq 24 and county eq 25;
id county;
choro county / nolegend statistic=first annotate=anno;
run;
quit;

More ...  rather than this ...

data student_list (keep=pidm address city state zip2);
set student_list;
by pidm;
if street_address_1 ~=' ';
address = street_address_1;
zip2=input(substr(zip,1,5),5.0);
run;

data student_list;
set student_list;
by pidm;
zip=zip2;
run;

how about (don't need those BY statements since you are not doiung anything that requires them) ...

data student_list (keep=pidm address city state zip);
set student_list (rename=(zip=temp street_address1=address));
if ^missing(address);
zip = input(temp,5.);
run;

Test ...

* zip as character, length 10;

data x;

zip = '12203-2007';

run;

* zip as numeric with only first 5 digits;;

data y;

set x (rename=(zip=temp));

zip = input(temp,5.);

put zip=;

run;

LOG ...

809  data y;
810  set x (rename=(zip=temp));
811  zip = input(temp,5.);
812  put zip=;
813  run;

zip=12203

[HarfordKaren]

I like the stars included with the map, but how do I add my data points along with the zip code stars to my map.

[MikeZdeb]

Send me a few lat/long pairs and I'l show you how to add them to the map ... msz03@albany.edu .

Do you want labels on the points?

[HarfordKaren]

I like the stars with the zip codes, but how can I add my points with the stars. I have the following code so how can I incorporate the two?

options nomlogic nomprint symbolgen nodate nonumber cleanup;

libname ora oracle dbprompt=yes path="irdb";

libname enroll "P:\Research\RESEARCH2\Internal\Ad_Hoc\data requests\Bus_Project";

libname pathout "C:\Geocode\Data";

proc import out=student_list

datafile="P:\Research\RESEARCH2\Internal\Ad_Hoc\data requests\Bus_Project\bus_students_2015_2.xlsx"

dbms=excel replace;

sheet="bus_students_2015";

run;

proc sort data=student_list;

by pidm;

run;

data student_list (keep=pidm address city state zip2);

set student_list;

by pidm;

if street_address_1 ~=' ';

address = street_address_1;

zip2=input(substr(zip,1,5),5.0);

if 21000<= zip2 <= 21099;

run;

data student_list;

set student_list;

by pidm;

zip=zip2;

run;

proc geocode

method=street

data=student_list

out=work.geocoded

lookupstreet=pathout.usm

attributevar=(TRACT);

run;

data anno;

set work.geocoded(rename=state=stname);

retain xsys ysys '2' when 'a';

function='label';

style='special';

text='J';

color='blue';

size=2;

/* convert the X and Y values from degrees to radians

to match the map data set */

x=atan(1)/45 * x;

y=atan(1)/45 * y;

/* adjust the hemisphere for the longitude to match the map data set */

x=-x;

run;

data harford;

set maps.county;

if state=24 and county=25;

run;

/* combine annotate data set with unprojected map data */

data all;

set harford anno;

run;

/* project the data */

proc gproject data=all out=allp dupok;

id state county;

run;

/* separate the annotate and map data sets */

data map anno;

set allp;

if when='a' then output anno;

else output map;

run;

/* set the map patterns to be empty */

pattern1 v=me c=black r=99;

Goptions reset=all border;

ods listing close;

ods html body="bus_map.html";

title1 "Map of Bus Riders within Harford County";

footnote1 j=r "2015 Bus Riders";

proc gmap data=map map=map;

id state county;

choro state / anno=anno nolegend;

run;

quit;

ods html close;

ods listing;

[MikeZdeb]

Hi... that SAS code you posted  uses a map from the old MAPS library.  If you use my example (new MAPSGFK library), the need to combine the annotate data set with the map, project them, separate the map and annotate data set is no longer needed.

Like I said, send me some lat/long pairs that you want on the map (or the whole data set of points) ... msz03@albany.edu. 

[MikeZdeb · Accepted solution]

Hi. Here's SAS code that uses tow annotate data sets to show the zip centroids plus your geocoded data.  There are not many comments so ask questions if there's code you do not understand. 

I have one comment ... rather than use the map data sets in the MAPS library, start uing the maps in the MAPSGFK library.  As you can see, if you use MAPSGFK.US_COUNT IES, you no longer have to go through the older method of combining annnotate data with the map/project the combined data/separate the files.  Also, you do not have to convert your annotate data lat/long in degrees to radians (HOORAY).

Here's the code and the maps are attached.

* read Harford County data ... create annotate data set ;
libname hc 'z:\temp\geocoded_data.xls';

data hc;
length color style $20;
retain xsys ysys '2' hsys '3' function 'label' color 'green' size 2
       style 'wingdings' text '6c'x when 'a' state 24;
do until (last);
   set hc.'geocoded_data$'n (keep=x y _score_) end=last;
   output;
end;
* add location of Harford CC (from Google Maps);
y = 39.5583912; x =-76.2842196;
style = 'wingdings 3';
text  = '70'x;
size  = 6;
output;
run;

libname hc clear;

* project data set hc;
proc gproject data=hc out=anno_hc dupok
  parmin=mapsgfk.projparm parmentry=us_counties;
  id state;
run;

* find lat/long of zips in Harford County ... create annotate data set;
data centroids;
retain xsys ysys '2' hsys '3' function 'label' color 'yellow' size 3
       style 'wingdings' text 'ab'x when 'a';
set sashelp.zipcode (keep=state county x y);
where state eq 24 and county eq 25;
run;

* project data set centroids;
proc gproject data=centroids out=anno_zip dupok
  parmin=mapsgfk.projparm parmentry=us_counties;
  id state;
run;

* graphics options ... output to a PNG file;
goptions reset=all gunit=pct ftext='calibri' dev=png gsfname=gout xpixels=1500 ypixels=1000;

filename gout 'z:\harford.png';

* draw a map ... two annotate data sets used;
proc gmap data=mapsgfk.us_counties map=mapsgfk.us_counties anno=anno_zip;
where state eq 24 and county eq 25;
id county;
choro county / nolegend statistic=first anno=anno_hc;
run;
quit;

* add a LEGEND to describe the points;

legend1 origin=(3,15)pct value=none shape=bar(.00001,.00001)pct mode=share
label=(h=2.5
       j=l f='wingdings 3' '70'x f='calibri' ' HARFORD CC'
       j=l f='wingdings'   'ab'x f='calibri' ' ZIP CENTROIDS'
       j=l f='wingdings'   '6c'x f='calibri' ' GEOCODED (GREEN SCORE 50+, RED SCORE <50)');
   
filename gout 'z:\harford_leg.png';

* change annotate colors based on gecoding score;

data anno_hc;
set anno_hc;
select;
   when (_score_ lt 50) color='red';
   when (_score_ ge 50) color='green';
   otherwise;
end;
run;

* draw a map ... two annotate data sets used;
proc gmap data=mapsgfk.us_counties map=mapsgfk.us_counties anno=anno_zip;
where state eq 24 and county eq 25;
id county;
choro county / statistic=first legend=legend1 anno=anno_hc;
run;
quit;