Together at Last: Spatial Analysis and SAS Mapping? Not quite yet.

[genemroz]

Esteemed Advisers:

I'm trying to implement an analog to the spatial mapping program that was developed and demonstrated by D. Massengill and A. Kolovos (SAS Global Forum 2012 Paper 284-2012, "Together at Last: Spatial Analysis and SAS Mapping").

I've been successful in replicating their results using their example dataset and codes. I have not succeeded when using my own dataset and their codes.

Specifically, the logfile shows multiple messages of the type:

NOTE: PROBLEM IN OBSERVATION 1 -
DATA SYSTEM REQUESTED, BUT VALUE IS NOT ON GRAPH 'X' -23707.7
DATA SYSTEM REQUESTED, BUT VALUE IS NOT ON GRAPH 'Y' 9944.658

X and Y variables are the appropriate longitude and latitude values as called for by the code. But the associated numbers in the logfile are meaningless to me.

Attached is the code that produces the result above.

I'm hoping that someone familiar with the paper by Massengill and Kolovos might take a look at what's happening here and can help me out.

Thanks in advance for any assistance,

Regards,

Gene

[ballardw · Accepted solution]

The message you are getting is typically associated with requesting a plot point, such as annotate data, that is outside the bounds of the MAP data set. You can check the bounds on your map data set by

Proc means data=<your map dataset name goes here> max min;

   var x y;

run;

Which gets for this example:

Variable	Label	Maximum	Minimum
X Y	Projected Longitude from Radians Projected Latitude from Radians	0.0351414 0.1421963	-0.3721635 -0.1938140

One suspects the map data set they used in the original example either had a different projection or missed the step to project your data from lat/long to x/y or convert to radians. You may have to watch the sign of the results when you get done.

Which will display the maximum and minimum values of the x and y variables. One strongly suspects that min of x is larger than -23707 and that the max of Y is less than 9944 in the map set you are using.

Run the same code on your annotate data set

[genemroz]

Yes indeed, failure to properly project coordinates was the source of the problem.

 Thanks for your help,

Gene

[GraphGuy]

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;