Wednesday, 10 August 2016

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:

"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."

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:


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:

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


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:

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
p1.number <> p2.number
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. 

1 comment:

  1. I think you could also just used SqlSheet that comes with Autodesk Infrastructure Administrator (C:\Program Files\Autodesk\Autodesk Infrastructure Administrator 20xx\bin) to performe all SQL Statements once connected to your SQLite-File.

    Kind regards,