BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
genemroz
Quartz | Level 8

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

 

 

1 ACCEPTED SOLUTION

Accepted Solutions
ballardw
Super User

 

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

View solution in original post

3 REPLIES 3
ballardw
Super User

 

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
Quartz | Level 8

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

 

 Thanks for your help,

 

Gene

GraphGuy
Meteorite | Level 14

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;

map_grid.png

sas-innovate-2024.png

Available on demand!

Missed SAS Innovate Las Vegas? Watch all the action for free! View the keynotes, general sessions and 22 breakouts on demand.

 

Register now!

How to Concatenate Values

Learn how use the CAT functions in SAS to join values from multiple variables into a single value.

Find more tutorials on the SAS Users YouTube channel.

Click image to register for webinarClick image to register for webinar

Classroom Training Available!

Select SAS Training centers are offering in-person courses. View upcoming courses for:

View all other training opportunities.

Discussion stats
  • 3 replies
  • 475 views
  • 2 likes
  • 3 in conversation