Here's a slightly modified version of the code, that uses lat/long values until the end, and then projects them to projected x/y values at the very end. The original SGF paper example used maps.states and deals with converting between degrees and radians, whereas I use mapsgfk.us_states and keep the coordinates in degrees until I project them.
%macro make_anno(annodata,inputdata,
rwidth,rdepth,lcolor,fcolor,solid,center,tipvar);
data &annodata;
length function $8 color $ 12 style $20 ;
retain xsys ysys '2' hsys '3' when 'a' size .7;
/* default grid block size variables, if not passed in (in degrees) */
width=1; depth=1;
%if (&inputdata ^= ) %then %do;
set &inputdata;
%end;
%if (&solid = 'Y') %then %do;
style='solid';
%end;
%else %do;
style='empty';
%end;
line=1;
/* optional arguments */
%if (&rwidth ^= ) %then %do;
width=&rwidth;
%end;
%if (&rdepth ^= ) %then %do;
depth=&rdepth;
%end;
/* optional tool tips */
%if (&tipvar ^= ) %then %do;
html= 'title=' || quote("&tipvar" ||'=' || trim(left(&tipvar)) );
%end;
start_long =long;
start_lat= lat;
/*** create the rectangle ***/
/* coordinates of the rectangle */
function='POLY'; color=&fcolor;
long=start_long-width; lat=start_lat+depth; output; /* start position */
function='POLYCONT';
long=start_long+width; lat=start_lat+depth; output;
long=start_long+width; lat=start_lat-depth; output;
long=start_long-width; lat=start_lat-depth; output;
run;
%mend;
/* get min & max values, to use for scaling colors */
%macro get_minmax(inds,varname);
%global minval;
%global maxval;
data tempmm; set &inds END=lastob;
retain maxval -999999999;
retain minval 999999999;
if &varname <= minval then
minval=&varname;
else if &varname >= maxval then
maxval=&varname;
if (lastob = 1) then do;
output;
call symput("minval",minval);
call symput("maxval",maxval);
end;
run;
%mend;
%macro draw_map(mapdata,respdata, annodata, id, choro, name, dev, mapcolor, htmlvar);
GOPTIONS reset=all;
GOPTIONS CBACK=white;
GOPTIONS DEVICE=&dev xpixels=700 ypixels=500;
ODS LISTING CLOSE;
%if ((&mapcolor ne ) and (&mapcolor = "N")) %then %do;
pattern v=me c=CX000000 r=100;
%end;
%else %do;
pattern v=s c=CXE9E8DC r=100;
%end;
%if (&htmlvar ^= ) %then %do;
%let htmlv=html=&htmlvar;
%end;
%else %do;
%let htmlv=;
%end;
goptions border;
proc gmap map=&mapdata anno=&annodata data=&respdata all;
id &id;
choro &choro /nolegend &htmlv;
run;
quit;
%mend;
data my_data;
infile datalines dsd truncover;
input orientation:$10. Target:$15. k:32. station:$8. lat:8.1 long:8.1 xpos1:32. ypos1:32. ksum:32. ksumk:32.;
format lat 8.1 long 8.1;
datalines;
45/35,360/-30/100,20,USL000,34.5,-110.7,360,-30,3,23
45/35,390/-30/100,18,USL000,34.5,-110.4,390,-30,3,21
45/35,420/-60/100,20,USL000,34.2,-110.1,420,-60,4,24
45/35,450/-60/100,20,USL000,34.2,-109.7,450,-60,4,24
45/35,450/-30/100,19,USL000,34.5,-109.7,450,-30,3,22
45/35,480/-90/100,19,USL000,34.0,-109.4,480,-90,4,23
45/35,480/-60/100,19,USL000,34.2,-109.4,480,-60,4,23
45/35,480/-30/100,19,USL000,34.5,-109.4,480,-30,3,22
45/35,510/-90/100,19,USL000,34.0,-109.1,510,-90,4,23
45/35,540/-120/100,17,USL000,33.7,-108.8,540,-120,4,21
;;;;
run;
%let oildata=my_data;
%let predvar=ksumk;
%get_minmax(&oildata,&predvar);
/* set the color for each data point */
data &oildata; set &oildata;
length hx $2;
/*make negative numbers green to black */
/*the more negative the number, the greener */
if (&predvar < 0) then do;
cnum=(&predvar * (255/&minval)); /*scale colors */
/*Make sure we stay within the threshold values*/
if (cnum < 0) then cnum=0;
else if (cnum > 255) then cnum=255;
hx=put(cnum,hex2.);
/* Transparency for 9.3 and up*/
color="A00"||trim(left(hx))||"0099";
/* Pre-9.3 cannot use transparency ****
color="CX00"||trim(left(hx))||"00";
***/
end;
/* Make Positive values black to Red */
/* Larger numbers are more red */
else do;
cnum=(&predvar * (255/ &maxval));
/*Make sure we stay within the threshold values*/
if (cnum <= 0) then cnum=0;
else if (cnum > 255) then cnum=255;
hx=put(cnum,hex2.);
/* Alpha-transparent colors */
color="A"||trim(left(hx))||"000099";
end;
run;
%make_anno(annogrid,&oildata,.15,.15,color,color,'Y','N',&predvar);
run;
data mystates;
set mapsgfk.us_states (where=(statecode in ('TX' 'OK' 'KS' 'CO' 'NM' 'AZ' 'UT' 'NV' 'CA')));
run;
/* project the map and the annotate data using the same parameters */
proc gproject data=mystates out=mystates latlong eastlong degrees project=hammer parmout=parms;
id id;
run;
proc gproject data=annogrid out=annogrid latlong eastlong degrees parmin=parms parmentry=mystates;
id;
run;
%draw_map(mystates, mystates, annogrid,state,state,oil,png,"N");
run;
... View more
