In ArcGIS, there is a built in command to find all the pixels within a distance of feature (point), and then report back how many pixels in each category.

I find that many people need to do this, to get things like percentages of landcover around study sites, etc. Well, I set out to do that same thing in manifold. I originally used the SQL approach, which has been posted here on the forum before (using contains or touches with new points for the surface cells). This methods works well, but only for small rasters.

I was given a raster with over 1 billion pixels, and the sql approach was waaayyy to slow. So i turned to UI scripting. Using transfer selection (with buffers made on the fly from the points) was much faster, but my script still had a bottleneck; after doing the transfer selection, i used a query to find the selected pixels in the raster to move out to a table for summarizing. This process took over a minute for each point, which was just too long.

Then i realized i could use copy/paste to get out those pixels, and then using sql to get the heights would be snappy. I've loaded up the map to my S3 account, you can access it from here

http://manifold.s3.amazonaws.com/TabulateAreas.map

or see just the script at the bottom

Inside is a raster with 62 million pixels, a pretty good test sized surface for testing. The process of selecting the pixels and copying them over in a summary table (which gets transform/pivot at the end) now only takes about 20 seconds a point on this test surface.

Apparently the arcgis script could do 400 point summaries in about 1 hour. My estimates are that the manifold script could now do that many in about 2 hours. This was down from about 8-10 hours before!

Oh, and make sure to read the comments at the top of script that say to turn off the autorefresh view on the map, this lets the script go faster, because just updating the view on the raster takes forever-

/**

* Analyze Pixel Selections

* This script selects pixels with a buffer of a point

* and then summarizes the numbers of pixels in each value category (pivot query)

*

* Note: Uses UI scripting, so the map component will be automatically opened

* Turn autorefresh view off, so the script doesn't have to wait for visual changes to complete

*

* TODO - should the preSummaryTable also be removed or cleared upon script completion?

* @author chrismarx

* @version 11.19.08

*/

function Main() {

//helper functions and variables

var alert = function(msg){Application.MessageBox(msg, "Alert"); };

var prompt = function(msg){var input = Application.InputBox(msg,"Prompt",""); return input; };

var app = Application;

var doc = Document;

var compset = doc.ComponentSet;

//specify components and options

//hardcoded the names just for this test script

var drawingName = "Fl_test_sites"; //prompt("Please enter name of drawing");

var mapName = "Map"; //prompt("Please enter name of map");

var surfaceName = "Fl_mosaic_utm"; //prompt("Please enter name of surface");

var surface = compset.item(surfaceName);

surface.selectNone(); //ensure there is no selection already in raster

var bufferSize = .01;

//ui scripting

var ui = Application.UserInterface, dialog, controls;

//check that the help menu is closed

if((function isHelpMenuClosed(ui){

try{

var dialog = ui.modalDialog;

alert("UI scripts cannot run while the Manifold User Manual is open. Please close and try again.");

return true;

} catch(e){ return false; }

})(ui)){

return;

}

//create temp query object

var query = doc.newQuery("Temp Query");

//make sure object is selected

var map = compset.item(mapName);

map.open();

//add temporary buffer layer (make sure points columns are copy/copy

var drawing = compset.item(drawingName);

var coordSys = drawing.coordinateSystem;

var columns = drawing.ownedTable.columnSet;

//get user colunns and intrinsic columns in separate arrays

var columnNames = function(){

var allColumnsNames = [], allColumnsTypes = [], userColumnsNames = [], userColumnsTypes = [];

for(var i=0;i<columns.count;i++){

var col = columns.item(i);

if(col.name.indexOf("(I)") !== -1){

allColumnsNames.push(col.name);

allColumnsTypes.push(col.type);

} else {

userColumnsNames.push(col.name);

userColumnsTypes.push(col.type);

}

}

return {allNames:allColumnsNames,allTypes:allColumnsTypes,userNames:userColumnsNames,userTypes:userColumnsTypes};

}();

//create temp buffer drawing

var tempBufferDrawing = doc.newDrawing("TempBufferDrawing",coordSys);

//add user created columns

var tempBufferColumnset = tempBufferDrawing.ownedTable.columnSet;

for(var i=0;i<columnNames.userNames.length;i++){

var col = tempBufferColumnset.newColumn();

col.name = columnNames.userNames[i];

col.type = columnNames.userTypes[i];

tempBufferColumnset.add(col);

}

//run create buffer query

query.text = "Insert into [" + tempBufferDrawing.name + "] ([" + columnNames.userNames.join("],[") + "],[geom (i)]) " +

"select [" + columnNames.userNames.join("],[") + "], buffer(p.[geom (i)]," + bufferSize + ") from " + drawingName + " p";

query.run();

//add the tempbuffer to the map (the map component should have the surface in it)

map.layerset.add(doc.newLayer(tempBufferDrawing));

//make sure raster is the selected layer (to access raster ui commands)

var layersControlSet = app.UserInterface.Panes("Layers").ControlSet;

layersControlSet.item("ListViewLayers").Text = surfaceName;

layersControlSet.item("ViewLayersSwitchTo").Push();

//create a summary table

//TODO check if this table is already made?

var preSummaryTable = doc.newTable("PreSummaryTable");

var preSummaryTableColumnSet = preSummaryTable.columnSet;

var colID = preSummaryTableColumnSet.newColumn();

colID.name = "ID";

colID.type = ColumnTypeFloat32;

var colHeight = preSummaryTableColumnSet.newColumn();

colHeight.name = "Height";

colHeight.type = ColumnTypeFloat32;

//loop through each buffer(records), make selection = true, transfer selection, run summary query, clear selection, repeat

var tempBufferRecords = tempBufferDrawing.ownedTable.recordset;

for(var i=0;i<tempBufferRecords.count;i++){

tempBufferDrawing.selectNone(); //clear selection

var recordID = tempBufferRecords.item(i).data("ID");

query.text = "update [" + tempBufferDrawing.name + "] set [selection (i)] = true where id = " + recordID

query.run();

//transfer selection

ui.InvokeCommand("MapTransferSelection")

dialog = ui.modalDialog;

controls = dialog.controlSet;

controls.item("ComboBoxModify").text = surfaceName;

controls.item("MapTransferSelectionSelectNone").push();

controls.item("ListViewComponents").itemSet.item(tempBufferDrawing.name).checked = true;

dialog.accept();

//use copy paste to get the records, rather than slow sql

try{

//there is no way to target last added component, so check for naming conflicts

var newSurface = compset.item(surfaceName + " 2");

alert("There is already a duplicate of this surface ending in 2. Please remove and restart");

compset.remove(query);

compset.remove(tempBufferDrawing);

return;

} catch(e){}

//if there are no naming conflicts, then continue...

surface.copy(true);

doc.paste();

//TODO check if the selection has any rows

query.text = "select " + recordID + " as id, [height (i)] from [" + surfaceName + " 2];";

query.table.copy();

preSummaryTable.paste(false);

//delete the temp surface selection copy

compset.remove(compset.item(surfaceName + " 2"));

//deselect all pixels

surface.selectNone();

}

//create the pivot query and open up the result table

try{

var pivotQuery = compset.item('PivotQuery'); //do we already have the component

} catch(e){

var pivotQuery = doc.newQuery("PivotQuery");

}

pivotQuery.text = "TRANSFORM count(*) select id from [" + preSummaryTable.name + "] GROUP BY id pivot [height (i)]";

pivotQuery.table.open();

//delete temp components

compset.remove(query);

compset.remove(tempBufferDrawing);

}