/*******************************************************************\
| 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;
|
