/*******************************************************************\
| 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. |
\*******************************************************************/
/********************************************************************
* Centroid Level Geocoding *
* Example 1 from SUGI 30 Presentation: *
* "Cheap Geocoding: SAS/GIS and Free TIGER Data" *
********************************************************************/
/********************************************************************
* 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;
/********************************************************************
* 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 multiplier */
Y = Y * 10000000; /* Adjust Y by spatial entry units multiplier */
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';
/********************************************************************
* Determine which states contain customers. *
********************************************************************/
proc freq data = CENTROID.CUSTOMERS noprint;
tables state / out = _states_ (label = 'States containing customers');
run;
/********************************************************************
* Generate a where clause to subset the map data for only the *
* states containing customers. The clause is placed in the macro *
* variable &CoverWhere for use in later PROC GIS code. *
********************************************************************/
data _null_;
set _states_ (rename = (state = curFIPS)) end = last;
length FIPS /* FIPS codes for states */
coverWhere $ 256; /* Where-clause to subset map */
retain FIPS '';
curState = fipstate(curFIPS);
/*****************************************************************
* If the first state, initialize the where clause string. *
*****************************************************************/
if _n_ = 1 then FIPS = trim(left(curFIPS));
/*****************************************************************
* If not the first state, append it to the string. *
*****************************************************************/
else FIPS = compress(FIPS || "," || curFIPS);
if last then do;
/**************************************************************
* The last state. Assemble where clause string. The *
* 'stateL' and 'stateR' expressions subset the state chains. *
* The 'node=1' expression displays the imported addresses. *
**************************************************************/
coverWhere = 'stateL in(' || trim(FIPS) ||
')|stateR in(' || trim(FIPS) ||
')|node=1';
/**************************************************************
* If the where clause is too long, print warning and use a *
* where clause that will display the whole map. *
**************************************************************/
if length(coverWhere) > 200 then do;
put 'WARNING: Where clause maximum of 200 characters exceeded.';
put 'WARNING: Map not subset by where. Coverage will display entire country.';
coverWhere = '1';
end;
/**************************************************************
* Put where clause in macro variable for PROC GIS use. *
**************************************************************/
call symput( 'coverWhere', trim(coverWhere) );
output;
stop;
end;
run;
/********************************************************************
* Uncomment this assignment statement to display all states in *
* the map. This will replace the where clause generated above. *
********************************************************************/
* %let CoverWhere=1;
/********************************************************************
* 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 */
/*****************************************************************
* Use the &CoverWhere macro variable string to subset the *
* spatial data and display features in only the wanted states. *
*****************************************************************/
coverage replace ALL / /* Coverage entry name */
where = "&CoverWhere" /* Subsetting clause */
spatial = USA /* Spatial entry name */
des = 'Coverage for geocoded address map';
/*****************************************************************
* 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';
|
