Showing posts with label spatial query. Show all posts
Showing posts with label spatial query. Show all posts

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:

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



Wednesday, 2 September 2015

Map - and GIS analysis

There was a posting at CAD.DE regarding a specific task to be accomplished by using Map.

Here is a short description:

A drawing file contains two different sets of point features.

point feature set 1 - blocks, including an ID (block attribute), labelled as "b1"
point feature set 2 - points, including an ID (Map object data table) labelled as "s1"

Here is a screenshot:





The task was to create a list where each point of feature class "s1" is linked to the closest point in feature class "b1".

Several suggestions were made how to solve the issue with Map - but either it involved laborious manual work or wouldn't always return the closest point. One suggestion was to use a spatial database to get the desired result.

I then decided to try whether that would be possible without installing any software and using only freely available tools.

Here is a short description of what I did using SQLite and Spatialite:

(1) export both point data sets (including block attributes and Map object data table) to SQLite using _mapexport


(2) download spatialite_gui.exe (http://www.gaia-gis.it/gaia-sins/index.html) - doesn't require an installation


(3) in spatial_lite connect to SQLite file
        Spatialite recognises the Map/FDO SQLite data model and geometry columns. There are accessible as "virtual tables" meaning that they cannot  be modified. It also seems that you cannot run spatial queries on the virtual tables - at least I didn't work for me. But I'm not experienced with spatial lite and might have just done something wrong.


(4) copy data from virtual tables to normal SQLite tables using CTAS (create table as select)
        With spatialite you need to create a geometry column in a very special way. First you have to create the table, afterwards you add the geometry column using the following statement: SELECT AddGeometryColumn(...). As I created my two tables with CTAS I did not get "proper" geometry columns. But for the task in question it didn't matter. The only drawback I noticed with my approach was that I couldn't export to SHP. This requires proper geometry columns.


(5) looking around I found some SQL statements for SQLite for similar tasks, I tried and played around a bit and found the following statement to return correct results:

SELECT
s1.featid as sn_id, b1.nr as b_nr, min(ST_Distance(s1.geom, b1.geom))  as dist, x(s1.geom) as xgeom , y(s1.geom) as ygeom
FROM b1, s1
group by s1.featid
ORDER BY s1.featid ;

(6) I created a new table using the select statement above


(7) to perform a visual check I wanted to load the result into Map again. As I had modified the SQLite database - by adding tables - Map refused to connect to the file. Anyway - even I had been able to load the file into Map it wouldn't have recognized the spatial_lite geometry column as Map/FDO and spatialite use different structures to store geometries. I therefore exported the table as CSV file.

(8) Now as I thought it would take too long to convert the CSV file into a spatial feature class using Map I just loaded the CSV file into QGIS and exported it as SHP data set which I then connected to in Map.

(9) overlaying the SHP file and the original drawing file I could visually check the correctness of the result.








Tuesday, 14 July 2015

Map workflow - separate files for specific area

Situation: ca 7000 tiled raster files in one folder but only a subset for a specific area is required, SHP file containing tiles and filenames and SHP file with area of interest available

raster tiles (green), area of interest (blue)



workflow:

1 - load both SHP files into Map
2 - spatial query to select tiles for area of interest


raster tiles for area of interest


3 - export attribute data for selection


shp with tiles (selection for area of interest) - export attribute data (filenames) 

4 - open attribute data in Excel
5 - extract exact filename with Excel and add addtional parameters (DOS commands), as result you need something like this: 

...
copy 6895_2595_r3x.tif c:\temp\dom
copy 6895_2600_r3x.tif c:\temp\dom
copy 6895_2605_r3x.tif c:\temp\dom
...


further processing in Excel

6 - copy result into new textfile, save as BAT file into folder with raster files and execute



What happens if you do not have a SHP file containg tiles and filenames?

This might be an option to create one with Map:

- load raster files into Map (command: _mapiinsert), choose option to show frame only (do not load thousands of raster files, only let Map display the frames)


_mapiinsert - show frames only


- add AutoCAD layer with image frames to DisplayManager (query current drawing)

- add a label to the newly add layer, choose ".IMAGENAME" as label text



label added to iamge frames - path to file is used as text


- export as SHP (command: _mapexport), export option "Text" : you will get a Point SHP with text insertion points as points and label text as attribute

- use point shape for spatial query