/*******************************************************************\
| 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. |
\*******************************************************************/
/********************************************************************
* Randomly Locate Points About Centroids *
* Example 3 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';
/********************************************************************
* Macro to randomly distribute points about centroids. *
********************************************************************/
%macro ShiftXY( DS, /* Data set with X/Y vars */
maxDistFt, /* Max distance to shift point */
mult ); /* Multiplier to adjust X/Y values */
data &DS (drop=pi twopi lat1 lon1 distFt distRad tc
latNew LonNew temp seed);
set &DS;
/**************************************************************
* Set initial random number seed and constants. *
**************************************************************/
retain seed -1 pi 3.14159 TWOpi 6.2831853;
/**************************************************************
* Compute distance to move point within allowable max. *
**************************************************************/
call ranuni(seed, temp);
distFt = temp * &maxDistFt;
/**************************************************************
* Convert distance from feet to radians. *
**************************************************************/
distRad = (DistFt / 6076.12) * (pi / (180 * 60));
/**************************************************************
* Convert X/Y from lat/lon to radians. *
**************************************************************/
lat1 = y * pi / 180 / &mult;
lon1 = x * pi / 180 / &mult;
/**************************************************************
* Set random direction to move point. *
**************************************************************/
call ranuni(seed, temp);
tc = temp * TWOpi;
/**************************************************************
* Compute new Y value. *
**************************************************************/
latNew = arsin( sin(lat1)*cos(distRad) +
cos(lat1)*sin(distRad)*cos(tc) );
/**************************************************************
* Compute new X value. *
**************************************************************/
if cos(latNew) = 0 then
lonNew = lon1; /* endpoint is a pole */
else
lonNew = mod( lon1-arsin(sin(tc)*sin(distRad) /
cos(latNew))+pi,2*pi ) - pi;
/**************************************************************
* Convert X/Y from radians back to decimal degrees. *
**************************************************************/
y = latNew * 180 / pi * &mult;
x = lonNew * 180 / pi * &mult;
run;
quit;
%mend ShiftXY;
/********************************************************************
* Create customer address data set. Repeat first observation *
* ten times to place ten points at the same location. *
********************************************************************/
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
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
123 E. Wine Street J. Cheever Loophole Mullins SC 29574 985.33
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 X = . or Y = . then
output CENTROID.MISSING; /* Customers with bad ZIP Code */
else
output CENTROID.CUSTOMERS; /* Found customer locations */
run;
/********************************************************************
* Randomly scatter geocoded points about ZIP centroid. *
********************************************************************/
%ShiftXY( CENTROID.CUSTOMERS, /* Point data to be shifted */
3000, /* Max distance to move a point */
1 ); /* Coordinate unit multiplier */
/********************************************************************
* Copy the base map to the 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 to determine the its *
* coordinate system unit multiplier. *
********************************************************************/
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 in the Output Window to find the Spatial Storage Projection *
* coordinate system and units multipler. That unit multiplier *
* (10000000) is used to adjust the customer coordinate values. *
********************************************************************/
/********************************************************************
* Adjust address coordinates to match map data. *
********************************************************************/
data CENTROID.CUSTOMERS;
set CENTROID.CUSTOMERS;
X = X * 10000000; /* Adjust X by spatial units multipler */
Y = Y * 10000000; /* Adjust Y by spatial units multipler */
run;
/********************************************************************
* Set parameters to batch import the geocoded addresses as points. *
********************************************************************/
%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 */
/********************************************************************
* Invoke the batch import program. *
********************************************************************/
dm 'af c=SASHELP.GISIMP.BATCH.SCL';
/********************************************************************
* Determine which states contain customers. The base map data *
* will be subset to display only these states. *
********************************************************************/
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;
/********************************************************************
* Use the GIS Procedure to modify the map layers and add a theme *
* to the address points using customer revenue. *
********************************************************************/
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 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. Use the magnifying glass zoom tool to *
* zoom in on the point in northeastern South Carolina. You may *
* have to zoom in more than once. As you zoom, you'll see that *
* there are multiple points there. There are the ten observations *
* that were repeated in the customer address data set. Note that *
* they are not all at the same location but are randomly scattered *
* about that ZIP centroid. See the SUGI presentation paper for *
* other information on using the map. *
********************************************************************/
dm 'gis map = CENTROID.USA.STATE';
|
