Tuesday 23 August 2016

Python - update georeferencing information in multiple files

You need to install GDAL Python bindings first.

https://pythongisandstuff.wordpress.com/2011/07/07/installing-gdal-and-ogr-for-python-on-windows/

Here is a basic script as template:

 import os  
 import fnmatch  
 from datetime import date  
 from osgeo import gdal  
 from gdalconst import *  
 # logfile-path  
 logfile = "c:/temp/geotifheader_update.log"  
 #path to GeoTIFs  
 tif_filefolder = "S:/MG_DATEN/Rasterdaten"  
 #shift for all images - we just add a fixed shift to all images   
 shift_x = 2000000.950  
 shift_y = 999999.800  
 # get all TIF files in folder, including subfolders  
 def list_files(dir):  
   r = []  
   for root, dirs, files in os.walk(dir):  
     #for name in files:  
     for name in fnmatch.filter(files, "*.tif"):  
       r.append(os.path.join(root, name))  
   return r  
 #write message to log file  
 def writelog(message):  
   with open(logfile, 'a') as myfile:  
     myfile.write(message)  
 #init log file  
 def initlog():  
   with open(logfile, 'a') as myfile:  
     myfile.write(str(date.today())+"\n")  
 #update georeferencing information for each image  
 def setHeader(filepath):  
   dataset = gdal.OpenEx(filepath, GA_Update)  
   geotransform = dataset.GetGeoTransform()  
   if not geotransform is None:  
     print "**file**: ", filepath  
     print 'Origin = (', geotransform[0], ',', geotransform[3], ')'  
     print 'Pixel Size = (', geotransform[1], ',', geotransform[5], ')'  
     ulx = geotransform[0];  
     uly = geotransform[3];  
     #simple check whether image should be updated or not  
     if (ulx > 1000000 or ulx == 0 or uly > 1000000 or uly == 0):  
       m = "current insertion point out of bounds: " + str(ulx) + "," + str(uly)+"\n"  
       print m  
       writelog(m)  
     else:  
       ulx = ulx + shift_x  
       uly = uly + shift_y  
       lrx = ulx + (dataset.RasterXSize * geotransform[1])  
       lry = uly + (dataset.RasterYSize * geotransform[5])  
       gt = [ulx, (lrx - ulx) / dataset.RasterXSize, 0, uly, 0, (lry - uly) / dataset.RasterYSize]  
       dataset.SetGeoTransform(gt)  
       m = "image has been updated.\n"  
       print m  
       writelog(m)  
   else:  
     m = "not a geotiff.\n "  
     print m  
     writelog(m)  
 files = list_files(tif_filefolder)  
 initlog()  
 for file in files:  
   writelog("-----------------------\n")  
   writelog(file+"\n")  
   setHeader(file)  
   writelog("-----------------------\n")  
 '''  
 for single file:  
 setHeader("S:/GIS/DTMAV/Relief/relief_dtmav30.tif")  
 '''  
 # https://gis.uster.ch/dokumentation/datenkonvertierung/gdal-tools  
 # set path=C:\Program Files\GDAL;C:\Python27  

Thursday 18 August 2016

MapGuide/AIMS - modifying maps and layers as XML

Often I export MapGuide map definitions or layers as XML in order to add or remove content. Most often it is faster than using Studio. The modified XML file is uploaded with MapAgent and "SetResource". Specifically map definitions contain hundreds of layers and often layers need to be added or removed. I just realized that instead of deleting layers and uploading the map definition I can just put the layers in questions within XML comment tags and upload the file again. In case I have made a mistake and "deleted" a layer which I didn't want to delete I can go back to my XML file, remove the comment tags and upload the file again. Unfortunately the comments are filtered out when uploading - so when you save the map definition again as XML file from Studio all comments are gone. 
By the way - when you put comments into Map layer definition files and open and save the DisplayModel with Map the comments are retained.




 </MapLayer>  
 -->  
  <MapLayer>  
  <Name>GR_Baum</Name>  
  <ResourceId>Library://FS_INFO/WT_GR_I/Layers/GR_Baumkataster/GR_Baum.LayerDefinition</ResourceId>  
  <Selectable>false</Selectable>  
  <ShowInLegend>true</ShowInLegend>  
  <LegendLabel>GR_Baum</LegendLabel>  
  <ExpandInLegend>false</ExpandInLegend>  
  <Visible>true</Visible>  
  <Group>GR_Baumkataster</Group>  
  </MapLayer>  
  <!--  
  <MapLayer>  
  <Name>SW_Projektperimeter</Name>  
  <ResourceId>Library://WMS_IN/IS/Layers/SW_Projektperimeter.LayerDefinition</ResourceId>  
  <Selectable>false</Selectable>  
  <ShowInLegend>true</ShowInLegend>  
  <LegendLabel>SW_Projektperimeter</LegendLabel>  
  <ExpandInLegend>false</ExpandInLegend>  
  <Visible>false</Visible>  
  <Group>Leitungskataster</Group>  
  </MapLayer>  
  -->  
  <MapLayer>  

Saturday 13 August 2016

Lost in translation - even worse....

Just created network deployment for Map 2017 and noticed that ReCap description hasn't changed much (see screenshot). If you do not speak german than you really miss out on a ridiculous translation. 

ReCap - "Realitätserfassung in einer Vorbereitungsumgebung" - good luck with that!
link


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.