Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Manipulation du fichier de formes (vecteur)

This section introduces geospatial analysis of shapefiles with gdal, ogr, and osr. For interactive reading and executing code blocks Binder and find geo01-shp.ipynb, or install Python and JupyterLab locally.

Watch this section and the Python tutorials in video formats

Watch this section as a video on the @Hydro-Morphodynamics channel on YouTube.

Charger un fichier de formes existant

Le module d’OSGeo ogr est une excellente source de manipulation des shapefiles. Après l’importation des bibliothèques, le bon pilote pour gérer les shapefiles ("ESRI Shapefile") avec Python doit être appelé. Ensuite, avec l’objet pilote (ogr.GetDriverByName("SHAPEFILE")), nous pouvons ouvrir un shapefile (objet shp_driver.Open("SHAPEFILE")), qui contient des informations sur la couche. C’est précisément l’information de la couche (c.-à-d. les références aux attributs shapefile) avec laquelle nous voulons travailler. Par conséquent, nous devons activer un objet shapefile avec shp_dataset.GetLayer().

from osgeo import ogr
shp_driver = ogr.GetDriverByName("ESRI Shapefile")
shp_dataset = shp_driver.Open("geodata/shapefiles/cities.shp")
shp_layer = shp_dataset.GetLayer()

Créer un nouveau fichier de formes

ogr permet également de créer un nouveau point, ligne, ou polygone shapefile. Le bloc de code suivant définit une fonction pour créer un shapefile, où l’argument optionnel overwrite est utilisé pour contrôler si un shapefile existant avec le même nom doit être écrasé (par défaut: True). La commande shp_driver.CreateDataSource(SHP-FILE-DIR) crée un nouveau shapefile et le reste de la fonction ajoute un calque au shapefile si les paramètres optionnels de mot-clé layer_name et layer_type sont fournis. Les deux mots clés optionnels doivent être strings, où layer_name peut être n’importe quel nom pour le nouveau calque. layer_type doit être soit "point", "line", ou "polygon" pour créer un point, (poly)line, ou polygone shapefile, respectivement. La fonction utilise le geometry_dict dictionnaire pour attribuer le correct ogr.SHP-TYPE au calque. Il existe d’autres options pour étendre la fonction create_shp(...) listée sur pcjericks’ Github pages.

from osgeo import ogr
import os


def create_shp(shp_file_dir, overwrite=True, *args, **kwargs):
    """
    :param shp_file_dir: STR of the (relative shapefile directory (ends on ".shp")
    :param overwrite: [optional] BOOL - if True, existing files are overwritten
    :kwarg layer_name: [optional] STR of the layer_name - if None: no layer will be created (max. 13 chars)
    :kwarg layer_type: [optional] STR ("point, "line", or "polygon") of the layer_name - if None: no layer will be created
    :output: returns an ogr shapefile layer
    """
    shp_driver = ogr.GetDriverByName("ESRI Shapefile")

    # check if output file exists if yes delete it
    if os.path.exists(shp_file_dir) and overwrite:
        shp_driver.DeleteDataSource(shp_file_dir)

    # create and return new shapefile object
    new_shp = shp_driver.CreateDataSource(shp_file_dir)

    # create layer if layer_name and layer_type are provided
    if kwargs.get("layer_name") and kwargs.get("layer_type"):
        # create dictionary of ogr.SHP-TYPES
        geometry_dict = {"point": ogr.wkbPoint,
                         "line": ogr.wkbMultiLineString,
                         "polygon": ogr.wkbMultiPolygon}
        # create layer
        try:
            new_shp.CreateLayer(str(kwargs.get("layer_name")),
                                geom_type=geometry_dict[str(kwargs.get("layer_type").lower())])
        except KeyError:
            print("Error: Invalid layer_type provided (must be 'point', 'line', or 'polygon').")
        except TypeError:
            print("Error: layer_name and layer_type must be string.")
        except AttributeError:
            print("Error: Cannot access layer - opened in other program?")
    return new_shp

The create_shp function is also provided with flusstools (flusstools.geotools.create_shp()) and aids to create a new shapefile (make sure to get the directory right):

a_new_shp_file = create_shp(r"" + os.getcwd() + "/geodata/shapefiles/new_polygons.shp", layer_name="basemap", layer_type="polygon")

# release data source
a_new_shape_file = None

Shapefiles peut également être créé et dessiné dans QGIS et les figures suivantes guident à travers la procédure de création d’un fichier de forme polygone. Nous n’aurons plus besoin du shapefile résultant dans cette section, mais plus tard pour le faire interagir avec les ensembles de données raster.

The first step to making a shapefile with QGIS is obviously to run QGIS and create a new project. The following example uses water depth and flow velocity raster data as background information to delineate the so-called morphological unit of slackwater. Both the water depth and flow velocity rasters are part of the River Architect sample datasets (precisely located at RiverArchitect/SampleData/01_Conditions/2100_sample/). After downloading the sample data, they can be imported in QGIS by dragging the files from the Browser panel into the Layers panel. Then:

qgis create new shapefile

Figure 1:QGIS : Créer un calque > Nouveau calque...

qgis define new shapefile

Figure 2:Définir le nouveau calque Shapefile

qgis shapefile toggle-editing

Figure 3:Activer l’édition de shapefile

qgis draw polygon shapefile

Figure 4:Dessiner les polygones dans un fichier de forme en QGIS.

Obtenir et définir les projections de Shapefile

La terminologie utilisée dans les fichiers .prj d’un shapefile correspond aux définitions de la section geospatial data. En Python, des informations sur le système de coordonnées sont disponibles via shp_layer.GetSpatialRef() de la bibliothèque ogr:

shp_srs = shp_layer.GetSpatialRef()
print(shp_srs)
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]]

Cette définition GEOGCS du fichier de forme ci-dessus correspond au texte bien connu d’Esri (WKT). Comme le format shapefile a été développé par Esri, le format WKT d’Esri (esriwkt) doit être utilisé dans les fichiers .prj. Le Open Geospatial Consortium (OGC) utilise un texte bien connu dans ses définitions EPSG:XXXX (par exemple, disponible à spatialreference.org).

GEOGCS["WGS 84",
       DATUM["WGS_1984", SPHEROID["WGS84", 6378137, 298.257223563, AUTHORITY["EPSG", "7030"]], AUTHORITY["EPSG","6326"]],
       PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
       UNIT["degree",0.01745329251994328, AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

Pour redéfinir ou définir le système de coordonnées d’un shapefile, nous pouvons utiliser spatialreference.org par défaut urllib bibliothèque.

import urllib

# function to get spatialreferences with epsg code
def get_esriwkt(epsg):    
    # usage get_epsg_code(4326)
    try:
        with urllib.request.urlopen("https://spatialreference.org/ref/epsg/{0}/esriwkt/".format(epsg)) as response:
            return str(response.read()).strip("b").strip("'")
    except Exception:
        pass
    try:
        with urllib.request.urlopen("https://spatialreference.org/ref/sr-org/epsg{0}-wgs84-web-mercator-auxiliary-sphere/esriwkt/".format(epsg)) as response:
            return str(response.read()).strip("b").strip("'")
        # sr-org codes are available at "https://spatialreference.org/ref/sr-org/{0}/esriwkt/".format(epsg)
        # for example EPSG:3857 = SR-ORG:6864 -> https://spatialreference.org/ref/sr-org/6864/esriwkt/ = EPSG:3857
    except Exception:
        print("ERROR: Could not find epsg code on spatialreference.org. Returning default WKT(epsg=4326).")
        return 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295],UNIT["Meter",1]]'

Cette fonction peut être utilisée pour créer un nouveau fichier de projection :

# open the hypy-area shapefile
shp_file = "hypy-area"

# create new .prj file for the shapefile (.shp and .prj must have the same name)
with open("geodata/shapefiles/{0}.prj".format(shp_file), "w") as prj:
    epsg_code = get_esriwkt(4326)
    prj.write(epsg_code)
    print("Wrote projection file : " + epsg_code)
Wrote projection file : GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Une alternative hors ligne pour générer un fichier .prj est la bibliothèque osr qui vient avec gdal:

from osgeo import osr

def get_wkt(epsg, wkt_format="esriwkt"):
    default = 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295],UNIT["Meter",1]]'
    spatial_ref = osr.SpatialReference()
    try:
        spatial_ref.ImportFromEPSG(epsg)
    except TypeError:
        print("ERROR: epsg must be integer. Returning default WKT(epsg=4326).")
        return default
    except Exception:
        print("ERROR: epsg number does not exist. Returning default WKT(epsg=4326).")
        return default
    if wkt_format=="esriwkt":
        spatial_ref.MorphToESRI()
    # return a nicely formatted WKT string (alternatives: ExportToPCI(), ExportToUSGS(), or ExportToXML())
    return spatial_ref.ExportToPrettyWkt()

Transformer (ré-projeter) un Shapefile

Une nouvelle projection peut être nécessaire, par exemple, pour utiliser un fichier de formes à EPSG:4326 (par exemple, créé avec QGIS) dans une application Web-GIS (par exemple, des cartes de rue ouvertes) qui utilise habituellement EPSG:3857. Pour appliquer une projection différente aux objets géométriques d’un fichier de forme, il ne suffit pas de réécrire le fichier .prj. Par conséquent, l’exemple suivant montre la re-projection du fichier countries.shp shapefile à partir du Natural Earth Quick start kit. Le workflow suivant effectue la reprojection:

  • Le shapefile à transformer se trouve dans le sous-répertoire geodata/shapefiles/countries.shp et l’ouvre dans le script Python comme décrit ci-dessus.

  • Lire et identifier le système de référence spatiale utilisé dans le fichier de formes d’entrée:

    • Créer un objet de référence spatiale avec in_sr = osr.SpatialReference(str(shapefile.GetSpatialRef())).

    • Détecter le système de référence spatiale au format EPSG avec AutoIdentifyEPSG().

    • Affecter le système de référence spatiale formaté EPSG à l’objet de référence spatiale du fichier de forme d’entrée (ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))).

  • Créer la référence spatiale de sortie avec out_sr = osr.SpatialReference() et appliquer le code EPSG cible (out_sr.ImportFromEPSG(3857)).

  • Créer un objet de transformation de coordonnées (coord_trans = osr.CoordinateTransformation(in_sr, out_sr)) qui permet de re-projecter les objets géométriques plus tard.

  • Create the output shapefile, which will correspond to a copy of the input shapefile (use the above-defined create_shp function with layer_name="basemap" and layer_type="line").

  • Copier les noms de champs et les types du fichier de forme d’entrée :

    • Lire la couche d’attribut à partir des définitions de couche du fichier d’entrée avec in_lyr_def = in_shp_lyr.GetLayerDefn()

    • Itérer à travers les définitions de champ et les ajouter à la couche de shapefile de sortie (out_shp_lyr)

  • Itérer à travers les caractéristiques géométriques dans le fichier de forme d’entrée:

    • Utilisez les définitions de calque du nouveau fichier de formes (output) (out_shp_lyr_def = out_shp_lyr.GetLayerDefn()) pour ajouter les objets de géométrie transformés plus tard.

    • Définir une variable d’itération in_feature comme une instance de in_shp_lyr.GetNextFeature.

    • Dans une boucle while, activez chaque géométrie (geometry = in_feature.GetGeometryRef()) dans le shapefile d’entrée, transformez le geometry (geometry.Transform(coord_trans)), convertissez-la en ogr.Feature() avec la méthode SetGeometry(geometry), copiez les définitions de champs (nested for-loop), et ajoutez la nouvelle fonctionnalité à la couche shapefile de sortie (out_shp_lyr_def.CreateFeature(out_feature)).

    • À la fin de la boucle while-loop, recherchez la prochaine fonctionnalité dans les attributs du fichier de formes d’entrée avec in_feature = in_shp_lyr.GetNextFeature()

  • Débloquez (relevez) tous les calques et les shapefiles en écraseant les objets avec None (rien n’est écrit dans n’importe quel fichier tant que ces variables existent !).

  • Attribuer la nouvelle projection EPSG:3857 en utilisant la fonction get_wkt définie ci-dessus.

from osgeo import ogr
from osgeo import osr

shp_driver = ogr.GetDriverByName("ESRI Shapefile")

# open input shapefile and layer
in_shp = shp_driver.Open(r"" + os.path.abspath("") + "/geodata/shapefiles/countries.shp")
in_shp_lyr = in_shp.GetLayer()

# get input SpatialReference
in_sr = osr.SpatialReference(str(in_shp_lyr.GetSpatialRef()))
# auto-detect epsg
in_sr.AutoIdentifyEPSG()
# assign input SpatialReference
in_sr.ImportFromEPSG(int(in_sr.GetAuthorityCode(None)))

# create SpatialReference for new shapefile
out_sr = osr.SpatialReference()
out_sr.ImportFromEPSG(3857)

# create a CoordinateTransformation object
coord_trans = osr.CoordinateTransformation(in_sr, out_sr)

# create output shapefile and get layer
out_shp = create_shp(r"" + os.path.abspath("") + "/geodata/shapefiles/countries-web.shp", layer_name="basemap", layer_type="line")
out_shp_lyr = out_shp.GetLayer()

# look up layer (features) definitions in input shapefile
in_lyr_def = in_shp_lyr.GetLayerDefn()
# copy field names of input layer attribute table to output layer
for i in range(0, in_lyr_def.GetFieldCount()):
    out_shp_lyr.CreateField(in_lyr_def.GetFieldDefn(i))

# instantiate feature definitions object for output layer (currently empty)
out_shp_lyr_def = out_shp_lyr.GetLayerDefn()

# iteratively append all input features in new projection
in_feature = in_shp_lyr.GetNextFeature()
while in_feature:
    # get the input geometry
    geometry = in_feature.GetGeometryRef()
    # re-project (transform) geometry to new system
    geometry.Transform(coord_trans)
    # create new output feature
    out_feature = ogr.Feature(out_shp_lyr_def)
    # assign in-geometry to output feature and copy field values
    out_feature.SetGeometry(geometry)
    for i in range(0, out_shp_lyr_def.GetFieldCount()):
        out_feature.SetField(out_shp_lyr_def.GetFieldDefn(i).GetNameRef(), in_feature.GetField(i))
    # add the feature to the shapefile
    out_shp_lyr.CreateFeature(out_feature)
    # prepare next iteration
    in_feature = in_shp_lyr.GetNextFeature()

# release shapefiles and layers
in_shp = None
in_shp_lyr = None
out_shp = None
out_shp_lyr = None

# create .prj file for  output shapefile for web application references
with open(r"" + os.path.abspath('') + "/geodata/shapefiles/countries-web.prj", "w+") as prj:
    prj.write(get_wkt(3857))

En cas de succès, la séquence de code in_sr.AutoIdentifyEPSG() retourne 0 (c.-à-d. qu’elle a reconnu avec succès le numéro EPSG), mais malheureusement, de nombreux numéros EPSG ne sont pas connus à AutoIdentifyEPSG(). Dans le cas où AutoIdentifyEPSG() ne fonctionnait pas correctement, la méthode ne retourne pas 0, mais 7, par exemple. Une solution pour la fonctionnalité limitée de srs.AutoIdentifyEPSG() est srs.FindMatches. srs.FindMatches retourne un matching srs_match d’une base de données plus grande, qui est quelque peu imbriquée; par exemple, utilisez:

matches = srs.FindMatches()

Ensuite, matches ressemble à ceci : [(osgeo.osr.SpatialReference, INT)]. Par conséquent, une solution complète pour srs.AutoIdentifyEPSG() (ou in_sr.AutoIdentifyEPSG() dans le bloc de code ci-dessus) ressemble à ceci:

# set epsg and create spatial reference object
epsg = 3857
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)

# identify spatial reference
auto_detect = srs.AutoIdentifyEPSG()
if auto_detect != 0:
    srs = srs.FindMatches()[0][0]  # Find matches returns list of tuple of SpatialReferences
    srs.AutoIdentifyEPSG()  # Re-perform auto-identification

Ajouter des champs et des fonctions de point à un fichier de formes

Une fonction shapefile peut être un point, une ligne, ou un polygone, qui a des attributs de champ (par exemple, "id"=1 pour décrire qu’il s’agit d’un polygone numéro 1 ou associé à un id bloc 1). Les attributs de champ peuvent être plus qu’un simple ID et inclure, par exemple, la zone de polygone ou les étiquettes de ville comme dans l’exemple ci-dessus illustré cités.shp (shp_driver.Open("geodata/shapefiles/cities.shp")).

Pour créer un shapefile, nous pouvons utiliser la fonction create_shp (ou flusstools.geotools.create_shp(shp_file_dir="")) introduite ci-dessus et définir la projection avec la fonction get_wkt() (également introduite ci-dessus). Le bloc de code suivant montre l’utilisation des deux fonctions pour créer un rivière. shp point shapefile qui contient trois points situés sur trois rivières en Europe centrale. Le bloc de code illustre également la création d’un champ dans la table des attributs et la création de trois points. Voici comment le code fonctionne:

  • Le shapefile est situé dans la variable rivers_pts. Notez que le layer_type détermine déjà le type de géométrie qui peut être utilisé dans le shapefile. Par exemple, ajouter une ligne ou une fonction de polygone à une couche ogr.wkbPoint entraînera un message ERROR 1.

  • La variable basemap (couche) est attribuée à la variable lyr = river_pts.GetLayer().

  • Un champ de type string est ajouté et annexé au tableau d’attributs :

    • instantanéez un nouveau champ avec field_gname = ogr.FieldDefn("FIELD-NAME", ogr.OFTString) (rappel : le nom du champ peut ne pas avoir plus de 10 caractères !)

    • ajouter le nouveau champ au shapefile avec lyr.CreateField(field_gname)

    • Les autres types de champ que OFTString peuvent être : OFTInteger, OFTReal, OFTDate, OFTTime, OFTDateTime, OFTBinary, OFTIntegerList, OFTRealList, ou OFTStringList.

  • Ajouter trois points stockés dans pt_names = {RIVER-NAME: (x-coordinate, y-coordinate)} dans une boucle sur les clés du dictionnaire :

    • pour chaque nouveau point, créez une fonctionnalité en tant qu’enfant des définitions de calque avec feature = ogr.Feature(lyr.GetLayerDefn())

    • définir la valeur du nom du champ pour chaque point avec feature.SetField(FIELD-NAME, FIELD-VALUE)

    • créer une chaîne du nouveau point au format WKT avec wkt = "POINT(X-COORDINATE Y-COORDINATE)"

    • convertir le point formaté WKT en une géométrie de type point avec point = ogr.CreateGeometryFromWkt(wkt)

    • définir le nouveau point comme la géométrie de la nouvelle fonctionnalité avec feature.SetGeometry(point)

    • ajouter la nouvelle fonctionnalité au calque avec lyr.CreateFeature(feature)

  • Débloquer (libérer) le shapefile en écraseant la variable lyr et river_pts avec None.

shp_dir = r"" + os.path.abspath('') + "/geodata/shapefiles/rivers.shp"
river_pts = create_shp(shp_dir, layer_name="basemap", layer_type="point")

# create .prj file for the shapefile for web application references
with open(shp_dir.split(".shp")[0] + ".prj", "w+") as prj:
    prj.write(get_wkt(3857))

# get basemap layer
lyr = river_pts.GetLayer()

# add string field "rivername"
field_gname = ogr.FieldDefn("rivername", ogr.OFTString)
lyr.CreateField(field_gname)

# names and coordinates of central EU rivers in EPSG:3857 WG84 / Pseudo-Mercator
pt_names = {"Aare": (916136.03, 6038687.72),
            "Ain": (623554.12, 5829154.69),
            "Inn": (1494878.95, 6183793.83)}

# add the three rivers as points to the basemap layer
for n in pt_names.keys():
    # create Feature as child of the layer
    pt_feature = ogr.Feature(lyr.GetLayerDefn())
    # define value n (river) in the rivername field
    pt_feature.SetField("rivername", n)
    # use WKT format to add a point geometry to the Feature
    wkt = "POINT(%f %f)" % (float(pt_names[n][0]), float(pt_names[n][1]))
    point = ogr.CreateGeometryFromWkt(wkt)
    pt_feature.SetGeometry(point)
    # append the new feature to the basemap layer
    lyr.CreateFeature(pt_feature)
    
# release files
lyr = None
river_pts = None

Le fichier rivers.shp shapefile résultant peut être importé dans QGIS avec un Modèle numérique de terrain (MNT) de la Natural Earth Quick start kit.

river map shapefile

Figure 5:Les noms de la rivière se lient aux points d’un fichier point de forme montré dans QGIS.

Multiligne (Polyline) Shapefile

Comme pour la création et l’ajout de points à un nouveau shapefile, une ligne (multi) (ou polyline) peut être ajoutée à un shapefile. La fonction create_shp crée un shapefile multiligne lorsque le type de calque est défini comme "line" (flusstools.geotools.create_shp(shp_file_dir="", layer_type="line")). Le système de coordonnées est créé avec la fonction get_wkt() définie ci-dessus.

Le bloc de codes suivant utilise les coordonnées des villes le long du Rhin, stockées dans un dictionnaire* nommé station_names. Les noms des villes ne sont pas utilisés et seules les coordonnées sont jointes à line.AddPoint(X, Y). Comme avant, un champ est créé pour donner un nom à la rivière. La fonctionnalité de ligne réelle est à nouveau créée comme un enfant du calque avec line_feature = ogr.Feature(lyr.GetLayerDefn()). L’exécution de ce bloc de codes produit une ligne qui suit approximativement le Rhin entre la France et l’Allemagne.

shp_dir = r"" + os.path.abspath("") + "/geodata/shapefiles/rhine_proxy.shp"
rhine_line = create_shp(shp_dir, layer_name="basemap", layer_type="line")

# create .prj file for the shapefile for web application references
with open(shp_dir.split(".shp")[0] + ".prj", "w+") as prj:
    prj.write(get_wkt(3857))

# get basemap layer
lyr = rhine_line.GetLayer()

# coordinates for EPSG:3857 WG84 / Pseudo-Mercator
station_names = {"Basel": (844361.68, 6035047.42),
                 "Kembs": (835724.27, 6056449.76),
                 "Breisach": (842565.32, 6111140.43),
                 "Rhinau": (857547.04, 6158569.58),
                 "Strasbourg": (868439.31, 6203189.68)}

# create line object and add points from station names
line = ogr.Geometry(ogr.wkbLineString)
for stn in station_names.values():
    line.AddPoint(stn[0], stn[1])

# create field named "river"
field_name = ogr.FieldDefn("river", ogr.OFTString)
lyr.CreateField(field_name)

# create feature, geometry, and field entry
line_feature = ogr.Feature(lyr.GetLayerDefn())
line_feature.SetGeometry(line)
line_feature.SetField("river", "Rhine")

# add feature to layer
lyr.CreateFeature(line_feature)

lyr = None
rhine_line = None

Le fichier rhine_proxy.shp shapefile résultant peut être importé en QGIS avec un Modèle numérique de terrain (MNT) et les villes point shapefile à partir du Natural Earth Quick start kit.

rhine line shapefile

Figure 6:La ligne approximant le Rhin en QGIS.

Polygon Shapefile

Les polygones sont des patchs de surface qui peuvent être créés point par point, ligne par ligne, ou à partir d’une définition "Multipolygon" WKB. Lorsque nous créons des polygones à partir de points ou de lignes, nous voulons créer une surface et c’est pourquoi le type de géométrie correspondant est appelé wkbLinearRing pour construire des polygones à partir de points ou de lignes (plutôt que wkbPoint ou wkbLine, respectivement). Le bloc de code suivant présente un exemple pour la construction d’un polygone shapefile qui délimite le laboratoire hydraulique de l’Université de Stuttgart. La différence entre l’exemple ci-dessus pour créer un shapefile de ligne est:

  • La projection est EPSG:4326.

  • Les coordonnées de points sont générées à l’intérieur d’un objet ogr.wkbLinearRing étape par étape plutôt que dans une boucle sur les entrées dictionary.

  • Noms de fichiers, de variables et de champs.

shp_dir = r"" + os.path.abspath("") + "/geodata/shapefiles/va4wasserbau.shp"
va_geo = create_shp(shp_dir, layer_name="basemap", layer_type="polygon")

# create .prj file for the shapefile for GIS map applications
with open(shp_dir.split(".shp")[0] + ".prj", "w+") as prj:
    prj.write(get_wkt(4326))

# get basemap layer
lyr = va_geo.GetLayer()

# create polygon points
pts = ogr.Geometry(ogr.wkbLinearRing)
pts.AddPoint(9.103686, 48.744251)
pts.AddPoint(9.104689, 48.744198)
pts.AddPoint(9.104667, 48.743960)
pts.AddPoint(9.103557, 48.744009)

# create polygon geometry
poly = ogr.Geometry(ogr.wkbPolygon)
# build polygon geometry from points
poly.AddGeometry(pts)

# add field to classify building type
field = ogr.FieldDefn("building", ogr.OFTString)
lyr.CreateField(field)

poly_feature_defn = lyr.GetLayerDefn()
poly_feature = ogr.Feature(poly_feature_defn)
poly_feature.SetGeometry(poly)
poly_feature.SetField("building", "Versuchsanstalt")

lyr.CreateFeature(poly_feature)

lyr = None
va_geo = None

Construire Shapefile à partir de JSON

Le chargement des données géométriques à partir de variables définies en ligne est lourd dans la pratique, où les données géospatiales sont souvent fournies sur des plates-formes publiques (p. ex. utilisation des sols ou couverture). L’exemple suivant utilise un fichier JSON généré avec des données de services cartographiques provenant de Institut d’État de l’environnement, de l’étude et de la conservation de la nature de Baden-Wurtemberg (LUBW), où les nœuds polygonaux sont stockés au format de géométrie polygonale WKT ("MultiPolygon (((node1_x node1_y, nodej_x, nodej_y, ... ...)))"):

  • Le fichier JSON (téléchargez hq100-dreisam.json et sauvegardez-le dans un sous-répertoire appelé /geodata/json/) est lu avec Pandas et le shapefile est créé, comme avant, avec la fonction create_shp.

  • La projection est EPSG:25832.

  • Deux champs sont ajoutés sous forme de

    • "tbg_name" (le nom de chaîne original des polygones dans les données LUBW), et

    • "area" (un champ de nombre réel, dans lequel la zone polygonale est calculée en m2 en utilisant polygon.GetArea()).

  • Les géométries polygonales sont dérivées des définitions formatées par WKT dans le champ "wkt_geom" de l’objet pandas dataframe dreisam_inundation.

import pandas as pd

# get data from json file
dreisam_inundation = pd.read_json(r"" + os.path.abspath("") + "/geodata/json/hq100-dreisam.json")

# create shapefile
shp_dir = r"" + os.path.abspath('') + "/geodata/shapefiles/dreisam_hq100.shp"
dreisam_hq100 = create_shp(shp_dir, layer_name="basemap", layer_type="polygon")

# create .prj file for the shapefile for GIS map applications
with open(shp_dir.split(".shp")[0] + ".prj", "w+") as prj:
    prj.write(get_wkt(25832))

# get basemap layer
lyr = dreisam_hq100.GetLayer()

# add string field "tbg_name"
lyr.CreateField(ogr.FieldDefn("tbg_name", ogr.OFTString))

# add string field "area"
lyr.CreateField(ogr.FieldDefn("area", ogr.OFTReal))

for wkt, tbg in zip(dreisam_inundation["wkt_geom"], dreisam_inundation["TBG_NAME"]):
    # create Feature as child of the layer
    feature = ogr.Feature(lyr.GetLayerDefn())
    # assign tbg_name
    feature.SetField("tbg_name", tbg)
    # use WKT format to add a polygon geometry to the Feature
    polygon = ogr.CreateGeometryFromWkt(wkt)
    # define default value of 0 to the area field
    feature.SetField("area", polygon.GetArea())

    feature.SetGeometry(polygon)
    # append the new feature to the basemap layer
    lyr.CreateFeature(feature)

lyr = None
dreisam_hq100 = None

Aussi Les données de GeoJSON peuvent être utilisées pour créer un ogr.Geometry avec ogr.createFromGeoJson(FILENAME):

from osgeo import ogr
geojson_data = """{"type":"Point","coordinates":[1013452.282805,6231540.674235]}"""
point = ogr.CreateGeometryFromJson(geojson_data)
print("X=%d, Y=%d (EPSG:3857)" % (point.GetX(), point.GetY()))
X=1013452, Y=6231540 (EPSG:3857)

Calculer les attributs géométriques

Le bloc de code ci-dessus illustre l’utilisation de polygon.GetArea() pour calculer la zone polygonale à m2^2. La bibliothèque ogr fournit beaucoup d’autres fonctions pour calculer les attributs géométriques des fonctionnalités et voici un résumé:

  • Unifier plusieurs polygones
    wkt_... = ...
    polygon_a = ogr.CreateGeometryFromWkt(wkt_1)
    polygon_b = ogr.CreateGeometryFromWkt(wkt_2)
    polygon_union = polygon_a.Union(polygon_b)

  • Intersecter deux polygones
    polygon_intersection = polygon_a.Intersection(polygon_b)

  • Enveloppe d’un polygone
    env = polygon.GetEnvelope()
    print("minX: %d, minY: %d, maxX: %d, maxY: %d" % (env[0],env[2],env[1],env[3])

  • Coque convex (surface enveloppante) de plusieurs géométries (points, lignes, polygones)
    all_polygons = ogr.Geometry(ogr.wkbGeometryCollection)
    for feature in POLYGON-SOURCE-LAYER: all_polygons.AddGeometry(feature.GetGeometryRef())
    convexhull = all_polygons.ConvexHull()
    Save convexhull to shapefile (use the create_shp function as shown in the above examples or read more at pcjericks’ Github pages)
    Tip: To create a tight hull (e.g., of a point cloud), look for concavehull functions.

  • Longueur (d’une ligne)
    wkt = "LINESTRING (415128.5 5320979.5, 415128.6 5320974.5, 415129.75 5320974.7)"
    line = ogr.CreateGeometryFromWkt(wkt)
    print("Length = %d" % line.Length())

  • Zone (d’un polygone): polygon.GetArea() (voir l’exemple ci-dessus)

  • Exemple pour calculer coordonnées centroïdes de polygones.

Exportation vers d’autres formats

Les exemples ci-dessus traitent uniquement des fichiers .shp, mais d’autres formats peuvent être utiles (par exemple, pour créer des applications Web ou exporter vers Google Earth). À cette fin, les paragraphes suivants illustrent la création de fichiers GeoJSON et KML. Plusieurs autres conversions peuvent être effectuées, non seulement entre les formats de fichiers, mais aussi entre les types de fonctionnalités. Par exemple, les polygones peuvent être créés à partir de nuages point (entre autres avec la méthode ConvexHull mentionnée ci-dessus). Les étudiants intéressés peuvent en savoir plus sur les conversions dans Michael Diener’s Python Geospatial Analysis Cookbook.

GéoJSON

Les fichiers GeoJSON peuvent être facilement créés comme avant, même sans activer un pilote :

triangle = ogr.Geometry(ogr.wkbLinearRing)
triangle.AddPoint(-11717151.498691, 2356192.894805)
triangle.AddPoint(-11717120.446149, 2355586.175893)
triangle.AddPoint(-11719392.059083, 2354012.050842)

polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(triangle)

with open(r"" + os.path.abspath('') + "/geodata/geojson/pitillal-triangle.geojson", "w+") as gjson:
    gjson.write(polygon.ExportToJson())

Pour une gestion de fichier plus robuste et la définition d’une projection, activez le pilote ogr.GetDriverByName("GeoJSON"). Ainsi, la création et la manipulation de fichiers GeoJSON fonctionnent de la même manière que les gestionnaires de fichier de forme montrés ci-dessus.

gjson_driver = ogr.GetDriverByName("GeoJSON")

# make spatial reference
sr = osr.SpatialReference()
sr.ImportFromEPSG(3857)

# create GeoJSON file
gjson = gjson_driver.CreateDataSource("pitillal-full.geojson")
gjson_lyr = gjson.CreateLayer("pitillal-full.geojson", geom_type=ogr.wkbPolygon, srs=sr)

# get layer feature definitions
feature_def = gjson_lyr.GetLayerDefn()
# create new feature
new_feature = ogr.Feature(feature_def)
# assign the triangle from the above code block
new_feature.SetGeometry(polygon)
# add new feature to Layer
gjson_lyr.CreateFeature(new_feature)

# close links to data sources
gjson = None
gjson_lyr = None

KML (Google Earth)

Pour afficher les fonctions point, ligne ou polygone dans Google Earth, les fonctionnalités peuvent être branchées dans le KML (Keyhole Markup Language), similaire à la création d’un fichier GeoJSON, et avec la fonction geometry.ExportToKML:

triangle = ogr.Geometry(ogr.wkbLinearRing)
triangle.AddPoint(-11717151.498691, 2356192.894805)
triangle.AddPoint(-11717120.446149, 2355586.175893)
triangle.AddPoint(-11719392.059083, 2354012.050842)

polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(triangle)

#geojson = poly.ExportToJson()
with open(r"" + os.path.abspath("") + "/geodata/pitillal-triangle.kml", "w+") as gjson:
    gjson.write(polygon.ExportToKML())

De plus, comme pour les fichiers GeoJSON et les shapefiles, les fichiers KML peuvent être générés de manière plus robuste (avec une projection définie par exemple). Tout ce que vous devez faire est d’initier le pilote KML (kml_driver = ogr.GetDriverByName("KML")) et de définir une source de données KML (kml_file = kml_driver.CreateDataSource(FILENAME.KML)).

Géoanalyses
Bibliothèques géospatiales libres de Python
Géoanalyses
Manipulation des grilles