Maps Online >Resources SAS logo
Feedback   

Sample Data
Sample Programs
Tools
Useful Links

 /*******************************************************************\
| Copyright (C) 1998 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.                              |
\*******************************************************************/


/* Re-projecting US, USCITY and USCENTER Datasets */

/*******************************************************************
|                                                                 |
| This program was used to reproject US, USCITY and USCENTER data |
| sets from the STATES data set and the LONG/LAT variables in     |
| the USCITY and USCENTERS data sets; prior to release 8.2.       |
|                                                                 |
| It was necessary to re-project these datasets with the purchase |
| of new city data.                                               |
| Please be aware that this program cannot be used ASIS. Any new  |
| information being used should first be processed.               |
|                                                                 |
| This program uses datasets from SAS MAPS library and new city   |
| information from RAWCITY dataset.In the end, this job produces  |
| three work data sets: US,USCITY, USCENTER and graphics output.  |
| In the mean time the following naming convention is used for    |
| work ds:                                                        |
|                                                                 |
| Prefix - Indicates geographic area                              |
|     US   - All 50 states                                        |
|     CONT - 48 states                                            |
|     AK   - Alaska                                               |
|     HI   - Hawaii                                               |
| Middle - Projection / reduction                                 |
|     U    - Unprojected                                          |
|     P    - Projected                                            |
|     R    - Projected and reduced                                |
|   blank  - Projected low density                                |
| Suffix - Indicates data's type(s)                               |
|     MAP  - Only map points                                      |
|     CTY  - Cities                                               |
|     CTR  - Centers                                              |
|     ALL  - Map, cities, and centers                             |
|                                                                 |
*******************************************************************/

/*******************************************************************
*    Combine STATES CITIES and CENTERS                            *
*******************************************************************/

data USUMAP;
set MAPS.STATES;
FLAG = 1;
if STATE = 02 then do;
if SEGMENT = 1 then output;
end;
else if STATE = 25 | STATE = 36 | STATE = 53 then do;
if SEGMENT = 1 then output;
end;
else output;
run;

/*********** Process new city information ***************************/
data USUCTY;
set RAW.RAWCITY;
D2R = atan(1) / 45;
retain D2R;
if LATD =. | LONGD = . then delete;
else do;
if LATS = . then LATS = 0;
if LONGS = . then LONGS = 0;
LONG = LONGD + LONGM/60 + LONGS/3600;
LAT  = LATD + LATM/60 + LATS/3600;
X = LONG * D2R;
Y = LAT  * D2R;
SEGMENT = 99;
FLAG = 2;
end;
drop ST LATD LATM LATS LONGD LONGM LONGS D2R;
run;

data USUCTR;
set MAPS.USCENTER;
D2R = atan(1) / 45;
retain D2R;
X = LONG * D2R;
Y = LAT  * D2R;
SEGMENT = 98;
FLAG = 3;
drop D2R;
run;

data USUALL;
set USUMAP USUCTY USUCTR;
run;

/*********************** Seperate Alaska and Hawaii *****************/

data CONTUALL AKUALL HIUALL;
set USUALL;
if STATE = 15 then output HIUALL;
else if STATE = 2 then  output AKUALL;
else output CONTUALL;
run;

/********************* Project **************************************/

proc gproject data=CONTUALL out=CONTPALL;
id STATE;
run;

proc gproject data=AKUALL out=AKPALL;
id STATE;
run;

proc gproject data=HIUALL out=HIPALL;
id STATE;
run;

/******************* Determine Scaling Factors **********************/

title 'Range of Alaska';
proc means noprint MIN MAX data=AKPALL;/* find range of proj. Alaska*/
var X Y;
output out=AKMEAN min=XMIN YMIN max=XMAX YMAX;
run;

title 'Range of Hawaii';
proc means noprint MIN MAX data=HIPALL; /*find range of proj.Hawaii */
var X Y;
output out=HIMEAN min=XMIN YMIN max=XMAX YMAX;
run;

title 'Range of Continental US';
proc means noprint MIN MAX data=CONTPALL;
var X Y;
output out=CONTMEAN min=CXMIN CYMIN max=CXMAX CYMAX;
run;

/********************************************************************
*            Move Alaska and Hawaii                                *
********************************************************************/

data AKPALL;
/********************************************************************/
/*              rescale Alaska to fit under California              */
/* --------------------------------------------------------         */
/*   Scale X to (CXMIN, CXMIN + 0.200 * (CXMAX - CXMIN))            */
/*   Scale Y to (CYMIN, CYMIN + 0.285 * (CYMAX - CYMIN))            */
/* --------------------------------------------------------         */
/********************************************************************/
set AKPALL;
if _n_ = 1 then do;
set CONTMEAN;
FACT1 =  CXMIN + 0.200 * (CXMAX - CXMIN);
FACT2 =  CYMIN + 0.285 * (CYMAX - CYMIN);
retain FACT1 FACT2;
set AKMEAN;
AX = (FACT1 - CXMIN) / (XMAX - XMIN);
AY = (FACT2 - CYMIN) / (YMAX - YMIN);
BX = CXMIN - AX * XMIN;
BY = CYMIN - AY * YMIN;
retain AX BX AY BY;
end;
X = X * AX + BX;
Y = Y * AY + BY;
drop FACT1-FACT2 BX BY AX AY XMAX XMIN YMAX YMIN;
drop CXMAX CXMIN CYMAX CYMIN;
run;

data HIPALL;

/*****************************************************************/
/*           rescale Hawaii to fit under New Mexico              */
/* -----------------------------------------------------         */
/*   Scale X to (CXMIN + 0.200 * (CXMAX - CXMIN),                */
/*               CXMIN + 0.300 * (CXMAX - CXMIN))                */
/*   Scale Y to (CYMIN, CYMIN + 0.143 * (CYMAX - CYMIN))         */
/* -----------------------------------------------------         */
/*****************************************************************/
set HIPALL;
if _n_ = 1 then do;
set CONTMEAN;
FACT1 =  CXMIN + 0.200 * (CXMAX - CXMIN);
FACT2 =  CXMIN + 0.300 * (CXMAX - CXMIN);
FACT3 =  CYMIN + 0.143 * (CYMAX - CYMIN) ;
retain FACT1 FACT2 FACT3;
set HIMEAN;
AX = (FACT2 - FACT1) / (XMAX - XMIN);
AY = (FACT3 - CYMIN) / (YMAX - YMIN);
BX = FACT1 - AX * XMIN;
BY = CYMIN - AY * YMIN;
retain AX BX AY BY;
end;
X = X * AX + BX;
Y = Y * AY + BY;
drop FACT1-FACT3 BX BY AX AY XMAX XMIN YMAX YMIN;
drop CXMAX CXMIN CYMAX CYMIN;
run;

/********************* Separate Data Sets ************************/

data AKPCTY AKPCTR AKPMAP;
set AKPALL;
if FLAG = 1 then output AKPMAP;
else if FLAG = 2 then output AKPCTY;
else output AKPCTR;
run;

data HIPCTY HIPCTR HIPMAP;
set HIPALL;
if FLAG = 1 then output HIPMAP;
else if FLAG = 2 then output HIPCTY;
else output HIPCTR;
run;

data CONTPCTY CONTPCTR CONTPMAP;
set CONTPALL;
if FLAG = 1 then output CONTPMAP;
else if FLAG = 2 then output CONTPCTY;
else output CONTPCTR;
run;

/********************* Assign Densities **************************/

title 'Reducing Maps';
proc greduce data=HIPMAP out=HIRMAP n1=114;
id STATE;
run;

proc greduce data=AKPMAP out=AKRMAP n1=109;
id STATE;
run;

proc greduce data=CONTPMAP out=CONTRMAP n1=84;
id STATE;
run;

/********************* Strip Maps to Low Density *****************/

data HIMAP;
set HIRMAP;
if DENSITY <= 1;
run;

data AKMAP;
set AKRMAP;
if DENSITY <= 1;
run;

data CONTMAP;
set CONTRMAP;
if DENSITY <= 1;
run;
/****************** Remove Alaska Crossover Points ******************/

data AKMAP;
set AKMAP;
if _n_ = 11 then delete;
if _n_ = 64 then delete;
if _n_ = 87 then delete;
run;

/********************* Recombine States ****************************/

data USMAP; /* add repostioned Alaska and Hawaii to US */
set CONTMAP AKMAP HIMAP;
run;

data USCTY; /* add repostioned Alaska and Hawaii to US */
set CONTPCTY AKPCTY HIPCTY;
run;

data USCTR;/* add repostioned Alaska and Hawaii to US */
set CONTPCTR AKPCTR HIPCTR;
run;

/********************* Sort Data Sets by STATE *********************/

proc sort data=USMAP out=USMAP;
by STATE;
run;

proc sort data=USCTY out=USCTY;
by STATE;
run;

proc sort data=USCTR out=USCTR;
by STATE descending OCEAN;
run;

/************************ Final datasets ***************************/

data US(label='Copyright(c) 1984 SAS Institute Inc.,USA');
set USMAP;
label X       = 'Projected Longitude in Radians';
label Y       = 'Projected Latitude in Radians';
label SEGMENT = 'State Segement Number';
label STATE   = 'State FIPS Code';
keep X Y SEGMENT STATE;
run;

data USCITY(label='Copyright(c) 1984 SAS Institute Inc.,USA');
set USCTY;
label X       = 'Projected Longitude in Radians';
label Y       = 'Projected Latitude in Radians';
label LONG    = 'Unprojected Longitude in Degrees';
label LAT     = 'Unprojected Latitude in Degrees';
label CITY    = 'City Name';
label CITYFIPS= 'City FIPS Code';
label STATE   = 'State FIPS Code';
label CAPITOL = 'Y or N';
label ALT     = 'Altitude in Feet';
label POP     = 'Approximate Population in Thousands';
keep X Y LONG LAT CITY STATE CAPITOL ALT POP CITYFIPS;
run;

data USCENTER(label='Copyright(c) 1984 SAS Institute Inc.,USA');
set USCTR;
label X       = 'Projected Longitude in Radians';
label Y       = 'Projected Latitude in Radians';
label LONG    = 'Unprojected Longitude in Degrees';
label LAT     = 'Unprojected Latitude in Degrees';
label CITY    = 'City Name';
label OCEAN   = 'Y/N';
label STATE   = 'State FIPS Code';
keep X Y LONG LAT CITY OCEAN STATE;
run;

title 'Contents of US Map';
proc contents data=US;
run;

title 'Contents of USCITY';
proc contents data=USCITY;
run;

title 'Contents of USCENTER';
proc contents data=USCENTER;
run;

/*******************************************************************/
/*                                                                 */
/*   Output the US map, USCENTER and USCITY data sets.             */
/*                                                                 */
/*******************************************************************/

filename GSASFILE 'us.gif';

goption dev=GIF
noprompt
rotate=LANDSCAPE
gaccess=gsasfile;

pattern v=E r=99;
footnote "As of &SYSDATE";

proc datasets dd=WORK mt=CATALOG;
delete US;
run;

ods listing close;
ods html file="us.html" (title="USA State Map")
nogtitle nogfootnote;

title f=NONE h=1 'USCENTER Data Set on US Map';
data CENTER;
length COLOR FUNCTION $ 8;
set USCENTER;
RETAIN FLAG COLOR 'black';
IF _N_ = 1
THEN FLAG = 0;
XSYS = '2';
YSYS = '2';
WHEN = 'A';
TEXT = FIPSTATE (STATE);
FUNCTION = 'LABEL   ';
SIZE=.65;
IF FLAG = 1
THEN DO;
FUNCTION = 'DRAW    ';
FLAG = 0;
END;
IF OCEAN = 'Y'
THEN DO;
POSITION = '6';
OUTPUT;
FUNCTION = 'MOVE    ';
FLAG = 1;
END;
OUTPUT;
RUN;

PROC GMAP DATA=US MAP=US gout=WORK.US;
ID STATE;
CHORO  STATE / NOLEGEND anno=CENTER name='CENTER'
des='Visual Centers of US Map';
RUN;

TITLE f=NONE h=1 'USCITY Data Set on US Map';
DATA CITY;
retain XSYS YSYS '2' WHEN 'A';
SET USCITY;
TEXT='*';
run;

proc gmap data=US map=US gout=WORK.US;
by STATE;
id STATE;
choro  STATE / nolegend anno=CITY name='CITY' des='US Cities';
run;

ods html close;
ods listing;