Monday, April 21, 2008

Quantum GIS: First Impressions

As I mentioned in the welcome to this blog, part of why I want to write here is to grow professionally by trying technologies outside the four walls of my FalconView / ArcGIS office.  Quantum GIS is a piece of GIS software not much different from FalconView in its mission.  It displays vector layers, raster maps and imagery; it also includes some task-oriented geoprocessing and editing capabilities.  Being a typical engineering type, I only spent a few minutes reading the instructions, too eager to dive right in and start trying the thing.  Here are my first impressions of Quantum GIS.  Please keep in mind that these are impressions coming from an experienced GIS user with no previous exposure to QGIS, so some of my facts about QGIS may be incorrect, but I probably represent a typical GIS user who would be considering this type of solution.  Here are my impressions from my short test drive.



Right off, I discovered support for drawing GeoTIFFs and Shapefiles.  The picture above and left is an image of the Georgia Tech campus drawn in QGIS.  (I included the same image in FalconView on the right just for kicks.)  The image below shows the Shapefile.  Right off I was encouraged: QGIS is easy and intuitive to use.  Note that I opened the GeoTIFF in an Ubuntu QGIS install and the Shapefile in a Windows install.  This is a nice thing about most open source software: the authors generally do a good job of achieving platform independence.

Notice that I don't have the Shapefile and the GeoTIFF drawing together in the same window.  After a little poking around, I wasn't able to get this to work properly.  I assume my difficulty was due to differing coordinate systems between the two files and my lack of experience setting up coordinate systems in QGIS, but I'm still not sure.  Since I was going for first impressions, I gave up trying the get them to draw together after a short time.  (ArcMap and FalconView generally transform different coordinate systems to draw together without much difficulty, though there is some danger to doing this.)  I did get WMS support working, but I could only get their sample WMS servers to draw.  I suspect that with some more experimentation (or actually reading the instructions), I could get support for different coordinate systems and WMS working better.

There are a few other features I quickly discovered that are worth mentioning.  As you would expect, QGIS includes good support for custom symbolization of vector data and exporting a map to an image file.  There is also a neat feature where you can reclaim real estate on the map by pressing CTRL + T / T to hide and show the toolbars.  Finally - and I've saved the best for last - if you're in the market for something to create and edit your own Shapefile, QGIS does this nicely.  Unlike ArcMap, where feature editing, though powerful, is about as natural as breathing underwater, QGIS makes it intuitive.  (Shameless plug: You can create and edit Shapefiles in FalconView via the GIS Editor, which is essentially the same as doing it in ArcMap.)

I very much want to revisit, in a future posting here, the advanced programmatic features of Quantum GIS.  There is support for an interactive Python console (I'm a Python novice, so maybe this will be the push I need to learn it better), and there is support for plugins to QGIS.  I spent about 30 minutes poking around on the QGIS web site, wiki, and Internet in general to learn more about plugin capabilities and how to write one, but I couldn't find anything useful.  Some quick e-mail correspondence with a QGIS developer pointed me here.  (And there we find yet another advantage to using this active open source project: it's not too hard to find help from the project volunteers.)

I'll sum up my initial impressions with Quantum GIS by saying that it's a maturing product worthy of applause.  It does most of what it's supposed to do with ease, and I suspect that with a little more experience, I could overcome the hurdles I experienced.  I enjoyed playing with the tool, which says a lot since I have a pretty short attention span.

If you want a quick demo of the basic capabilities, check out this video.

Thursday, April 10, 2008

Creating a Mosaic of FalconView DTED in the GIS Editor

In my previous post on the subject of DTED in the FalconView ArcGIS Editor, I mentioned that bringing DTED into the GIS Editor is a two step process.  The FalconView user first has the Editor create a feature class that shows all the DTED tiles in the FalconView map data manager; the second step is to select the area of interest and create a DTED mosaic from that area.  The previous post examined the creation of a feature class to show DTED coverage.  This post will examine the technical details of how the GIS Editor uses the geoprocessing capabilities of ArcGIS to add a nicely rendered DTED mosaic to the map.

(Please note that the code below is still somewhat experimental.  Comments and suggestions are welcome.  As FalconView 4.2 is still under development, there are some discussions within Georgia Tech and our user community as to how the final user experience for this will happen.)

1. Get a list of DTED tiles to add to the mosaic

The code block below shows how to enumerate all features the layer of FalconView DTED coverage, collecting a list of DTED tiles which will later be passed to the Mosaic geoprocessing tool.  Note that the only end result from this code block is sMosaicInput, the colon separated list of tiles to add to the mosaic (e.g. "C:\dted\w086\n35.dt1;C:\dted\w086\n34.dt1").  I have highlighted the two most interesting parts of the code: the part that builds the lookup table of selected features and the part that builds sMosaicInput from the features in the lookup table.
string sMosaicInput = null;

// ... some redacted code ...

// The block below populates sMosaicInput

try
{
// Create a lookup table of selected features

IFeatureSelection pFeatureSelection = (IFeatureSelection)pLayer;
Dictionary<int, object> dictSelectedIDs = new Dictionary<int, object>();

IEnumIDs pEnumIDs = pFeatureSelection.SelectionSet.IDs;
for (int id = pEnumIDs.Next(); id != -1; id = pEnumIDs.Next())
dictSelectedIDs[id] = null;

// ... some redacted code ...

// Enumerate all features, adding selected paths to the list

IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;
IFeatureCursor pFeatureCursor = pFeatureLayer.FeatureClass.Search(null, true);
int iPathField = pFeatureCursor.FindField("Path");

for (IFeature pFeature = pFeatureCursor.NextFeature(); pFeature != null; pFeature = pFeatureCursor.NextFeature())
{
if (dictSelectedIDs.ContainsKey(pFeature.OID))
{
string sPath = (string)pFeature.get_Value(iPathField);
if (sMosaicInput == null)
sMosaicInput = sPath;
else
sMosaicInput += ';' + sPath;
}
}
}
catch (Exception ex)
{
// ... some redacted code ...

}
2. Use geoprocessing tools to create the mosaic raster

Once we have compiled the string containing the paths of the DTED tiles to be added to the mosaic, creating the mosaic is a simple task requiring the use of two geoprocessing tools.  In the code below, the first highlighted section executes a CreateRasterDataset geoprocessor, creating a scratch raster dataset in a scratch file geodatabase.  The second highlighted lines execute the Mosaic geoprocessor.
// Create a new raster dataset to hold the mosaic

string sRasterName = "dted_mosaic" + Guid.NewGuid().ToString("N");
const string PIXEL_TYPE = "16_BIT_SIGNED"; // Matches DTED 1 & 2 in ArcMap
const int NUMBER_OF_BANDS = 1; // Matches DTED 1 & 2 in ArcMap

CreateRasterDataset createRasterDataset = new CreateRasterDataset(TempGeodatabaseUtils.TempFileGeodatabase, sRasterName, PIXEL_TYPE, NUMBER_OF_BANDS);

Geoprocessor geoprocessor = new Geoprocessor();
geoprocessor.Execute(createRasterDataset, null);
for (int j = 0; j < geoprocessor.MessageCount; j++)
   LoggingUtils.TraceMessage(geoprocessor.GetMessage(j));

// Create the mosaic

Mosaic mosaic = new Mosaic(sMosaicInput, TempGeodatabaseUtils.TempFileGeodatabase + '\\' + sRasterName);
mosaic.mosaic_type = "BLEND";
mosaic.nodata_value = -32767;  // NoData for DTED 1 & 2 in ArcMap

geoprocessor = new Geoprocessor();
geoprocessor.Execute(mosaic, null);
for (int j = 0; j < geoprocessor.MessageCount; j++)
   LoggingUtils.TraceMessage(geoprocessor.GetMessage(j));
3. Symbolize the mosaic and add it to the map

Once the mosaic is created, it must be symbolized before it is added to the FalconView map.  After some experimentation, we found that a simple color ramp shows the features of the DTED well.  (Note that a possible improvement to this would be to add a hillshade effect.)  The first section of highlighted code below shows how to create raster statistics programmatically.  The raster statistics will be used to determine the parameters of the renderer.  The highlighted section shows how to build a simple stretch color ramp.  The last lines show how to create a raster layer from the raw raster, apply the renderer, and add the raster to the map control.
// Set up a raster data object from the mosaic

IWorkspaceFactory pWorkspaceFactory = new FileGDBWorkspaceFactory();
IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(TempGeodatabaseUtils.TempFileGeodatabase, hWnd);
IRasterWorkspaceEx pRasterWorkspace = (IRasterWorkspaceEx)pWorkspace;
IRasterDataset pRasterDataset = pRasterWorkspace.OpenRasterDataset(sRasterName);

// Create raster statistics so we can find the highest and lowest points on the raster for the color ramp

IRasterBandCollection pRasterBandCollection = (IRasterBandCollection)pRasterDataset;
IEnumRasterBand pEnumRasterBand = pRasterBandCollection.Bands;
IRasterBand pRasterBand = pEnumRasterBand.Next();
pRasterBand.ComputeStatsAndHist();
IRasterStatistics pRasterStatistics = pRasterBand.Statistics;

// Create a renderer to be used with the raster layer

IRgbColor pFromColor = new RgbColorClass();
pFromColor.Red = 254;
pFromColor.Green = 254;
pFromColor.Blue = 255;

IRgbColor pToColor = new RgbColorClass();
pToColor.Red = 0;
pToColor.Green = 0;
pToColor.Blue = 255;

IAlgorithmicColorRamp pAlgorithmicColorRamp = new AlgorithmicColorRampClass();
pAlgorithmicColorRamp.FromColor = pFromColor;
pAlgorithmicColorRamp.ToColor = pToColor;
pAlgorithmicColorRamp.Algorithm = esriColorRampAlgorithm.esriCIELabAlgorithm;
pAlgorithmicColorRamp.Size = 50;

bool bOK;
pAlgorithmicColorRamp.CreateRamp(out bOK);
LoggingUtils.Assert(bOK, "Failed to create color ramp");

IRasterStretchColorRampRenderer pRasterStretchColorRampRenderer = new RasterStretchColorRampRendererClass();
pRasterStretchColorRampRenderer.BandIndex = 0;
pRasterStretchColorRampRenderer.LabelHigh = pRasterStatistics.Maximum.ToString();
pRasterStretchColorRampRenderer.LabelLow = pRasterStatistics.Minimum.ToString();
pRasterStretchColorRampRenderer.ColorRamp = pAlgorithmicColorRamp;

// Create the raster layer and add it to the map

IRasterLayer pRasterLayer = new RasterLayerClass();
pRasterLayer.CreateFromDataset(pRasterDataset);
pRasterLayer.Name = "FalconView DTED Mosaic";
pRasterLayer.Renderer = (IRasterRenderer)pRasterStretchColorRampRenderer;

pMapControl2.AddLayer(pRasterLayer, 0);







Wednesday, April 2, 2008

Consuming FalconView DTED in the GIS Editor

The FalconView 4.2 ArcGIS Editor will have the capability of pulling DTED from FalconView into an ESRI ArcMap document so that users can perform geoprocessing tasks with the terrain data. In this blog posting, I will examine this capability that will be in FalconView 4.2 and will give an overview of how this capability was implemented technically. The reader will learn how to retrieve coverage data from the FalconView map data manager and how to create a feature class from this coverage data.

Users add DTED tiles to their FalconView coverage database through the FalconView map data manager. Once the DTED is loaded into FalconView, users can draw elevation data visually, use elevation data in route planning, calculate terrain masking of radar systems, and so on. While FalconView has numerous uses for terrain elevation data, ArcGIS has a generalized geoprocessing framework which allows GIS intelligence officers to perform a wide variety of specialized calculations for their mission planning. The integration of FalconView and ArcGIS will provide a simplified framework for mission planning which requires geoprocessing with terrain data.

The screen capture to the left shows FalconView drawing DTED two ways (the DTED view is zoomed out to 25% to remove detail so that the screen capture is suitable for posting to a public forum.) The black and white background map is the normal DTED base map drawn by FalconView. The blue and white tile in the center is a mosaic of four DTED tiles being drawn by the GIS Editor. The FalconView base map is optimized for a viewing experience that shows the general layout of the terrain. The GIS Editor splotch in the center of the map is optimized for a viewing experience that highlights the high and low areas of the terrain.

For GIS Editor users, creating a mosaic of FalconView DTED is a two step process. First, the user creates a vector overlay of polygons which describe the DTED coverage available from FalconView. Next, the user selects the DTED tiles he wishes to turn into a mosaic and creates the mosaic. The coverage overlay and the mosaic are both automatically created by the GIS Editor and are stored in temporary geodatabases. The screen capture to the right shows the tiles of DTED available on my development system. Notice the little raster blotch in the Eastern US. This is a mosaic of four DTED tiles. (It may be necessary to click on the image to the right to see the detail of the overlay.)

Let's explore the technical steps necessary to bring this DTED from FalconView into ArcGIS. (In this blog post we will not examine the steps used to create the DTED mosaic, though we may explore that in a follow-up post.)

1. Creating an ArcGIS feature class to contain the FalconView DTED

ArcGIS always holds feature classes in some sort of geodatabase. In the case of the FalconView DTED, we load the coverage data into a scratch file geodatabase.
// Create the output workspaces

IWorkspaceFactory pOutWorkspaceFactory = new FileGDBWorkspaceFactory();
IWorkspace pOutWorkspace = pOutWorkspaceFactory.OpenFromFile(TempGeodatabaseUtils.TempFileGeodatabase, hWnd);
IFeatureWorkspace pOutFeatWorkspace = (IFeatureWorkspace)pOutWorkspace;
The temporary file geodatabase referred to above is created during the call to the TempFileGeodatabase property, if it didn't exist before.

After we get the reference to the output feature workspace, we must create the Fields that will be stored with the feature class. The fields we create include a shape field to hold the shape of the polygon, a path field which holds the path to the DTED tile, and a level field which indicates the DTED level (0, 1, 2, 3, etc.).
// Create the fields object which will hold the fields for the feature class. I have put much of the field creation
// logic in its own block to demonstrate that the only object that matters after the fields are built is pFields.

IFields pFields = new FieldsClass();

{
// Set up the polygon shape field for the feature class

IFieldsEdit pFieldsEdit = (IFieldsEdit)pFields;

IField pField = new FieldClass();
IFieldEdit pFieldEdit = (IFieldEdit)pField;
pFieldEdit.Name_2 = "Shape";
pFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;

IGeometryDef pGeometryDef = new GeometryDefClass();
IGeometryDefEdit pGeometryDefEdit = (IGeometryDefEdit)pGeometryDef;
pGeometryDefEdit.GeometryType_2 = esriGeometryType.esriGeometryPolygon;

ISpatialReferenceFactory2 pSpaRefFact2 = new SpatialReferenceEnvironmentClass();
IGeographicCoordinateSystem pGeoCoordSys = pSpaRefFact2.CreateGeographicCoordinateSystem(4326); // 4326 is esriSRGeoCS_WGS1984
ISpatialReference pSpatialReference = (ISpatialReference)pGeoCoordSys;

// Set the domain on the coordinate system
// This prevents the error, "The XY domain on the spatial reference is not set or invalid."

ISpatialReferenceResolution pSpatialReferenceResolution = (ISpatialReferenceResolution)pSpatialReference;
pSpatialReferenceResolution.ConstructFromHorizon();

pGeometryDefEdit.SpatialReference_2 = pSpatialReference;
pFieldEdit.GeometryDef_2 = pGeometryDef;
pFieldsEdit.AddField(pField);

// Add the Path field

pField = new FieldClass();
pFieldEdit = pField as IFieldEdit;
pFieldEdit.Name_2 = "Path";
pFieldEdit.Type_2 = esriFieldType.esriFieldTypeString;
pFieldsEdit.AddField(pField);

// Add the Level field

pField = new FieldClass();
pFieldEdit = pField as IFieldEdit;
pFieldEdit.Name_2 = "Level";
pFieldEdit.Type_2 = esriFieldType.esriFieldTypeInteger;
pFieldsEdit.AddField(pField);

// Add other required fields to the feature class

IObjectClassDescription pObjectClassDescription = new FeatureClassDescriptionClass();

for (int i = 0; i < pObjectClassDescription.RequiredFields.FieldCount; i++)
{
pField = pObjectClassDescription.RequiredFields.get_Field(i);
if (pFieldsEdit.FindField(pField.Name) == -1)
pFieldsEdit.AddField(pField);
}
}
Note the loop an the end of this block of code that ensures all required fields are added to the feature class. When creating a feature class from scratch, it is important to add all of the required fields to the class. In a debug session, I see this block of code skips adding the SHAPE field, which we manually created at the beginning of the block, but adds the required OBJECTID field. The image below is a snapshot of the attributes table of the new feature class as seen in the FalconView GIS Editor. If you click on the image, you will see the fields we added to the feature class, along with some of the data which we will add in the next step.
The next step is to create the feature class, which is relatively simple. The code below also shows how we set up some references we will use to add features to the feature class.
// Create the feature class

string sFeatureClassName = "dted" + Guid.NewGuid().ToString("N");
IFeatureClass pFeatureClass = pOutFeatWorkspace.CreateFeatureClass(
sFeatureClassName, pFields, null, null, esriFeatureType.esriFTSimple, "Shape", null);

// Get an insert cursor and a feature buffer

IFeatureCursor pFeatureCursor = pFeatureClass.Insert(true);
IFeatureBuffer pFeatureBuffer = pFeatureClass.CreateFeatureBuffer();

// Get the indicies for the fields of interest

int iPathFieldIndex = pFeatureBuffer.Fields.FindField("Path");
int iLevelFieldIndex = pFeatureBuffer.Fields.FindField("Level");
Note that creating a feature class and adding fields to the class can also be done using several geoprocessing objects, but, once you understand the steps above, it's simpler to do it as described above.

2. Pulling DTED coverage tiles from the FalconView Map Data Manager

The block of code below shows how to enumerate all FalconView map coverage, creating polygons from DTED tiles and adding them to our feature class. Note that there are three nested loops here. The outer loop enumerates over all of the data sources in the FalconView coverage database. The next inner loop enumerates over all of the DTED levels on the data source under consideration (we're only enumerating levels 1 through 3 for now). The innermost loop enumerates over all of the tiles of the particular DTED level on the particular data source. A polygon feature is created for each tile of DTED enumerated in the innermost loop and is added to the feature class using the insert cursor.
// Initialize the coverage rowset

ICoverageRowset coverageRowset = new CoverageRowsetClass();
coverageRowset.Initialize("Dted");

// Enumerate all data sources, cataloging DTED

const int FLUSH_INTERVAL = 100; // Experiment with this value as needed
int iRowsProcessed = 0;
Dictionary<string, object> dictFoldersAdded = new Dictionary<string, object>();

IDataSourcesRowset dataSourcesRowset = new DataSourcesRowsetClass();
dataSourcesRowset.SelectAll(1 /* Online Only */);
while (dataSourcesRowset.m_Identity > 0)
{
string sLocalFolderName = dataSourcesRowset.m_LocalFolderName;

// Enumerate all DTED on this data source

for (int iDtedLevel = 1; iDtedLevel <= 3; iDtedLevel++)
{
coverageRowset.SelectByGeoRectAndDS(dataSourcesRowset.m_Identity, iDtedLevel, -90, -180, 90, 180);

while (coverageRowset.m_LocationSpec != "")
{
// Get the path to this DTED tile and skip it if it's already added

string sPath = sLocalFolderName + '\' + coverageRowset.m_LocationSpec;

if (dictFoldersAdded.ContainsKey(sPath.ToLower()))
goto next_row; // Sometimes can happen when map data server acts up
dictFoldersAdded[sPath.ToLower()] = null;

// Get the bounds of this DTED tile

double dllLat, dllLon, durLat, durLon;
coverageRowset.GetBounds(out dllLat, out dllLon, out durLat, out durLon);

// Create a polygon from the coverage rectangle

IPolygon pPolygon = new PolygonClass();
IPoint pPoint = new PointClass();
IPointCollection pPointCollection = pPolygon as IPointCollection;
object missing = Type.Missing;

pPoint.X = dllLon;
pPoint.Y = dllLat;
pPolygon.FromPoint = pPoint;
pPoint.X = durLon;
pPointCollection.AddPoint(pPoint, ref missing, ref missing);
pPoint.Y = durLat;
pPointCollection.AddPoint(pPoint, ref missing, ref missing);
pPoint.X = dllLon;
pPointCollection.AddPoint(pPoint, ref missing, ref missing);
pPolygon.Close();

// Add the polygon to the feature class

pFeatureBuffer.Shape = pPolygon;
pFeatureBuffer.set_Value(iPathFieldIndex, sPath);
pFeatureBuffer.set_Value(iLevelFieldIndex, iDtedLevel);
pFeatureCursor.InsertFeature(pFeatureBuffer);

// Flush the data from time to time

if (++iRowsProcessed % FLUSH_INTERVAL == 0) pFeatureCursor.Flush();

// Move to the next DTED coverage rowset

next_row: coverageRowset.MoveNext();
}
}

// Move to the next FalconView data source

dataSourcesRowset.MoveNext();
}
The feature cursor is flushed from time to time during its population and again after all of the DTED tiles have been added to the feature class. It's important after this step to manually release all of your COM objects so that the lock on the feature class in the geodatabase is removed. The last step is to create a feature layer from the feature class. (All of these steps are shown below.) The feature layer may be directly added to the map.
// Flush the feature cursor and release COM objects

pFeatureCursor.Flush();

Marshal.ReleaseComObject(pFeatureBuffer);
Marshal.ReleaseComObject(pFeatureCursor);
Marshal.ReleaseComObject(pFeatureClass);
Marshal.ReleaseComObject(pOutFeatWorkspace);
Marshal.ReleaseComObject(pOutWorkspace);
Marshal.ReleaseComObject(pOutWorkspaceFactory);

// Return a DTED feature layer

IWorkspaceFactory pWorkspaceFactory = new FileGDBWorkspaceFactoryClass();
IWorkspace pWorkspace = pWorkspaceFactory.OpenFromFile(TempGeodatabaseUtils.TempFileGeodatabase, hWnd);
IFeatureWorkspace pFeatureWorkspace = (IFeatureWorkspace)pWorkspace;
pFeatureClass = pFeatureWorkspace.OpenFeatureClass(sFeatureClassName);

IFeatureLayer pFeatureLayer = new FeatureLayerClass();
pFeatureLayer.Name = FV_DTED_LAYER_NAME;
pFeatureLayer.FeatureClass = pFeatureClass;