Heatmap with geographic features?

[genemroz]

I'm using Heatmap to display the distribution of meteor observations over a range of latitudes and longitudes. I'd like to include some geographic markers on the map to show the locations of a few major towns and highways to help the viewer view the distribution of meteor observations within a geographic context. How can I do that?

Thanks,

Gene

[DanH_sas]

What version of SAS are you using? If you are not sure, run this statement and post the response from your log:

%put &sysvlong;

I ask this question because it could greatly impact the answer to your question.

Thanks!

Dan

[genemroz]

Thanks for the prompt response!  I hope this is the right answer:

ISAS Studio University Edition.

Release: 3.8 (Basic Edition)

SAS release 9.04.01M6P11072018

Thanks,

Gene

[DanH_sas]

Good, you should be able to PROC SGMAP to get what you want. Here is an example (with non-sensical data :-)) that should produce what you want:

proc sgmap plotdata=sashelp.class;
openstreetmap;
scatter x=weight y=height / markerattrs=(symbol=circlefilled) transparency=0.5;
run;

I added the filled circle shape and transparency to help see concentrations of sighting at/near the same lat/long. If you update to the latest University edition, there are a number of new enhancements to PROC SGMAP, including COLORRESPONSE support for SCATTER and BUBBLE plots. You can use that option to help simulate your heatmap . let me know if you have any more questions about it.

Thanks!

Dan

[Reeza]

Example Two in this paper has both markers and a filled in choro map which may be what you're after.

https://www.sas.com/content/dam/SAS/support/en/sas-global-forum-proceedings/2019/3314-2019.pdf

The documentation contains several more examples:
https://documentation.sas.com/?cdcId=pgmsascdc&cdcVersion=9.4_3.5&docsetId=grmapref&docsetTarget=n18...

[GraphGuy]

Can you post an example that's similar to the map you'd like to create?

[genemroz]

Robert,

Thanks for responding.  I've looked at many of the excellent graphs on your excellent blog site and haven't quite found one that looks like what I'm envisioning--which is a semi-transparent heatmap that overlays a geographical map.  Alternatively, I'd be happy if I could plot a few reference markers indicating the locations of cities onto the Heatmap.  I'll work on preparing an example of the heatmap that displays the observation count and an example of a scatter map that displays the locations of the observations.  I'll post them here today when I have them prepared.  The desire is to merge those together somehow.  

Thanks,

Gene

[genemroz]

Robert,

Attached is a Heatmap that shows the count of meteor observations (obs) by a network of cameras over New Mexico for the month of November 2020.  The intersection of the reference lines denotes the location of Albuquerque (ABQ).  I'd like to show the locations of other cities, preferably as some labelled icon, within the latitude/longitude range of the map.  I would also like to show the location of the  two major interstate highways that bisect the area N/S (I-25) and E/W(I-40) as that would help the viewer see the coverage of the camera network in a geographic context.  I can plot the  lat/long locations of the heatmap elements onto a esri or openstreet map but don't know how to get the "obs" count displayed on that map.

Hope this helps and thanks again for your offer to help.

Gene

[DanH_sas]

Can you use the DATALABEL option to display the "obs" count?

[GraphGuy · Accepted solution]

This is not 'exactly' what you asked for, but I think it will be a good starting point...

I put together this little example that might help. Although sgmap doesn't allow you to overlay a 'heatmap', sgplot (which supports heatmaps) does let you plot geographical/map polygons, and other geographically visual reference points.

In this example, I use sgplot to draw county polygons for New Mexico, and overlay some city markers, and then overlay a heatmap (using fake/random data for the demo). If you have lat/long points along the roads, you could also annotate or overlay them as lines.

data my_map; set mapsgfk.us_counties (where=(density<=3 and statecode='NM'));
length id_plus_segment $50;
id_plus_segment=trim(left(id))||'_'||trim(left(segment));
run;

/* generate some random/fake data, for the demo */
data grid_data;
do grid_lat=32.5 to 35 by .5;
 do grid_long=-108.5 to -106 by .5;
  count=round(ranuni(123)*5);
  output;
  end;
 end;
run;

data city_data; set mapsgfk.uscity (where=(statecode='NM' and pop_type^='under10k'));
city_long=long;
city_lat=lat;
if capital='Y' then citylabel=city;
run;

data all_data; set my_map grid_data city_data;
run;

title h=12pt "SGPlot: 'Heatmap' grid overlayed on map polygons";
proc sgplot data=all_data noborder noautolegend;
polygon x=long y=lat id=id_plus_segment / tip=none
   nofill outline lineattrs=(color=black);
scatter x=city_long y=city_lat / markerattrs=(color=blue)
   datalabel=citylabel datalabelpos=top datalabelattrs=(color=blue);
heatmapparm x=grid_long y=grid_lat colorresponse=count /
   outline outlineattrs=(color=graycc) transparency=.2
   colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);
text x=grid_long y=grid_lat text=count;
gradlegend;
xaxis display=(nolabels);
yaxis display=(nolabels);
run;

[genemroz]

Robert,

Thanks so much for taking the time to work out this solution-it looks very promising.  I tried running your code but SAS University Edition kicked back an error message:  "The SAS system library MAPSGFK may not be reassigned".  I don't know what that means nor how to fix it.  Any ideas?  Thanks!

Gene

[GraphGuy]

Sound like you probably don't have the map polygons (which are shipped with the SAS/Graph product, which I don't think is included with the free University Edition).

You could try "proc print data=mapsgfk.us_counties; run;" to see if  you have the US county map dataset.

[genemroz]

Hmmm, so close yet so far.  You are correct that the free SAS University Edition does not include SAS/Graph.  I need a Plan B.

Gene

[genemroz]

Answering my own post:  Plan B is to sign up for SAS OnDemand for Academics which includes SAS/Graph.

I tried Robert Allison's code there and it ran fine.  I'll have to do some tweaking of it to fit my particular circumstances but this appears to be a viable path forward so I will accept his solution.

Thanks, Robert, for your help.

Gene

[genemroz]

I modified Robert Allison's code and introduced my own data.  Instead of the HEATMAPPARM statement which is used for discrete data, I used the HEATMAP statement to bin the continuous data. See code below. The map (attached) has become distorted.  I don't know why?  Any ideas?

data my_map; set mapsgfk.us_counties (where=(density<=3 and statecode='NM'));
length id_plus_segment $50;
id_plus_segment=trim(left(id))||'_'||trim(left(segment));
run;

/* generate some random/fake data, for the demo */
/*data grid_data;
do grid_lat=31.5 to 37 by .25;
 do grid_long=-109 to -103 by .25;
  count=round(ranuni(123)*5);
  output;
  end;
 end;
run;*/

data city_data; set mapsgfk.uscity (where=(statecode='NM' and pop_type^='under10k'));
city_long=long;
city_lat=lat;
citylabel=city;
/*if capital='Y' then citylabel=city;*/
run;

data all_csrnmdata; set my_map city_data dec20.csrnm;
run;

 

title h=12pt "SGPlot: 'Heatmap' grid overlayed on map polygons";
proc sgplot data=all_csrnmdata noborder noautolegend;
polygon x=long y=lat id=id_plus_segment / tip=none
   nofill outline lineattrs=(color=black);
scatter x=city_long y=city_lat / markerattrs=(color=blue)
   datalabel=citylabel datalabelpos=right datalabelattrs=(color=blue);
heatmap x=lng1 y=lat1  / name="heatmap" nxbins=30 nybins=30
outline outlineattrs=(color=graycc) transparency=.6
   colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);
gradlegend 'HeatMap' ;
/*heatmapparm x=lng1 y=lat1 colorresponse=count /
   outline outlineattrs=(color=graycc) transparency=.6
   colormodel=(cxfeebe2 cxfbb4b9 cxf768a1 cxc51b8a cx7a0177);*/
/*text x=grid_long y=grid_lat text=count;remove code that puts count data in each heatmap element*/
gradlegend;
xaxis display=(nolabel);
yaxis display=(nolabel);
run;