How to add a Google Maps styled layer in Qgis 2

I have been struggling for days to make a very simple and light base-map of the UK showing town circle markers, town names sized according to their population, a couple of country borders and nothing more than a plain white background for lands above sea-level.

After using numerous data sources (Natural Earth, Ordnance Survey, Open Gazetteers, Open Street Map vectors, OSM-GB WMS layers - none of which worked well) I finally came across an acceptable solution...

...which is to make a Google Maps Styled and to tweak Qgis Openlayers Plugin in order to load such a customised layer.


Firstly, make your styled Google Map. You can find mine on Snazzymaps. It is very easy to make your own, just have a go on the Styled Maps Wizard.

Once you have finished your customisations, export the style Javascript Array and place it as the value of the stylez variable in this Html.


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <title>OpenLayers Google Streets Layer</title>
    <link rel="stylesheet" href="qgis.css" type="text/css">
    <link rel="stylesheet" href="google.css" type="text/css">
    <script src="http://maps.google.com/maps/api/js?v=3.3&amp;sensor=false"></script>
    <script src="OpenLayers.js"></script>
    <script src="OlOverviewMarker.js"></script>
    <script type="text/javascript">
        var map;
        var loadEnd;
        var oloMarker; // OpenLayer Overview Marker
        function init() {
            map = new OpenLayers.Map('map', {
                theme: null,
                controls: [],
                units: "m",
                maxResolution: 156543.0339,
                maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34)
            });

            var stylez = [ { featureType: "all", elementType: "all", stylers: [ { visibility: "simplified" } ] } ];
            
            var styledMapType = new google.maps.StyledMapType(stylez, styledMapOptions);           
            
            var gmap = new OpenLayers.Layer.Google(
                "Google Custom",
                { type: 'styled' }
            );
            map.addLayer(gmap);
            
            gmap.mapObject.mapTypes.set('styled', styledMapType);
            gmap.mapObject.setMapTypeId('styled');

            loadEnd = false;
            map.events.register('movestart', map, function() {
                loadEnd = false;
            });
            google.maps.event.addListener(gmap.mapObject, "tilesloaded", function() {
                // wait for tiles to fade in completely
                setTimeout(function() {
                  loadEnd = true;
                },
                200);
            });

            map.setCenter(new OpenLayers.LonLat(0, 0), 2);
            oloMarker = new OlOverviewMarker(map, getPathUpper(document.URL) + '/x.png');
            
        }
    </script>
  </head>
  <body onload="init()">
    <div id="map"></div>
  </body>
</html>


Save this code as google_custom.html in this folder:
C:\Users\youruser\.qgis2\python\plugins\openlayers_plugin\html\

For the sophisticated users:
%userprofile%\.qgis2\python\plugins\openlayers_plugin\html\

After that, enter folder C:\Users\youruser\.qgis2\python\plugins\openlayers_plugin\
and open openlayers_plugin.py with a text editor.

Go around line 124, where you find the #Layer comment and add this line to the layers section, just below the Google satellite entry:

    self.olLayerTypeRegistry.add( OlLayerType(self, 'Google Custom', 'google_icon.png', 'google_custom.html', True) )

Save and close the file and restart Qgis. A new Google custom item appears in the Openlayers plugin Menu.

Enjoy!

Vogliamo anche l'Italia nel registro INSPIRE



Hashtag: #italy4INSPIRE

Premessa

INSPIRE prevede che ogni Stato Membro fornisca almeno un endpoint nazionale per il discovery di metadati.

Ad oggi, la maggior parte degli Stati Membri (23 su 28) ha soddisfatto questo requisito registrando il proprio riferimento nazionale nel geoportale INSPIRE: http://inspire-geoportal.ec.europa.eu/INSPIRERegistry/

In particolare, come si può vedere, alcuni paesi hanno registrato più di un endpoint, come l'Austria, il Belgio e la Lettonia: è infatti possibile registrarne anche più di uno per paese

A differenza di ciò, l'Italia non ha ancora alcun endpoint registrato per il servizio di discovery.

Per questa registrazione è necessaria una semplice comunicazione (email) del National Contact Point INSPIRE (o di qualcuno delegato dal NCP) indirizzata a EC/EEA INSPIRE Team (env-inspire@ec.europa.eu) ed per conoscenza JRC (michael.lutz@jrc.ec.europa.eu).

Domanda

Perché il servizio CSW realizzato da RNDT non è ancora stato registrato come endpoint italiano?

Dal punto di vista normativo, sia il recepimento della Direttiva INSPIRE (Dlgs. 32/2010) che il Codice dell'Amministrazione Digitale riportano che RNDT è il riferimento nazionale in questo contesto: "Il repertorio nazionale dei dati territoriali, [...] costituisce il catalogo nazionale dei metadati relativi ai set di dati territoriali" (Dlgs. 32/2010, art.5) [1].

Dal punto di vista tecnico-operativo i test effettuati nel luglio 2013 e gennaio 2014 dal Joint Research Centre della Commissione Europea (su richiesta dell'Agenzia per l'Italia Digitale) hanno dimostrato che il servizio CSW del RNDT e la quasi totalità dei metadati raccolti sono perfettamente conformi a quanto previsto dai Regolamenti 1205/2008 (metadati) e 976/2009 (servizi di rete) della Commissione Europea, nonché alle relative Technical Guidelines (1.2 del 2010 per i metadati, e 3.1 del 2011 per i servizi di discovery).

In particolare il test effettuato a gennaio 2014 ha riportato 4412 metadati "passed" e 412 "passed with warnings" su un totale di 5540 metadati sottoposti ad harvesting (nel RNDT i metadati disponibili sono 6143).

Il livello di conformità rispetto a INSPIRE è quasi totale per i metadati di dataset e serie (4415 su 4462).
Questo è un risultato importante ed è da notare che risulta essere migliore rispetto ai risultati ottenuti da altri Stati Membri.

Il report completo è disponibile a questo indirizzo: http://inspire-geoportal.ec.europa.eu/resources/sandbox/INSPIRE-dc160d85-7f54-11e3-9486-d8d3855bd8fc_20140117-095358/services/1/PullResults/


Sottolineiamo che è importante che la registrazione del servizio sia fatta al più presto perché:
  1. la disponibilità dei metadati italiani nel catalogo europeo serve a dare visibilità alle informazioni territoriali esistenti in Italia, il tutto proiettato a
    1. supportare le politiche ambientali nazionali e comunitarie
    2. favorire la conoscenza e la promozione del nostro territorio;
  2. l'iniziale disponibilità di metadati potrà innescare un processo virtuoso spingendo gli enti pubblici di ogni livello a conferire i metadati all'RNDT per far conoscere le attività dell'amministrazione su scala internazionale;
  3. per incentivare la realizzazione di servizi innovativi da parte di professionisti, consulenti e PMI locali da offrire agli enti locali sulla base della disponibilità di dati;
  4. per istanziare il ruolo del "nodo" Italia all'interno della rete;
  5. per dare riconoscimento e visibilità alle persone che, su scala diversa, hanno attivamente operato per la realizzazione dell'infrastruttura e dei servizi.

Conclusioni

Alla luce di queste considerazioni, esortiamo il NCP INSPIRE italiano a comunicare al più presto al JRC l'indirizzo del servizio CSW di RNDT affinché questo venga registrato come primo endpoint italiano in INSPIRE.



Firmatari (in ordine alfabetico)
  • Giovanni Allegri
  • Roberto Angeletti, ExportToCanoma blog
  • Andrea Antonello
  • Fulvio Ananasso, Stati generali dell'innovazione
  • Associazione italiana per l'informazione geografica libera - GFOSS.it
  • Associazione Stati Generali dell'Innovazione
  • Ugo Bonelli, Stati generali dell'innovazione
  • Andrea Borruso
  • Stefano Campus
  • Giovanni Ciardi
  • Piergiorgio Cipriano
  • Bruno Conte, Stati generali dell'innovazione, Social4Social
  • Simone Cortesi
  • Laura Criscuolo
  • Antonio D'Argenio, Nadir
  • Margherita Di Leo
  • Alessio Di Lorenzo
  • Gianfranco Di Pietro, Geofunction
  • Antonio Falciano
  • Sergio Farruggia, Stati Generali dell'Innovazione, AMFM GIS Italia
  • Daniela Ferrari
  • Maurizio Foderà, Kartoblog
  • Marco Fratoddi, Stati generali dell'innovazione
  • Antonio Fregoli, MNDAssociation
  • Pietro Blu Giandonato
  • Cesare Gerbino
  • Simone Giannecchini
  • Jacopo Grazzini
  • Nicola Guarino, ISTC-CNR
  • Giuseppe Iacono, Stati generali dell'innovazione
  • Carlo Infante, Stati generali dell'innovazione, Urban Experience
  • Andrea Latino, Stati generali dell'innovazione
  • Simone Lella
  • Walter Lorenzetti, gis3w
  • Jody Marca
  • Flavia Marzano, Stati Generali dell'Innovazione e Rete WISTER
  • Giacomo Martirano, Epsilon Italia, coordinatore progetto smeSpire
  • Stefania Morrone, Epsilon Italia
  • Lorenzo Orlando, Stati generali dell'innovazione
  • Alessandro Oggioni
  • Mariella Pappalepore, Planetek Italia
  • Stefano Parodi, GeoWebLog
  • Lorenzo Perone
  • Emma Pietrafesa, Stati generali dell'innovazione (Rete WISTER)
  • Renzo Provedel, Stati generali dell'innovazione, SOSLOG
  • Angelo Quaglia
  • Morena Ragone, Stati generali dell'innovazione
  • Paolo Russo, Stati generali dell'innovazione
  • Alessandro Sarretta
  • Monica Sebillo, AMFM GIS Italia
  • Gian Bartolomeo Siletto
  • Claudia Spinnato, Consorzio TICONZERO
  • Lorenzino Vaccari, Provincia Autonoma Trento
  • Franco Vico, AMFM GIS Italia
  • Fabio Vinci, Epsilon Italia
  • Massimo Zotti

Se volete essere aggiunti in questa lista, inserite il vostro nome e/o la vostra affiliazione come commento.



[1] Inoltre, il Decreto 10 novembre 2011 relativo alle regole tecniche del RNDT, emanato dal Ministro per la Pubblica Amministrazione e l'Innovazione e dal Ministero per l’Ambiente e la Tutela del Territorio e del Mare, dispone che il RNDT, parte integrante dell'infrastruttura nazionale, eroghi i servizi di ricerca (art. 2) e prevede la pubblicazione dei metadati nel RNDT, assicurando il rispetto degli adempimenti di cui al Regolamento (CE) n. 1205/2008 e al D. Lgs. n. 32/2010 (DM art. 4)

PostgreSQL: Comparing sets of (similar) street names

Have you ever had the need to join tables by a street name string? Maybe to assign a street code from one table to another. A common problem here is that names could be written in different ways when they actually refer to the same street, for example Via Giacomo Leopardi should match Via Leopardi G.Queen Victoria Street should match Victoria St. and so on...

This join can hardly be performed with EQUAL, LIKE or REGEXP operators. How can we overcome this unequality problem for a proper join?

I made up a short query that did the trick for me (*).

First, you have to import the Trigram module into your database. The Trigram module offers string similarity operators and functions. Login as db superuser on the local host and insert the module like this:

  • PostgreSQL 9.1+
    you might locate and execute pg_env script before proceeding

    $ psql -U postgres -h localhost mydb
    mydb=# CREATE EXTENSION pg_trgm;
    mydb=# \q

  • PostgreSQL 8+
    you might locate and execute pg_env script before proceeding

    $ pg_config --sharedir
    /usr/share/postgresql/8.4
    $ cd /usr/share/postgresql/8.4/contrib
    $ psql -U postgres -h localhost mydb <pg_trgm.sql

Now, consider the following example as an explanation.

Table 1: list of street names to be coded
select * from table1;


Table 2: official street name list and street codes
select * from table2;

To get the right comparison I simply executed the next query. The underlying reasoning is this: for each street name in table 1, scan all the street name items in table 2 and suggest the best fitting in terms of similarity, meaning the first ranked in the trigram similarity function. 

select * from
   (select *, rank() over (partition by st_name_table1 order by similarity desc) from
      (select *, similarity(st_name_table1, st_name_table2) from
          table1, table2
      ) q0
   ) q1
where rank = 1
order by similarity desc


And this is the result (actually I removed from the real output the exactly identical pairs).
You get a table with all the entries of table 1 (non coded) each one related to the best fitting in table 2 (coded streets), just like a left join.
  • green: successfully matched (astonishing!)
  • orange: matched, but you might be careful
  • red: table 1 street names don't get a valid companion, should remain uncoded

Caution! This is suitable for quite small data tables (approximately <1000 rows). Otherwise, you have to set up column indexes and try other ways to get faster performance.

(*) these notes are intended for a PostgreSQL EnterpriseDB regular install or a linux distribution install

Python: shapefile merger utility

This Python script is a commandline shapefile merger. Many options provided to perform quickly the desired merge.
Features:
- accepts wildcards
- recursive (walk through paths)
- list of shapefiles to exclude
- adds fields filled with the source filename and the source feature id (possibility to reuse existing)
- force output result to bi-dimensional geometry type
- dbf schemes union or intersection

Download shapemerger.py or copy the code below.

To run the script you need Python with Ogr/Gdal Python bindings.
On Windows, I suggest to install it from OsGeo4W setup:
- run the setup and choose Advanced Install
- step forward until Select Packages form
- from Commandline_Utilities select and install: gdal, python (meta package), shapelib
- from Libs select and install: gdal-python
- click next and finish installation

Open OsGeo4W shell, move where you downloaded shapemerger.py and type

python shapemerger.py

If everything is ok, it will return the command usage.

Examples:
  • python shapemerger.py -o test.shp *.shp
    - merge all the shapefiles of the current directory in test.shp
    - the reference geometry type is extracted from the first file.
    - union of dbf schemes (default)
  • python shapemerger.py -o test.shp *.shp -i -r -x *_railways.shp
    - merge all the shapefiles of the current directory and subdirectories in test.shp
    - exclude all the shapefiles named as ..._railways.shp
    - intersection of dbf schemes
  • python shapemerger.py -s SHAPENAME -f SHAPEFID -e zones_*.shp
    - merge all the zones_...shp in shapemerge.shp (default name of merge file)
    - rename source filenames field to: SHAPENAME, source feature ids field to: SHAPEFID
    - use existing values of fields SHAPENAME and SHAPEFID if found in the input shapefiles


#!/usr/bin/python
# shapemerger.py
# v.1.1 by toni @ http://furiousgis.blogspot.com
#
# shapefile merger utility 

import os, sys, glob
try:
    from osgeo import ogr, gdal
except:
    try:
        import ogr, gdal
    except:
        print 'Error: shapemerge.py cannot find python OGR and GDAL modules'
        print '       for handling geometries. OsGeo4W python is recommended.'
        sys.exit(1)

def Usage():
    print
    print 'Shapefile merger'
    print 'Usage: shapemerger.py\t[-o <output.shp>] [-n] [-i] [-s <fieldname>] [-f <fieldname>]'
    print '\t\t\t[-e] [-2] [-r] <shapefiles> [-x <exclude shapefile list>]'
    #print '\t\t\t[-e] [-2] [-r] [-t plg|arc|pt|plg2D|arc2D|pt2D|plg3D|arc3D|pt3D]'
    #print '\t\t\t<shapefiles> [-x <exclude shapefile list>]'
    print 
    print '  <shapefiles>: input shapefile list (accepts wildcards)'
    print
    print '  -o:\toutput shapefile name (default: shapemerge.shp)'
    print '  -n:\tdo not add fields for the source filenames and fids'
    print '  -i:\tdbf schemes intersection (default: union)'
    print '  -s:\tcustom name for the source filenames field (default: SOURCESHP)'
    print '  -f:\tcustom name for the source fids field (default: SOURCEFID)'
    print '  -e:\tuse existing source filenames and fids fields, if any'
    print '  -2:\tforce output geometry type to 2D'
    print '  -x:\tlist of shapefiles to exclude (accepts wildcards)'
    print '  -r:\trecursive (grabs every path of input files)'
    #print '  -t:\tfilter input shapefiles by geometry type'
    sys.exit(1)

strCurrentPath = os.path.abspath('')
os.chdir(strCurrentPath)

strOutDirname = ''
strMergeName = 'shapemerge'
strOutBasename = 'shapemerge.shp'

pOgrShpDriver = ogr.GetDriverByName('Esri Shapefile')

lsFilenames=[]
lsExcludeFilenames=[]
lsVerifiedFilenames=[]
lsMergeFieldsToRetain=[]
lsMergeFieldNames=[]
lsMergeFieldOgrDefns=[]
#lsGeomTypeFilters=['PLG','ARC','PT','PLG2D','ARC2D','PT2D','PLG3D','ARC3D','PT3D']

bFlagExclude = 0
bFlagAddSourcefields = 1
bFlagUseExistingSfields = 0
bFlagIntersection = 0
bFlagRecursive = 0
bFlagGeomTFilter = 0
bFlag2D = 0

bGeomTypeRefIsSet = None

strFilenameFieldname = 'SOURCESHP'
strFidFieldname = 'SOURCEFID'

i = 1

# parse commandline (first check for recursive flag)
while i < len(sys.argv): 
    arg = sys.argv[i]
    if arg.upper() == '-R':
        bFlagRecursive = 1
        del sys.argv[i]
        break
    i = i + 1

i = 1
    
while i < len(sys.argv): 
    arg = sys.argv[i]

    if arg.upper() == '-N':
        bFlagAddSourcefields = 0
        bFlagExclude = 0

    elif arg.upper() == '-I':
        try:
            strGdalVersion = int(gdal.VersionInfo()[:2])
        except:
            strGdalVersion = 0
        
        if (strGdalVersion < 19):
            print 'Warning: dbf schemes intersection supported only with gdal >= 1.9.0, ignoring.'
            bFlagIntersection = 0
        else:
            bFlagIntersection = 1
        
        bFlagExclude = 0

    elif arg.upper() == '-S':
        i = i + 1
        strFilenameFieldname = ''.join(sys.argv[i].split())[:10].upper() #upper, tronca al decimo, toglie spazi
        bFlagExclude = 0

    elif arg.upper() == '-F':
        i = i + 1
        strFidFieldname = ''.join(sys.argv[i].split())[:10].upper() #upper, tronca al decimo, toglie spazi
        bFlagExclude = 0
        
    elif arg.upper() == '-E':
        bFlagUseExistingSfields = 1
        bFlagExclude = 0

    elif arg.upper() == '-2':
        bFlag2D = 1
        bFlagExclude = 0

    #if arg.upper() == '-T':
    #    i = i + 1
    #    strGeomTFilter = sys.argv(i)
    #    
    #    if (not (strGeomTFilter in lsGeomTypeFilters)):
    #        print "Warning: invalid geometry type declared, ignoring." 
    #    else:
    #        bFlagGeomTFilter = 1
    #    
    #    bFlagExclude = 0
        
    elif arg.upper() == '-X':
        bFlagExclude = 1

    elif arg.upper() == '-O':
        i = i + 1
        bFlagExclude = 0

        strOutDirname, strOutBasename = os.path.split(sys.argv[i])
        strMergeName, strExtension = os.path.splitext(strOutBasename)
            
    else:
        # expands wildcards: exclude, input list/recursive input list
        if (bFlagRecursive):
            strInputShapePath, strInputShapeBasename = os.path.split(arg)
            strInputShapePath = strInputShapePath or '.'
            
            lsInputShapeRoots=[]
            lsInputShapeRoots.extend(glob.glob(strInputShapePath)) #expands paths with wildcards

        if (bFlagExclude and bFlagRecursive):
            for strRootPath in lsInputShapeRoots:
                for strRoot,strDirs,strFiles in os.walk(strRootPath):
                    lsExcludeFilenames.extend ( glob.glob(strRoot + os.sep + strInputShapeBasename) )
        
        elif (bFlagExclude and (not bFlagRecursive)):
            lsExcludeFilenames.extend ( glob.glob(arg) )
        
        elif (bFlagRecursive):
            for strRootPath in lsInputShapeRoots:
                for strRoot,strDirs,strFiles in os.walk(strRootPath):
                    lsFilenames.extend ( glob.glob(strRoot + os.sep + strInputShapeBasename) )
        else:
            lsFilenames.extend ( glob.glob(arg) )
                
    i = i + 1

for strItem in lsExcludeFilenames:
    try:
        lsFilenames.remove(strItem)
    except:
        pass

if len(lsFilenames) == 0:
    Usage()
    sys.exit(1)

if not os.path.isdir(strOutDirname or '.'):
    print strOutDirname + ': invalid directory, outputting here'
    strOutDirname = ''

if os.path.isfile( os.path.join(strOutDirname, strOutBasename) ):
    strSuffix = ''
    j = 0
    while os.path.isfile( os.path.join(strOutDirname,strMergeName + strSuffix + '.shp') ):
        j = j + 1
        strSuffix = '_' + str(j)
  
    if (strMergeName != 'shapemerge'):
        print strOutBasename + ': file exists, outputting in ' + strMergeName + strSuffix + '.shp'
  
    if (j > 0):
        strMergeName = strMergeName + strSuffix
        strOutBasename = strMergeName + '.shp'  
  
# geometry type check / "source attributes" renaming
strSuffix = ''
j = 0 
for i,strFile in enumerate(lsFilenames):
    try:
        pOgrDatasource = pOgrShpDriver.Open(strFile)
        pOgrLayer = pOgrDatasource.GetLayer()
        pOgrLayerSchema = pOgrLayer.GetLayerDefn()
        
        if (not bGeomTypeRefIsSet):
            # picks up the first geometry type as reference
            intGeometryTypeRef = pOgrLayer.GetGeomType()
            
            if (bFlag2D):
                # fetch first feature to flatten geometry and grab TypeRef
                pOgrFeature = pOgrLayer.GetNextFeature()
                pOgrFeatureGeometry = pOgrFeature.GetGeometryRef()
                pOgrClonedGeometry = pOgrFeatureGeometry.Clone()
                pOgrClonedGeometry.FlattenTo2D()
                intGeometryTypeRef = pOgrClonedGeometry.GetGeometryType() 
                #pOgrFeature.Destroy()
                #pOgrClonedGeometry
            
            bGeomTypeRefIsSet = i # register index of the reference shapefile
                
        #check if layers the picked geometry type (ignore 3D flag)
        if pOgrLayer.GetGeomType() != intGeometryTypeRef:
            if ((pOgrLayer.GetGeomType() == ogr.wkbPolygon) or (pOgrLayer.GetGeomType() == ogr.wkbPolygon25D)) and ((intGeometryTypeRef == ogr.wkbPolygon) or (intGeometryTypeRef == ogr.wkbPolygon25D)) \
                or ((pOgrLayer.GetGeomType() == ogr.wkbLineString) or (pOgrLayer.GetGeomType() == ogr.wkbLineString25D)) and ((intGeometryTypeRef == ogr.wkbLineString) or (intGeometryTypeRef == ogr.wkbLineString25D)) \
                or ((pOgrLayer.GetGeomType() == ogr.wkbPoint) or (pOgrLayer.GetGeomType() == ogr.wkbPoint25D)) and ((intGeometryTypeRef == ogr.wkbPoint) or (intGeometryTypeRef == ogr.wkbPoint25D)) :
                    pass
            else:
                #print 'Error: shapefiles must have the same geometry type.'
                #print '\t' + os.path.basename(lsFilenames[0]) + ': ' + ogr.GeometryTypeToName(intGeometryTypeRef).upper()
                #print '\t' + os.path.basename(lsFilenames[i]) + ': ' + ogr.GeometryTypeToName(pOgrLayer.GetGeomType()).upper() 
                #sys.exit(1)
                print 'Warning: ' + os.path.basename( strFile ) + ' type ' + ogr.GeometryTypeToName(pOgrLayer.GetGeomType()).upper() + ' unfits first file type ' + ogr.GeometryTypeToName(intGeometryTypeRef).upper() + ' (' + os.path.basename(lsFilenames[bGeomTypeRefIsSet]) +'), skipping.'
                continue
                
        if ( (not bFlagUseExistingSfields) and bFlagAddSourcefields ):
            # check for "source fields" in current shapefile to eventually rename them
            while ( (pOgrLayerSchema.GetFieldIndex(strFilenameFieldname) != -1) or (pOgrLayerSchema.GetFieldIndex(strFidFieldname) != -1) ):
                j = j + 1
                strSuffix = '_' + str(j)
                strFilenameFieldname = strFilenameFieldname[:(10 - len(strSuffix))] + strSuffix
                strFidFieldname      = strFidFieldname[:(10 - len(strSuffix))] + strSuffix

        lsVerifiedFilenames.append( strFile )
        pOgrDatasource.Destroy()
        
    except:
        print 'Warning: ' + strFile + ' invalid shapefile, skipping.'
        #continue

if (len(lsVerifiedFilenames) == 0):
    print 'Error: all the input shapefiles are invalid.'
    sys.exit(1)
else:
    lsFilenames = lsVerifiedFilenames

if (j >0):
    print "Warning: \"source fieldnames\" renamed to: " +  strFilenameFieldname + ", " + strFidFieldname

# building attributes
if (bFlagAddSourcefields):
    # "source fields" attributes
    lsMergeFieldOgrDefns.append( ogr.FieldDefn(strFilenameFieldname, ogr.OFTString) )
    lsMergeFieldOgrDefns[0].SetWidth(255)
    lsMergeFieldOgrDefns.append( ogr.FieldDefn(strFidFieldname, ogr.OFTInteger) )
    
    lsMergeFieldNames.append(strFilenameFieldname)
    lsMergeFieldNames.append(strFidFieldname)


for i,strFile in enumerate(lsFilenames):
    # shapefiles attributes
    
    if ( i>0 and bFlagIntersection and ( (bFlagAddSourcefields and (len(lsMergeFieldNames)==2)) or (len(lsMergeFieldNames)==0) ) ):
        # no more scheme intersection to do, exit loop 
        break
    
    strFullname = strFile
    strDirname, strBasename = os.path.split(strFullname)
    strShapename, strExtension = os.path.splitext(strBasename)
    
    pOgrShpDatasource = pOgrShpDriver.Open (strFullname)
    pOgrShpLayer = pOgrShpDatasource.GetLayer()
    pOgrShpLayerSchema = pOgrShpLayer.GetLayerDefn()
     
    for j in range( pOgrShpLayerSchema.GetFieldCount() ):
        # parse shapefile attributes
        pOgrShpFielddefn = pOgrShpLayerSchema.GetFieldDefn(j)
        strShpFieldname = pOgrShpFielddefn.GetNameRef().upper()
        
        try:
            iMergeListsFieldIdx = lsMergeFieldNames.index(strShpFieldname)
        except:
            iMergeListsFieldIdx = -1

        if (iMergeListsFieldIdx == -1):
            
            if ( (not bFlagIntersection) or (bFlagIntersection and (i == 0)) ):
                # field doesn't exist: create it. If flag intersection is on, parse only the first file
                pOgrFielddefn = ogr.FieldDefn( strShpFieldname, pOgrShpFielddefn.GetType() )
                pOgrFielddefn.SetWidth( pOgrShpFielddefn.GetWidth() )
                pOgrFielddefn.SetPrecision( pOgrShpFielddefn.GetPrecision() )                
                pOgrFielddefn.SetJustify( pOgrShpFielddefn.GetJustify() )                                
                
                lsMergeFieldOgrDefns.append( pOgrFielddefn )
                lsMergeFieldNames.append( strShpFieldname )
                
                # add field to the retain list for intersection
                #lsMergeFieldsToRetain.append( lsMergeFieldNames.index(strShpFieldname) )
                lsMergeFieldsToRetain.append( len(lsMergeFieldNames) - 1 )
        
        else:
            # field exist, check if properties should be enlarged
            pOgrFielddefn = lsMergeFieldOgrDefns[iMergeListsFieldIdx]
            
            if pOgrFielddefn.GetWidth() < pOgrShpFielddefn.GetWidth(): 
                pOgrFielddefn.SetWidth( pOgrShpFielddefn.GetWidth() )
            
            elif pOgrFielddefn.GetPrecision < pOgrShpFielddefn.GetPrecision:
                pOgrFielddefn.SetPrecision( pOgrShpFielddefn.GetPrecision() )
                                    
            if (bFlagIntersection):
                # add field to the retain list for intersection
                lsMergeFieldsToRetain.append(iMergeListsFieldIdx)
            
    if ( bFlagIntersection ):
        # performing dbf intersection
        
        if (bFlagAddSourcefields):
            lsMergeFieldsToRetain.append( lsMergeFieldNames.index(strFilenameFieldname) )
            lsMergeFieldsToRetain.append( lsMergeFieldNames.index(strFidFieldname) )
            
        lsMergeFieldsToDelete = range( len(lsMergeFieldNames)-1, -1, -1) # reverse: from last index to 0
        
        for iIdx in lsMergeFieldsToRetain:
            try:
                lsMergeFieldsToDelete.remove(iIdx)
            except:
                pass
       
        for iIdx in lsMergeFieldsToDelete:
            # add or remove fields items
            del lsMergeFieldOgrDefns[iIdx]
            del lsMergeFieldNames[iIdx]            

    pOgrShpDatasource.Destroy()
    lsMergeFieldsToRetain = []

# generate the merge layer
strMergeShpFullname = os.path.join(strOutDirname, strOutBasename)

pOgrMergeShpDatasource = pOgrShpDriver.CreateDataSource( strMergeShpFullname )
pOgrMergeShpLayer = pOgrMergeShpDatasource.CreateLayer( strMergeName, geom_type = intGeometryTypeRef )

for i in lsMergeFieldOgrDefns:
    # create the fields
    pOgrMergeShpLayer.CreateField(i,1)

pOgrMergeShpLayerSchema = pOgrMergeShpLayer.GetLayerDefn()

# Loop for populating the merge shapefile
for i,strFile in enumerate(lsFilenames):
    strFullname = strFile
    strDirname, strBasename = os.path.split(strFullname)
    strShapename, strExtension = os.path.splitext(strBasename)
    
    pOgrShpDatasource = pOgrShpDriver.Open (strFullname)
    pOgrShpLayer = pOgrShpDatasource.GetLayer()
    pOgrShpLayerSchema = pOgrShpLayer.GetLayerDefn()
            
    pOgrShpLayer.ResetReading()
    
    pOgrShpFeature = pOgrShpLayer.GetNextFeature()    
    
    bFlagExistsFilenameFld = pOgrShpLayerSchema.GetFieldIndex(strFilenameFieldname)
    bFlagExistsFidFld = pOgrShpLayerSchema.GetFieldIndex(strFidFieldname)
    
    while pOgrShpFeature:
        pOgrMergeShpfeature = ogr.Feature(pOgrMergeShpLayerSchema)
        pOgrMergeShpfeature.SetFrom(pOgrShpFeature,1)
        
        if (bFlag2D):
            # performing 2d conversion
            pOgrMergeGeometry2d = pOgrMergeShpfeature.GetGeometryRef().Clone()
            pOgrMergeGeometry2d.FlattenTo2D()
            pOgrMergeShpfeature.SetGeometryDirectly(pOgrMergeGeometry2d)
        
        if (bFlagAddSourcefields):
            # "source fields" values, if any and if they should be populated
            if ( (bFlagUseExistingSfields and bFlagExistsFilenameFld == -1) or (not bFlagUseExistingSfields) ):
                #pOgrMergeShpfeature.SetField( strFilenameFieldname, strBasename )
                pOgrMergeShpfeature.SetField( strFilenameFieldname, strFullname )
            if ( (bFlagUseExistingSfields and bFlagExistsFidFld == -1) or (not bFlagUseExistingSfields) ):
                pOgrMergeShpfeature.SetField( strFidFieldname, pOgrShpFeature.GetFID() )
            
        pOgrMergeShpLayer.CreateFeature(pOgrMergeShpfeature)
        
        pOgrMergeShpfeature.Destroy()
        pOgrShpFeature.Destroy()
        pOgrShpFeature = pOgrShpLayer.GetNextFeature()
            
    pOgrShpDatasource.Destroy()

# flush all to disk before quitting
pOgrMergeShpDatasource.SyncToDisk()

print
print strMergeShpFullname
sys.exit(0)

Qgis: connecting a Geomedia mdb warehouse

...an open-source driver for Geomedia has been released for gvSIG. It's the extMDB extension and works pretty well. Anyway, let's talk about Qgis...

Some tricky ways to display the content of a Geomedia Access Mdb in a Windows installed Qgis.
First of all, you need Gdal driver to be >= 1.9.0 release, so if you have an old install of Qgis, please download Osgeo4W setup and install or update Qgis.

When ready,
  1. Try to load your data
    Add a regular Vector Layer and leave File as source type (do not switch to database), browse for your Geomedia mdb file and select it.

    If everything is ok, Qgis will ask you which layers to load and their spatial reference, then you'll get the feature classes on the map.

If Qgis fails to access to your file, follow next steps.
This implies to modify the metadata structure, so backup your file first. Unfortunately, you're going to see most of (but not all) the database tables.
  1. First workaround
    Open Geomedia mdb with Microsoft Access
    • delete feature classes with no records (pay attention not to delete the Geomedia metadata tables)
    • delete those records in table GeometryProperties where GeometryType = 5
    • delete those records in table GFeatures where PrimaryGeometryFieldName is empty

Not enough? Go forward
  1. Second workaround
    Again, Open Geomedia mdb with Microsoft Access and perform those Sql queries in this order, deleting all the resulting records
    • SELECT GFeatures.* from
      GFeatures WHERE GFeatures.FeatureName in
      (SELECT GFeatures.FeatureName
      FROM (SELECT FieldLookup.FeatureName
      FROM FieldLookup INNER JOIN GeometryProperties ON FieldLookup.IndexID = GeometryProperties.IndexID
      GROUP BY FeatureName HAVING Count(FeatureName)>1) AS QUERY INNER JOIN GFeatures ON QUERY.FeatureName = GFeatures.FeatureName);


    • Select GeometryProperties.* from
      GeometryProperties WHERE GeometryProperties.IndexID IN
      (SELECT GeometryProperties.IndexID
      FROM ((SELECT FieldLookup.FeatureName
      FROM FieldLookup INNER JOIN GeometryProperties ON FieldLookup.IndexID = GeometryProperties.IndexID
      GROUP BY FeatureName HAVING Count(FeatureName)>1) AS QUERY INNER JOIN FieldLookup ON QUERY.FeatureName = FieldLookup.FeatureName) INNER JOIN GeometryProperties ON FieldLookup.IndexID = GeometryProperties.IndexID);


    • SELECT GeometryProperties.*
      FROM (SELECT FieldLookup.*, GFeatures.FeatureName
      FROM FieldLookup LEFT JOIN GFeatures ON FieldLookup.FeatureName = GFeatures.FeatureName
      WHERE (((GFeatures.FeatureName) Is Null))) AS Query1 INNER JOIN GeometryProperties ON Query1.IndexID = GeometryProperties.IndexID;


Now, you should be out of troubles.

Autocad: (re)Create closed boundaries for all hatches

26.03.2013 - updated release, completely rewritten, major bugfixing

This Autolisp routine is intended primarily for a Dwg to Gis conversion. See if it fits your needs.

It recreates boundaries around all the hatch patterns in a drawing at once. The recreated boundaries are drawn in new layers and consist only of single closed polylines (interpreted as polygons by Gis applications). Each curved element of the boundary is converted to a segmented polyline and finally all these pieces are combined into one.
The outcome might be not identical to the original hatch outline, since it has to be considered as an approximation. Anyway, this is exactly the same process used by common desktop Gis when importing drawings.

Tested on Autocad 2012.

Features:
  • creates single closed polyline boundaries (Gis polygons)
  • assigns the hatch color to the boundary
  • draws in new layers called as the source layer of the hatch, plus a custom suffix
  • deletes or keeps all the previous associated boundaries
  • handles every kind of known geometry in a boundary (elliptic arcs, circles...)

Download bulk-recreate-hatch-boundaries.lsp

In Autocad, load it from Tools menu > AutoLISP > Load application
Usage: type BULKRHB in the command-line


Example (click on the images to enlarge)

Hatch on layer "0" (without boundary)


All the hatch outline pieces (HATCHGENERATEBOUNDARY command result)


(BULKRHB performed: 3 closed polylines on layer "0_BOUNDARIES")


The difference between HATCHGENERATEBOUNDARY and BULKRHB in ArcMAP

Geomedia: Useful functional attributes

Useful functional attributes expressions.
Arguments separator could be different according to your international settings.

  • Distance (projected) from point to point
    There might be a bug in the DISTANCE function in case of coincident points, so you better workaround like this:
  • LENGTH( CREATEPOLYLINE( pointgeometry_1 ; pointgeometry_2 ) ; 1 )

  • Add leading zeros
    examples to get a four-digit string like 0001, 0012, 0123
    • to a numeric field:
      TEXT( numericfield ; "0000" )
    • to a text field all filled with numeric strings:
      TEXT( VALUE( textfield ) ; "0000" )
    • to a text field:
      REPT( "0"; 4 - LEN( textfield ) ) + textfield

  • Trim leading zeros
    MID( textfield ; SEARCH( "[^0]" ; textfield ) )

  • Point projected coordinates
    • north:
      X( geometryfield ; 1 )
    • east:
      Y( geometryfield ; 1 )

  • Convert text geometry to point geometry
    ORIGIN( textgeometryfield )

  • In analytical merge, count those records that fulfill a specified where condition
    example: count all the records where field is less than 20
    SUM( IF( numericfield < 20 ; 1 ; 0) )

  • Casting a numeric value to text
    TEXT( textfield ; "" )
    Fundamentals of text formatting: put a "0" to get a digit-value or a 0; put a "#" to get a digit-value or nothing; add a section in the mask for negative values. Examples:
    TEXT( 2.715 ; "0.00"  ) returns "2.72"
    TEXT( 2.71  ; "0.000" ) returns "2.710"
    TEXT( 2.71  ; "0.###" ) returns "2.71"
    TEXT( -2.7  ; "0.###;-0.000" ) returns "-2.700"


  • Rotate a point or a label by an attribute value in degrees
    SPIN( {text}pointgeometryfield ; numericfield )

  • Create a linestring from rows of points aggregated by an attribute value
    Use Analytical merge for this. Select the By Attribute option and choose the discriminating attribute in the list below. Then edit the Geometry attribute in the Output functional attributes box clicking on Properties. Delete the default MERGE(Input.Geometry) and type
    CREATEPOLYLINE( pointgeometryfield ; numericfield )
    numericfield is an optional but recommended attribute to set the order of the point sequence