Spatial Queries and Map

Once in while you have CAD or GIS data and you need to link features but there are no attributes to do that. Map's spatial functions such as MapOverlay cannot always be used to establish a spatial relationship. Here are two examples:

#1
"I have a field of points ... What I would like is to be able to analyse the neighbours within the footprint (1 meter) of each recording point to determine the range elevation differences. So for instance, the circled point in the snip, has two neighbours and the deviations in elevation from the point being analysed are +0.007 and -0.010m. I would like to somehow record these values. Some points won't have any close neighbours so they would have no records against them."

#2
There are 8000 weld joints along a pipe. Weld joints are marked by means of a block including some attributes. Next to the pipe are 800 drill points. They are also marked by means of a block including certain attributes. The task was find the closest drill point for each weld joint. The distance between weld joints and drill point is not uniform.

I don't know whether you could solve the tasks with out-of-the-box functionality in QGIS or ArcGIS. But as both support Python it is probably not so difficult to find the answers writing a few lines of code. Unfortunately Map doesn't offer any useful scripting - VBA and Lisp only support a subset of the Map API.

If you need to solves tasks like the ones above and you do not have anything else at hand than Map then could try your luck with SQLite. 

SQLite itself is a file based database. There are tools freely available to connect to SQLite files and to perform SQL queries. There is also an extension for spatial data which means you can perform spatial queries. Map can import and export SQLite Spatial which allows you to run spatial queries against Map data without having to install a full blown (spatial) database system such as Oracle Spatial, SQL Server or PostGIS.

When working with spatial data in SQLite you can use Spatialite_GUI:

Download:
http://www.gaia-gis.it/gaia-sins/
http://www.gaia-gis.it/gaia-sins/windows-bin-x86/spatialite_gui-4.3.0a-win-x86.7z

Unzip file using 7-zip and execute the program "spatialite_gui.exe" - no installation required. 

I'm going to describe how to export to SQLite and how to perform a simple spatial query by taking the data from example 1 above. The drawing contains CIVIL objects which I exported to SDF using CIVIL command EXPORTTOSDF. I then imported the SDF into Map (command: mapimport) including the attribute data. As starting point I got a drawing containing plain AutoCAD points with attribute data attached. 

1. open drawing file

2. check the AC points - they have Map object data attached, each point has a point number and an elevation value: 

AutoCAD points and some attribte data attached (Map object data table)

3. export points to SQLite Spatial file, command: mapexport

make sure that you export your points as Point feature class (not as Point+Line+Polygon), also export the attached attributes, you can also set a name for the table the features will be exported to:

MAPEXPORT - export options

   
4. run SpatialLiteGUI   

5. connect to SQLite file

You should see the following message:
FDO-OGR detected;activating FDO-OGR auto-wrapping...
- Virtual table: tablename
...

Spatialite_Gui message about Map/FDO geometries

Spatialite_gui supports "native" SQLite spatial geometries as well as Map (FDO) geometries. The message tells you, that Map/FDO geometries were found in the database and that Spatialite_GUI has created a wrapper around those geometries by means of a "virtual table". With that mechanism in place Map/FDO geometries can be used as they were "native" SpatiaLite geometries. 

6. the application shows all tables found in the SQLite file, you will not only see the table containing the points you just exported but also some further tables created by Map containing some Map specific metadata. Below "User data" you will see at least two tables. Table number one will have the same name as given in Mapexport dialog box, the second table will have the same name but with the prefix "fdo". The latter one is also marked with "chain" symbol - meaning it is "virtual table". Both tables have basically the same content - but only the second table can be used for spatial queries:

Spatialite_GUI

Perform a quick test and type the following SQL statement into the input area and then execute the query:

select * from mypoints

the result will look like that:

1 BLOB sz=32 UNKNOWN type -2.638000 5522
2 BLOB sz=32 UNKNOWN type -2.887000 5523
3 BLOB sz=32 UNKNOWN type -2.930000 5524

Query result

The GEOM column contains the geometry of the feature - but the content of the geometry column is just shown as "BLOB" and marked as "unknown" - the geometry in the original table is stored in a way SQL can not access it.

If you do the same against the second table the result will look slightly different:

select * from fdo_mypoints

1 BLOB sz=60 GEOMETRY -2.638000 5522
2 BLOB sz=60 GEOMETRY -2.887000 5523
3 BLOB sz=60 GEOMETRY -2.930000 5524

It is still only shown as "BLOB" but SQL also recognises that those "BLOBS" store some kind of Geometry. 

To access the geometry and to apply spatial queries you have to use the second table.
Here is another example - show the X and Y coordinates of the points:

select number, ST_X(geom) as X, ST_Y(geom) as Y from fdo_mypoints

5522 15950702.370320 2537019.106847
5523 15950700.846391 2537017.708832
5524 15950699.357462 2537016.129788

If you execute the SQL statement against the other table you won't get any useful result:

select number, ST_X(geom) as X, ST_Y(geom) as Y from mypoints

5522 NULL NULL
5523 NULL NULL
5524 NULL NULL

7. Before we can perform the spatial query we need to copy our (virtual) table. I don't know why but when I run the SQL statement (as shown in next step) against my fdo_mypoints table I get a wrong result (too little rows are returned - it seems the cross join doesn't work on the virtual table). To copy the table the following SQL statement needs to be executed:

CREATE TABLE new_table_name AS SELECT * FROM existing_table_name

CREATE TABLE fdo_mypoints_copied AS SELECT * FROM fdo_mypoints

Afterwards you need to refresh the tree view (context menu for root node >> "Refresh") in order to see the newly created table.

7. To link points within distance of 1m  we can use the following statement:

select 
p1.number as p1_num, 
p2.number as p2_num,
p1.elevation as p1_ele,
p2.elevation as p2_ele,
ST_Distance(p1.geom, p2.geom) as dist_poi,
p1.elevation-p2.elevation as diff_ele,
p1.geom as p1_geom
from fdo_mypoints_copied as p1 cross join fdo_mypoints_copied as p2
where 
p1.number <> p2.number
and
dist_poi < 1
order by p1_num

As all points are in one table we need a cross join to process a point against all other points in the same table. We then filter out all points which are equal (p1.number <> p2.number) and where the distance between points is >= 1.

Now we can check the result - if the query returns roughly the number of rows we would expect - then we are ready to export the result. Just right-click anywhere in the table view choose "Export ResultSet --> as Shapefile".

Keep in mind: column names in SHP file are limited to 10 characters, your select statement should set columns names accordingly.

8. Add SHP file to Map
When you open the table view for the SHP file in Map and you click on one of the (numeric) columns you might get an error message (data type 'FdoDataType_Int64' not supported) and the table view gets blank. Just right-click on layer in Display Manager and choose "Export layer data as SDF", then remove SHP layer from drawing and add SDF file instead.


It's not that complicated, is it? Some knowledge of SQL and spatial queries is required though. There are also some peculiarities with SQLite / SQLite-Spatial but as soon as you aware of them the whole process can be done in a few minutes. 

Map API - MgFeatureReader - improve speed

The developer example for reading features with MgFeatureReader looks like that:

 while (featureReader.ReadNext())  
    {  
      featureCount++;  
      propertyCount = featureReader.GetPropertyCount();  
      for (int i = 0; i < propertyCount; i++)  
      {  
        propertyName = featureReader.GetPropertyName(i);  
        propertyType = featureReader.GetPropertyType(propertyName);  
        PrintPropertyValueFromReader(featureReader, propertyType, propertyName, ref sb);  
      }  
    }  

(form MapGuide Api Reference).

At least with Map this will be slow. Performance can be improved significantly if  the GetPropertyName and GetPropertyType calls are removed from within the FOR loop.

Example: a small application in order to displays feature attributes in a DataGrid. 

DataGrid

The sample SDF file contains 110'000 features and has 11 columns (text and numbers).

Using code similar to the one given above it takes 19 sec to populate and display the datagrid.

By removing the two calls the process takes less then 1 second and is even slightly faster then opening  Map's data table (mapdatatable):


Befehl: MAPDATATABLE2
 Finished in 0 seconds.
Befehl: MAPDATATABLE3

 Finished in 19 seconds.

All it takes is to get property name and property type for each column only once and store this information in an array (colProps) before iterating over the FeatureReader:

 ...  
  while (featureReader.ReadNext())  
       {  
         counter++;  
         newRow = dt.NewRow();          
         for (int i = 0; i < colProps.Length; i++)  
         {  
           propertyName = colProps[i].PropertyName;  
           propertyType = colProps[i].PropertyType;  
                          ....  

Map 2013, SP2

Convert Map/FDO lines to ACAD lines

There are situations when you have Map /FDO features but want to use an AutoCAD function, such as dimensions. With Map Industry Models you can apply dimensions on IM features but it is less comprehensive than AutoCAD dimensions. So I wrote - on request - a little tool to convert Map /FDO lines to ACAD polylines "on the fly".

Today another use case for this came up. Sometimes lines from one IM feature source (such as road axis) need to be "copied" to a second feature class (such as cycle paths).  There is no "copy" for IM / Map / FDO geometries in Map. Now the workflow is:
- convert Map / FDO lines to polylines in your drawing ("on the fly")
- create new IM feature and use polyline as geometry input

Map 2013, SP2

Feature with wrong geometry type for feature class

We had one feature (Industry Model) which was not displayed in Map. 1 Click Maintenance reported an unsupported geometry type for a polygon feature class. A closer look unveiled that the feature had been migrated as line feature and not as polygon feature:
  

Line string feature in a polygon feature class (Map Industry Model)

The geometry itself was correct - so it could not be "repaired" in Oracle Spatial. Here is one way to solve the issue:


1. connect via FDO Oracle and add feature class to Map 
   FDO Oracle will create a layer containing lines and polygons - the feature in question will therefore be displayed in Map.

Adding table via FDO Oracle - line and polygon features are displayed

2. save layer as sdf, add sdf to Map and remove layer based on FDO Oracle from Display Manager
3. rectify geometry in SDF based layer - for instance via "extract geometry" 
4. open form for Industry model feature where geometry type is incorrect
5. link rectified geometry with feature - confirm "Feature has already a geometry - replace?" 

Map 2013, SP2