/*******************************************************************\
| Copyright (C) 2005 by SAS Institute Inc., Cary, NC, USA. |
| |
| SAS (R) is a registered trademark of SAS Institute Inc. |
| |
| SAS Institute does not assume responsibility for the accuracy of |
| any material presented in this file. |
\*******************************************************************/
/********************************************************************
* Locate Addresses With Bad ZIP Codes *
* Example 4 from SUGI 30 Presentation: *
* "Cheap Geocoding: SAS/GIS and Free TIGER Data" *
* *
* Note: This code is not exactly as described in the SUGI paper. *
* The code order was rearranged to go ahead and import the *
* missing address into the map with the other points. *
* A complete US map is generated instead of just the *
* southeastern states. *
********************************************************************/
/********************************************************************
* Set libref to location where your map will be created. *
********************************************************************/
libname CENTROID 'C:\Temp';
/********************************************************************
* Create a customer address data set. *
********************************************************************/
data CENTROID.CUSTOMERS (label = 'Customer Addresses');
input address $ 1-22 /* House address (not used) */
name $ 24-48 /* Customer name */
city $ 49-63 /* City name (not used) */
statecode $ 64-65 /* State abbreviation (not used) */
zip /* Customer ZIP Code (numeric) */
revenue; /* Total revenue for customer */
format revenue dollar9.0; /* Show revenue in dollars */
label revenue = 'Revenue'; /* Label to appear on map legend */
cust_ID = _n_; /* Assign customer ID number */
cards;
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
611 Jackson Street Cuthbert J. Twillie Vidalia GA 30474 2533.25
300 North Odum Street Prof. Quincy Wagstaff Pembroke NC 28372 1899.45
5575 Peachtree Pkwy. Kaspar Gutman Norcross GA 30022 33.44
1525 N. Shenandoah Ave Dr. Hugo Z. Hackenbush Front Royal VA 22630 987.02
2310 Cannon Bridge Rd Charlie Allnut Cordova SC 29039 1459.65
305 Mill Creek Drive Larson E. Whipsnade Fuquay Varina NC 27526 4566.98
1911 W. Main St. Guillermo Ugarte Richmond VA 23220 17.45
2313 Dale Ave. Dewey, Cheatem & Howe Roanoke VA 24013 989.20
7825 Old Middlesex Rd Capt. Geoffrey Spaulding Bailey NC 27807 1798.33
1935 E. Frederick St Joel Cairo Blacksburg SC 29702 2566.45
900 Stony Landing Rd Charles Blutoski Easley SC 29641 783.25
26 Ancho Street Brigid O'Shaughnessy Burlingame CA 94000 543.99
1373 London Bridge Rd Fred C. Dobbs Virginia Beach VA 23451 898.25
626 Mohican Trail Rufus T. Firefly Wilmington NC 28409 123.45
216 12th Street Otis B. Driftwood Lynchburg VA 24504 98.56
35 W. Chatham Str Elwood Blues Pittsboro NC 27312 661.13
100 W. Church Street Duke Mantee Cherryville NC 28021 4531.87
;
run;
/********************************************************************
* Reorder customer data by ZIP for merging. *
********************************************************************/
proc sort data = CENTROID.CUSTOMERS;
by ZIP;
run;
/********************************************************************
* Assign ZIP centroid coordinates to customer addresses. *
********************************************************************/
data CENTROID.CUSTOMERS (label = 'Customers Located by ZIP Centroid Geocoding')
CENTROID.MISSING (label = 'Customers Not Found by ZIP');
merge CENTROID.CUSTOMERS (in=a) /* Customer addresses */
SASHELP.ZIPCODE; /* ZIP and centroid X/Y data */
by ZIP; /* ZIP is the linkage variable */
if a; /* Keep only customer obs */
/*****************************************************************
* If this customer has a bad ZIP Code, set this obs aside. *
*****************************************************************/
if X = . or Y = . then output CENTROID.MISSING;
/*****************************************************************
* Otherwise this customer's ZIP Code was found and has valid *
* X/Y values. *
*****************************************************************/
else output CENTROID.CUSTOMERS;
run;
/********************************************************************
* At this point, you can see one unfound address is in the *
* CENTROID.MISSING data set. It has a bad ZIP Code value. *
* We will assign X/Y coordinates using its city/state values. *
********************************************************************/
/********************************************************************
* Sort coordinate data to compute mean for city/state. *
********************************************************************/
proc sort data = SASHELP.ZIPCODE /* Centroid X/Y data */
out = STATECITY; /* Ordered by state/city */
by state city;
run;
/********************************************************************
* Compute average X/Y for each city/state. *
********************************************************************/
proc means data = STATECITY noprint;
by state city;
var X Y;
output mean=
out=STATECITY (drop=_type_ _freq_
label='Average centroid X/Y for each city in SASHELP.ZIPCODE');
run;
/********************************************************************
* Reorder customer data by ZIP for merging. *
********************************************************************/
proc sort data = CENTROID.MISSING;
by state city;
run;
/********************************************************************
* Fill in state FIPS code values for customer addresses. *
********************************************************************/
data CENTROID.MISSING;
set CENTROID.MISSING;
state = stfips(statecode);
run;
/********************************************************************
* Assign city centroid to customers not found by ZIP. *
********************************************************************/
data _TEMP_ /* Found by city/state values */
CENTROID.MISSING (label='Customers Not Found');
length city $ 35; /* Equalize city var lengths */
merge CENTROID.MISSING (in=a) /* Customers with bad ZIPs */
STATECITY; /* City centroid locations */
by state city; /* Linkage variables */
if a; /* Keep only customer obs */
if X = . or Y = . then
output CENTROID.MISSING; /* Customers with bad city/state */
else
output _TEMP_; /* New found customer locations */
run;
/********************************************************************
* Append customers located by city/state to previously found. *
********************************************************************/
data CENTROID.CUSTOMERS;
set CENTROID.CUSTOMERS _TEMP_;
run;
/********************************************************************
* At this point, the address that was initially not found by ZIP *
* has been located by its city and state and appended to the *
* CENTROID.CUSTOMERS data set. We'll now add them all to a map. *
********************************************************************/
/********************************************************************
* Copy the base map to your location. *
********************************************************************/
proc gis; /* Invoke GIS Procedure */
copy MAPS.USA.STATE.GISMAP / /* Base map to copy */
destlib = CENTROID /* Copy it to this libref ... */
destcat = CENTROID.USA /* ... and to this catalog */
select = (_all_) /* Copy all map entries and data */
blank /* Clear internal path references */
replace; /* Overwrite existing map data */
run;
/********************************************************************
* Query the map's spatial catalog entry. *
********************************************************************/
proc gis;
spatial contents CENTROID.USA.USA;
run;
/********************************************************************
* The base map and the geocoded points MUST be on the same *
* coordinate system. We obtained coordinate values from the *
* SASHELP.ZIPCODE data set. The variable labels in that data set *
* show the X and Y values are degrees of longitude and latitude. *
* Look at the output from the GIS Procedure to find the map's *
* coordinate system and unit multipler. That multiplier is used *
* to adjust the customer coordinates values below. *
********************************************************************/
/********************************************************************
* Adjust address coordinates to match map data. *
********************************************************************/
data CENTROID.CUSTOMERS;
set CENTROID.CUSTOMERS;
X = X * 10000000; /* Adjust X by spatial entry units multipler */
Y = Y * 10000000; /* Adjust Y by spatial entry units multipler */
run;
/********************************************************************
* Set up macro variables required for batch importing. *
********************************************************************/
%let imp_type = GENPOINT; /* Import customers as points */
%let maplib = CENTROID; /* Map libref */
%let mapcat = USA; /* Map catalog name */
%let mapname = STATE; /* Map entry name */
%let spalib = CENTROID; /* Spatial data set libref */
%let spaname = USA; /* Spatial entry name */
%let cathow = UPDATE; /* Use existing entries */
%let spahow = APPEND; /* Add to spatial data sets */
%let nidvars = 0; /* Put points in one layer */
%let infile = CENTROID.CUSTOMERS; /* Geocoded address data set */
/********************************************************************
* Run the batch import program. *
********************************************************************/
dm 'af c=SASHELP.GISIMP.BATCH.SCL';
/********************************************************************
* Change map and layer parameters with the GIS Procedure. *
********************************************************************/
proc gis; /* Invoke Procedure */
catalog CENTROID.USA; /* Map location */
spatial USA; /* Spatial entry name */
/*****************************************************************
* Modify visual parameters of customer layer. *
*****************************************************************/
layer update CUSTOMERS / /* Layer entry name */
where = 'node = 1' /* Layer definition */
des = 'Geocoded addresses' /* Layer entry label */
default = (point = (color = cxff0000 /* Symbol color */
font = marker /* Symbol font */
character = 'V' /* Point symbol (star) */
size = 10)); /* Symbol height */
/*****************************************************************
* Add a theme using customer revenue to the layer. *
*****************************************************************/
layer update CUSTOMERS / /* Layer entry name */
theme = (create /* Make a new theme */
link = cust_rev /* Data set link name */
dataset = CENTROID.CUSTOMERS /* Theme data set */
datavar = cust_ID /* Linking variable */
themevar = Revenue /* Thematic variable */
nlevels = 3 /* Number of ranges */
point = ((level = first /* First theme range: */
color = cxff0000 /* Symbol color */
size = 6) /* Symbol height */
(level = last /* Third theme range: */
color = cx00a900 /* Symbol color */
size = 18) /* Symbol height */
blendcolor /* Color middle range */
blendsize)); /* Size middle range */
/*****************************************************************
* Modify visual parameters of the state layer. *
*****************************************************************/
layer update STATES / /* Layer entry name */
static; /* Do not show theme */
/*****************************************************************
* Modify the map visually. *
*****************************************************************/
map update STATE / /* Map entry name */
layers = (STATES CUSTOMERS) /* Layers to include */
layerson = (STATES CUSTOMERS) /* Layers to display */
cback = gray /* Background color */
legend = hideall /* Turn off legend */
action = (delete name=_all_) /* Remove old actions */
des = 'Geocoded customer map'; /* Label for entry */
/*****************************************************************
* Add a Spatial Information action to map. *
*****************************************************************/
map update STATE / /* Map entry name */
action = (create /* Make a new action */
name = Spatial /* Name of action */
type = spatial /* Action type */
when = immediate); /* Run upon selection */
/*****************************************************************
* Add labels the map. *
*****************************************************************/
maplabel create /
force /* Replace existing DS */
dataset = CENTROID.LABELS /* label data set name */
text = 'Customer Locations' /* Label text */
map = CENTROID.USA.STATE /* Map entry */
attach_to = window /* Pin label to window */
position = (top left) /* Label location */
color = cxea1010 /* Text color */
font = 'Swis721 Blk BT-12pt-Roman-Normal'; /* Text font */
maplabel create /
text = 'Geocoded with SAS' /* Label text */
map = CENTROID.USA.STATE /* Map entry */
attach_to = window /* Pin label to window */
position = (bottom right) /* Label location */
color = cxea1010 /* Text color */
font = 'Swis721 Blk BT-11pt-Roman-Normal'; /* Text font */
run;
quit;
/********************************************************************
* Open the map in SAS/GIS. See the SUGI presentation paper for *
* information on using the map. *
********************************************************************/
dm 'gis map = CENTROID.USA.STATE';
|
