Vecteurs

Lire un fichier vecteur

La lecture d’un fichier vecteur se fait en plusieurs étapes :

• ouverture du fichier avec le pilote, on obtient l’objet initial ;
• lecture de la couche de cet objet, on obtient l’objet layer ;
• lecture des objets (feature) de l’objet layer, on obtient l’objet feature.

Pour chaque objet on a un ensemble de méthodes liées à la classe à laquelle il appartient. Les méthodes peuvent être listées grâce à la méthode dir(objet). Exemple dir(ogr) mais vous pouvez aussi regarder la doc de l’API Python.

Un fichier vecteur est lu grâce à la méthode Open() de la classe OGRSFDriverRegistrar, exemple :

ds_ref = ogr.Open("/url/to/file.shp", update = 0)

Nous obtenons un objet initial, quelles sont les méthodes que nous pouvons lui appliquer ? Facile : dir(ds_ref) ou http://wiki.gloobe.org/gdal/classosgeo_1_1ogr_1_1DataSource.html

Vous devez ensuite récupérer la ou les couche(s) grâce à GetLayer(), pour un shape, le nom de la couche est la même que le nom du fichier sans l’extension sinon vous devez utiliser le numéro de la couche :

Layer_ref = ds_ref.GetLayer()

La variable Layer_ref contient maintenant l’objet layer. Si vous ne connaissez pas le numéro de la couche, vous pouvez utiliser la méthode GetLayerByName() pour récupérer l’objet layer.

Vous pouvez ainsi accéder aux informations de cette couche :

Layer_ref.GetName() #Ne devrait pas être possible, GetName() n'existe pas pour l'objet OGRLayer !
Layer_ref.GetFeatureCount()
feature = Layer_ref.GetFeature(1) # Récupère l'objet (feature) numéro 1

Vous pouvez trouvez plus d’information dans la classe OGRDataSource.

Lorsque vous avez l’objet feature, vous pouvez gérer vos objets géométriques avec la classe Feature.

import os, sys, string,  shutil
from osgeo import ogr

#Nom et géométrie d'un shapefile
ds_ref = ogr.Open('D:/monShapeFile.shp', update = 0)
Layer_ref = ds_ref.GetLayer()
print "Layer name : " + Layer_ref.GetName()
print "Number of features : " + str(Layer_ref.GetFeatureCount())

#Sélection sur un shapefile
query='SURFACE > 1000'
Layer_ref.SetAttributeFilter(query)

# Export de la sélection
driver = ogr.GetDriverByName('ESRI Shapefile')
ds_dest=driver.CreateDataSource('D:/maSelection.shp')
ds_dest.CopyLayer (Layer_ref ,  'maSelection')
ds_dest.Destroy()

ds_ref.Destroy()

Écrire un fichier vecteur

# Create a Shapefile
output="d:/temp/testogr.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(output):
    driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
layer = ds.CreateLayer(output, geom_type=ogr.wkbPolygon)

Création de shapefile à partir de coordonnées

Pour créer des objets géométriques à partir de coordonnées une possibilité est de passer par Wkt ([http://en.wikipedia.org/wiki/Well-known_text Well-known text]). Il s’agit d’une simple chaîne de caractère qui permet de définir le type d’objet et les coordonnées associées.

Exemple de Wkt :

wkt_point = 'POINT(6 10)'
wkt_ligne = 'LINESTRING(3 4,10 50,20 25)'
wkt_polygone = 'POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2, 3 2, 3 3, 2 3,2 2))'

Vous pouvez alors créer un objet Geometry en utilisant la méthode CreateGeometryFromWkt.

Exemple : création d’un shapefile de points à partir des coordonnées x y

#Création d'un shapefile de points
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource('D:/monShapefile.shp')
layer = ds.CreateLayer('monShapefile', geom_type=ogr.wkbPoint)

#
feature_def=layer.GetLayerDefn()
f = ogr.Feature(feature_def)

# Ajout du point
x = 731065
y = 2368493
wkt = 'POINT(%f %f)' % (x, y)
p = ogr.CreateGeometryFromWkt(wkt)
f.SetGeometryDirectly(p)
layer.CreateFeature(f)

f.Destroy()

Projections

Exemple : définir une projection (Lambert 2 étendu) à un shapefile :

# Crée le système de référence spatiale
lambert2e = osr.SpatialReference()
lambert2e.ImportFromEPSG(27572)

#Crée le fichier .proj
lambert2e.MorphToESRI()
wkt_proj = lambert2e.ExportToWkt()
prj_file = open('D:/monShapefile.prj', 'w')
prj_file.write(wkt_proj)
prj_file.close()

Les manipulations de systèmes de projections se font sur les objets Geometry. On peut soit utiliser la méthode TransformTo soit Transform.

• TransformTo suppose qu’un système de coordonnées a été défini pour l’objet Geometry et qu’un objet SpatialReference a été défini comme système de coordonnées d’arrivé.
• Transform suppose qu’un objet CoordinateTransformation ait été défini. Il n’est pas indispensable d’avoir défini le système de coordonnées de l’objet Geometry traité, il sera supposé que ce système correspond à celui de départ de l’objet CoordinateTransformation.

Il est conseillé d’utiliser plutôt Transform et donc de définir au préalable un objet CoordinateTransformation pour des modifications sur un grand nombre d’objets.

Comment Reprojeter des données

On créé deux objets Références Spatiales et pour chacun on importe les données de reprojection via un import à partir d’un code EPSG. Puis on créé un objet Géométrie à partir d’un WKT. On lui assigne un projection, puis on reprojète l’objet.

to_srs = osr.SpatialReference()
to_srs.ImportFromEPSG(4326)
from_srs = osr.SpatialReference()
from_srs.ImportFromEPSG(900913)

wkt = 'POINT(%f %f)' % (x, y)

pt = ogr.CreateGeometryFromWkt(wkt)
pt.AssignSpatialReference(from_srs)
pt.TransformTo(to_srs)
geom = pt.GetX(),pt.GetY()
return pt

Requêtes spatiales

Par bbox

Par objet géométrique

Par requête SQL