BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
PhilWood
Obsidian | Level 7

Hi all,

I have successfully created an annotated state map from geocoded address info using the states map. However, when I try to use the county map data set, it appears that the geocoding output does not match  what information in the county map dataset. Here is a small program to illustrate the difficulty. Code doesn't get past proc gproject. It seems that the length of x and y of the map do not match the length of x and y in the geocoded output. thanks! Phil

:

libname lookup "D:\Documents\desi2023\frlanalyses\geocodedata__2006__ZIP4_Geocode_Data\data";
data zips;
input c zip;

anno_flag=1l
cards;
1 65203
2 65211
3 65201
;
proc geocode /* Invoke geocoding procedure */
method=plus4 /* Specify geocoding method */
lookup=lookup.zip4 /* Lookup data from MapsOnline */
data=zips /* Input data set to geocode */
out=geocodedzips; /* Specify name of Output data set of locations */
run;

data mocounties;set maps.counties;if state=29;

data mocounties;set mocounties geocodedzips;run;
proc gproject data=mocounties out=MOCountiesP; id state ; run;

data my_map locations; set mocountiesP;
if anno_flag=1 then output locations;
else output my_map;
run;

data anno_locations; set locations;
if c=1 then color="red";
if c=2 then color="bib";
if c=3 then color="viyg";
if c=4 then color="bippk";
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; rotate=360; size=.5;
style='psolid'; output;
length html $300;
html='title='||quote(trim(left(citystate))||'0d'x||trim(left(put(zip,z5.))));
style='pempty'; color='deb'; output;
run;
pattern1 v=e;

proc gmap data=mymap map=mymap anno=anno_locations;
id county;
choro segment / coutline=black levels=1 nolegend
coutline=gray99;run;

1 ACCEPTED SOLUTION

Accepted Solutions
ballardw
Super User

I had to seriously modify the starting geocode data set because you didn't paste it into a text box.

 

This creates a map with those 3 ZIP codes appearing.

Note that some differences are because the MY_MAP created with the geocoded data is different name than the MYMAP from your first posted code.

data WORK.GEOCODEDZIPS(label='Geocoded 15May2024');
   infile datalines dlm=',' dsd  truncover;
   length y x M_OBS 8 _matched_ $ 16 anno_flag 8 state 5 ;
   input Y X M_OBS  _MATCHED_:$16. 
         anno_flag state c zip:5.;
   label Y="Geocoded Latitude" 
         X="Geocoded Longitude"
   ;
datalines;
38.928422641,-92.36843533,.,ZIP mean,1,29,1,65203
38.948231,-92.327335,11274177,ZIP,1,29,2,65211
38.933052625,-92.29140936,.,ZIP mean,1,29,3,65201
;

data mocounties;
   set MAPSGFK.US_COUNTIES;
   if statecode="MO";
run;
data mocounties2;
   set mocounties 
   work.geocodedzips (rename=(x=long y=lat));
run;

proc gproject data=mocounties2 out=MOCountiesp latlong eastlong degree; 
   id state ; 
run;
data my_map locations; set MOCountiesp;
   if anno_flag=1 then output locations;
   else output my_map;
run;

data anno_locations; 
   set locations;
   if c=1 then color="red";
   if c=2 then color="bib";
   if c=3 then color="viyg";
   if c=4 then color="bippk";
   xsys='2'; ysys='2'; hsys='3'; when='a';
   function='pie'; rotate=360; size=.5;
   style='psolid'; 
output; length html $300; html='title='||'0d'x||put(zip,z5. -L); style='pempty'; color='deb';
output; run; proc gmap data=my_map map=my_map anno=anno_locations; id county; choro segment / coutline=black levels=1 nolegend coutline=gray99; run; quit;

Personally for the example Zip codes I would have used some from different counties. The scale of statewide map is such that it may be impossible to see one or more zip code centers.

 

Your prior code for annotate data set couldn't quite be applied as there wasn't anything adding in the city variable to use in the title.

 

View solution in original post

10 REPLIES 10
ballardw
Super User

If your MAPS.COUNTIES data set is the on provided by SAS then there is no geographic information other than the projected boundary points provided by the X,Y variable.

 

If you have the MAPSGFK.US_COUNTIES data set that has Lat, Long variables likely to have a better shot for getting geocoded data to place correctly.

PhilWood
Obsidian | Level 7

That's a great help! thanks! the only problem now is that Gproject only outputs one data point instead of the three geocoded zipcodes in the small data set. I.e.:

libname lookup "D:\Documents\desi2023\frlanalyses\geocodedata__2006__ZIP4_Geocode_Data\data";
data zips;
anno_flag=1;
state=29;
input c zip;
cards;
1 65203
2 65211
3 65201
;
proc geocode /* Invoke geocoding procedure */
method=plus4 /* Specify geocoding method */
lookup=lookup.zip4 /* Lookup data from MapsOnline */
data=zips /* Input data set to geocode */
out=geocodedzips; /* Specify name of Output data set of locations */
run;
proc contents data=MAPSGFK.US_COUNTIES;run;
data mocounties;set MAPSGFK.US_COUNTIES;if statecode="MO";run;
data mocounties;set mocounties geocodedzips;run;
proc means data=mocounties;var state;run;
proc gproject data=mocounties out=MOCountiesp latlong eastlong degree; id state ; run;
data my_map locations; set MOCountiesp;
if anno_flag=1 then output locations;
else output my_map;
run;

NOTE: There were 11766 observations read from the data set WORK.MOCOUNTIESP.
NOTE: The data set WORK.MY_MAP has 11765 observations and 18 variables.
NOTE: The data set WORK.LOCATIONS has 1 observations and 18 variables.

ballardw
Super User

Can you provide the result of your GeocodedZips data set from Proc Geocode in the form of a working data step? We don't have your lookup data set to see what the result might be.

 

Instructions here: https://communities.sas.com/t5/SAS-Communities-Library/How-to-create-a-data-step-version-of-your-dat... will show how to turn an existing SAS data set into data step code that can be pasted into a forum code box using the </> icon or attached as text to show exactly what you have and that we can test code against. Note that the code box is important as the main message window will reformat pasted text and can result in code that won't run properly.

PhilWood
Obsidian | Level 7

thanks all. OK, I *think* I did this correctly:

data WORK.GEOCODEDZIPS(label='Geocoded 15May2024');
infile datalines dsd truncover;
input Y:32. X:32. M_OBS:32. _MATCHED_:$16. anno_flag:32. state:32. c:32. zip:32.;
label Y="Geocoded Latitude" X="Geocoded Longitude";
datalines;
38.928422641 -92.36843533 . ZIP mean 1 29 1 65203
38.948231 -92.327335 11274177 ZIP 1 29 2 65211
38.933052625 -92.29140936 . ZIP mean 1 29 3 65201
;;;;

or, alternatively, if you want to see the proc print. I notice it's got X and Y, but not longitude and latitude. :

Obs Y X M_OBS _MATCHED_ anno_flag state c zip123
38.9284-92.3684.ZIP mean129165203
38.9482-92.327311274177ZIP129265211
38.9331-92.2914.ZIP mean129365201

the Gmap call is: 

proc gproject data=mocounties out=MOCountiesp latlong eastlong degree; id state ; run;

so maybe it's confused given that geocode didn't produce longitude and latitude.

rerunning the proc gproject removing with latlong option, however, doesn't fix the problem. thanks for all the help! Phil

ballardw
Super User

I would suggest renaming the variables to match the lat and long in the map data set.

BEFORE combining with the data and doing the projection.

 

Look at the selected county data a check on the sign of the long, may need to insure that matches with the values in your map set.

PhilWood
Obsidian | Level 7

I did that- same result, alas- thanks for responding, though!

ballardw
Super User

I had to seriously modify the starting geocode data set because you didn't paste it into a text box.

 

This creates a map with those 3 ZIP codes appearing.

Note that some differences are because the MY_MAP created with the geocoded data is different name than the MYMAP from your first posted code.

data WORK.GEOCODEDZIPS(label='Geocoded 15May2024');
   infile datalines dlm=',' dsd  truncover;
   length y x M_OBS 8 _matched_ $ 16 anno_flag 8 state 5 ;
   input Y X M_OBS  _MATCHED_:$16. 
         anno_flag state c zip:5.;
   label Y="Geocoded Latitude" 
         X="Geocoded Longitude"
   ;
datalines;
38.928422641,-92.36843533,.,ZIP mean,1,29,1,65203
38.948231,-92.327335,11274177,ZIP,1,29,2,65211
38.933052625,-92.29140936,.,ZIP mean,1,29,3,65201
;

data mocounties;
   set MAPSGFK.US_COUNTIES;
   if statecode="MO";
run;
data mocounties2;
   set mocounties 
   work.geocodedzips (rename=(x=long y=lat));
run;

proc gproject data=mocounties2 out=MOCountiesp latlong eastlong degree; 
   id state ; 
run;
data my_map locations; set MOCountiesp;
   if anno_flag=1 then output locations;
   else output my_map;
run;

data anno_locations; 
   set locations;
   if c=1 then color="red";
   if c=2 then color="bib";
   if c=3 then color="viyg";
   if c=4 then color="bippk";
   xsys='2'; ysys='2'; hsys='3'; when='a';
   function='pie'; rotate=360; size=.5;
   style='psolid'; 
output; length html $300; html='title='||'0d'x||put(zip,z5. -L); style='pempty'; color='deb';
output; run; proc gmap data=my_map map=my_map anno=anno_locations; id county; choro segment / coutline=black levels=1 nolegend coutline=gray99; run; quit;

Personally for the example Zip codes I would have used some from different counties. The scale of statewide map is such that it may be impossible to see one or more zip code centers.

 

Your prior code for annotate data set couldn't quite be applied as there wasn't anything adding in the city variable to use in the title.

 

PhilWood
Obsidian | Level 7

thanks so much! the actual data set I'm using is much larger, but I wanted to just give a small example of the code for problem solving. This works fine! Phil

GraphGuy
Meteorite | Level 14

Here's how to do it without combining/projecting/separating the map & data, by using the somewhat new feature of saving the projection parameters for the map and then applying them to the annotate dataset:

 

data zips;
input c zip;
cards;
1 65203
2 65211
3 65201
;

proc geocode data=zips out=geocodedzips (rename=(x=long y=lat));
run;

data mocounties; set mapsgfk.us_counties (where=(statecode='MO'));
run;

proc gproject data=mocounties out=mocounties 
latlong eastlong degrees parmout=projparm;
id county;
run;

proc gproject data=geocodedzips out=geocodedzips
latlong eastlong degrees parmin=projparm parmentry=mocounties;
id;
run;

data anno_locations; set geocodedzips;
if c=1 then color="red";
if c=2 then color="bib";
if c=3 then color="viyg";
if c=4 then color="bippk";
xsys='2'; ysys='2'; hsys='3'; when='a';
function='pie'; rotate=360; size=.5;
style='psolid'; 
output;
length html $300;
html='title='||quote(trim(left(c))||'0d'x||trim(left(put(zip,z5.))));
style='pempty'; color='deb'; 
output;
run;

pattern1 v=e;

proc gmap map=mocounties data=mocounties anno=anno_locations;
id county;
choro segment / coutline=black levels=1 nolegend
coutline=gray99;

mo_map.png

SAS Innovate 2025: Register Now

Registration is now open for SAS Innovate 2025 , our biggest and most exciting global event of the year! Join us in Orlando, FL, May 6-9.
Sign up by Dec. 31 to get the 2024 rate of just $495.
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.

SAS Training: Just a Click Away

 Ready to level-up your skills? Choose your own adventure.

Browse our catalog!

Discussion stats
  • 10 replies
  • 1231 views
  • 2 likes
  • 4 in conversation