Suite

Comment compter le point shp à l'intérieur du polygone shp en utilisant Python?

Comment compter le point shp à l'intérieur du polygone shp en utilisant Python?


C'est la première fois que je viens ici, alors je vais essayer d'être le plus précis possible.

Ce que j'essaie de faire, c'est de découvrir quel polygone contient le plus de points et d'écrire le résultat dans un fichier texte.

Voici mon script : (tab25Select sont les polygones et chefLyr sont les points)

chefLyr = arcpy.MakeFeatureLayer_management(cheflieuIN, "chefLyr") avec arcpy.da.SearchCursor(tab25Select, ["Numero", "Name"]) en lignes : pour les lignes en lignes : chefTab25Max = 0 chefTab25 = arcpy.GetCount_management(arcpy. SelectLayerByLocation_management(chefLyr, "WITHIN", tab25Select)) arcpy.AddMessage(cheflTab25) pour chf dans cheflTab25 : if chf > chefTab25Max : chefTab25Max = cheflTab25 a = row[0] b = row[1] ficTxt.write("
 n Le polygone contenant le plus de points est :{0}".format(chefTab25Max)) ficTxt.close()

Si quelqu'un peut m'aider :)


Je peux réduire cela à 3 lignes de code, aucun curseur requis !

import arcpy arcpy.SpatialJoin_analysis("Site", "points","in_memory/points_SpatialJoin", "JOIN_ONE_TO_MANY", "KEEP_ALL", "", "INTERSECT") arcpy.Statistics_analysis("points_SpatialJoin", "in_memory/stats", " Join_Count SUM","Id")

Ensuite, triez simplement le tableau pour trouver le polygone avec le plus de points.


L'approche suivante utilise un curseur de recherche et un dictionnaire Python pour effectuer le workflow suivant :

  1. Sélectionnez des points dans chaque entité surfacique
  2. Mettre à jour le dictionnaire avec la clé (OID) et la valeur (nombre de points) pour chaque itération
  3. Trouvez la valeur maximale du point et l'OID correspondant et écrivez dans un fichier texte

import arcpy, os points = r'C:	empmytestpoints.shp' polys = r'C:	empmytestpolys.shp' arcpy.MakeFeatureLayer_management(points, "pointsLyr") data = {} avec arcpy .da.SearchCursor(polys, ["[email protected]", "[email protected]"]) comme curseur : pour la ligne dans le curseur : arcpy.SelectLayerByLocation_management("pointsLyr", select_features = row[1]) count = int(arcpy.GetCount_management( "pointsLyr").getOutput(0)) data[row[0]] = count # Écrire la sortie dans un fichier texte text_file = open(r"C:	empOutput.txt", "w") text_file.write ("OID = %s; count = %s" % (max(data, key = data.get), data[max(data, key = data.get)])) text_file.close()


Comment obtenir la liste des points à l'intérieur d'un polygone en python ?

J'ai beaucoup cherché et je n'ai pas trouvé de réponse pratique à ma question. J'ai un polygone. Par exemple:

Je veux obtenir une liste de tous les points à l'intérieur de ce polygone de frontière. J'ai beaucoup entendu parler des techniques de triangulation polygonale ou linéaire/inondation/intersection/. algorithmes de remplissage. mais je ne peux pas vraiment trouver un moyen efficace de mettre cela en œuvre. Ce poly est petit, imaginez un polygone avec 1 milliard de points. J'utilise maintenant le polygone de dessin PIL pour remplir le poly de couleur rouge et boucler à l'intérieur pour trouver des points rouges. C'est une technique horriblement lente :

Je cherche un moyen Pythonic de résoudre ce problème. Merci à tous.


Qu'est-ce qu'un point dans un polygone ?

Le point dans le polygone est assez simple.

Voici une carte des quartiers de New York, tels que le West Village, Astoria, Washington Heights, Tribeca, etc.

Dans quel quartier se trouve une adresse ?

C'est le type de question à laquelle nous pouvons répondre en utilisant le point dans la recherche de polygones.

Chaque quartier est représenté sur la carte comme un polygone, qui, si vous vous souvenez des cours de mathématiques, est un Figure en 2 dimensions avec des côtés droits, comme un triangle, un rectangle, un pentagone, un hexagone, etc.

UNE indiquer est une paire de coordonnées de latitude et de longitude - le point rouge sur la carte est Times Square avec les coordonnées -73.9855 et 40.7580.

Vous pouvez voir que le point rouge se situe dans le polygone entouré en vert, qui est le 'Midtown-Midtown Sud' quartier dans cet ensemble de données.

J'ai choisi la solution de facilité et j'ai recherché les coordonnées de Times Square sur Google, mais vous pouvez utiliser n'importe quelle adresse et il vous suffit de la géocoder pour obtenir les coordonnées de latitude et de longitude.


Obtenir des informations sur un fichier de formes

Vous pouvez télécharger un fichier de formes à partir du NHGIS. Pour cet exemple, nous utiliserons le shapefile des comtés américains en 1850.

Lorsque vous décompressez le répertoire qui contient le shapefile, vous remarquerez qu'il contient en fait plusieurs fichiers. Le fichier qui se termine par .shp est le fichier de formes lui-même, mais les autres fichiers sont également importants. Le fichier qui se termine par .prj contient les informations de projection pour le fichier de formes et le fichier qui se termine par .dbf contient les attributs associés aux polygones, aux points ou aux lignes du fichier de formes. (R vous permettra de lire ce fichier DBF seul si vous voulez juste ces données.)

Nous pouvons voir quel type d'informations est disponible dans le fichier de formes en utilisant la commande ogrinfo avec quelques indicateurs qui indiquent à la commande de nous donner toutes les informations du fichier. En supposant que vous soyez dans votre terminal dans le même répertoire que le fichier de formes, vous pouvez exécuter la commande suivante.

Cette commande retournera beaucoup d'informations :

Cette sortie nous indique que le fichier contient des polygones au lieu de points ou de lignes ( Geometry: Polygon ) et qu'il y a 1 632 polygones dans le fichier ( Feature Count: 1632 ). Nous pouvons également voir que le fichier est projeté en utilisant une projection à surface égale d'Albers. C'est un bon choix pour une projection, mais cela rendra le travail avec le fichier de formes dans R plus difficile. La sortie nous indique également les variables associées à chaque polygone, telles que STATE et COUNTY . Dans ce fichier, les informations sont principalement l'identification des polygones, mais d'autres fichiers de formes contiendront des données réelles. Notez que nous pourrions obtenir les mêmes informations dans R en utilisant le wrapper rgdal autour de la commande ogrinfo. La fonction ogrInfo() prend deux arguments : le chemin d'accès au répertoire contenant le fichier de formes et le nom du fichier de formes sans l'extension .shp.


ÉDITER: Après avoir lu les commentaires de @MrE et @MartinR, je propose maintenant cette méthode d'échantillonnage de rejet. Bien que cela puisse manquer fréquemment dans un polygone avec une grande boîte englobante par rapport à sa surface, par exemple. une étoile de Noël à 8 branches avec un petit cercle intérieur.

C'était mon idée originale, mais ça n'avait pas l'air si bien le matin.

Tout d'abord, découvrez comment sélectionner un point dans un polygone.

Ensuite, appelez simplement cela dans une compréhension de liste pour générer autant de points que vous le souhaitez.

Mon approche consiste à sélectionner x au hasard dans le polygone, puis à contraindre y .

Cette approche ne nécessite pas NumPy et produit toujours un point dans le polygone.

L'échantillonnage de rejet a été proposé dans les commentaires sur l'autre réponse. Le problème avec l'échantillonnage de rejet est que l'aire d'un polygone peut être une fraction arbitrairement petite de sa boîte englobante, par exemple :

L'échantillonnage de rejet prendra en moyenne $1sur ε$ tentatives pour générer un point d'échantillonnage à l'intérieur de ε_poly(ε) , et $1sur ε$ peut être arbitrairement grand.

Une approche plus robuste est la suivante :

Triangule le polygone et calcule l'aire de chaque triangle.

Choisissez le triangle $t$ contenant l'échantillon, en utilisant une sélection aléatoire pondérée par l'aire de chaque triangle.

Choisissez un point aléatoire $x, y$ uniformément dans le carré unitaire.

Si $x + y > 1$ , utilisez le point $1-x, 1-y$ à la place. L'effet de ceci est de s'assurer que le point est choisi uniformément dans le triangle rectangle unité avec les sommets $(0, 0), (0, 1), (1, 0)$

Appliquez la transformation affine appropriée pour transformer le triangle rectangle unitaire en triangle $t$ .

C'est bien si $k$ est petit, comme indiqué par l'OP. (Si $k$ est grand, alors vous voudriez vectoriser la construction des points dans NumPy.)

Des inquiétudes ont été exprimées dans les commentaires selon lesquelles la transformation affine entraînerait le rapprochement des points aléatoires le long d'un axe plutôt que de l'autre, en raison de la mise à l'échelle différente sur les deux axes. Mais ce n'est pas une bonne façon d'y penser, car les points les plus proches ne sont pas préservés par les transformations affines. (Un point peut avoir un voisin le plus proche différent après la transformation.) Au lieu de cela, n'oubliez pas que les transformations affines sont linéaires sur les deux axes, de sorte que la probabilité qu'un point aléatoire apparaisse dans une région reste proportionnelle à la superficie de la région après la transformation. Cela signifie que si les points étaient uniformément répartis dans le triangle rectangle unitaire, ils restent uniformément répartis dans le triangle transformé. Voici une démo rapide :

Si le chiffre ne vous convainc pas, voici un code qui imprime les décalages absolus moyens le long des deux axes entre les points aléatoires et leurs voisins les plus proches :

La sortie typique montre que les moyennes sur les deux axes sont très proches même si les coordonnées $x$ ont été mises à l'échelle et les coordonnées $y$ n'ont pas :


Lire et écrire Geojson

Travailler avec des données geojson. Par défaut, la source de données est créée en tant que WGS84.

  • Préparation des données
  • Lire toutes les données
  • Lecture et itération sur chaque fonctionnalité
  • Créer un nouveau fichier geojson
  • Lecture à partir du fichier geojson

Lire et écrire WKT

Travailler avec des données wkt. La collection de géométries et les géométries uniques sont prises en charge. Par défaut, la source de données est créée en tant que WGS84.

L'objet wkt a quelques paramètres :

as_geometry_collection: renvoie une collection de géométrie identique lorsque les données sont une seule géométrie par méthode collection.

srid: SRID initial pour WKT.

  • Lire toutes les données
  • Lecture et itération sur chaque géométrie
  • Création d'un nouveau fichier wkt
  • Lecture à partir du fichier wkt

Lire et écrire Shapefile

Travailler avec un fichier de formes en lecture et en écriture. Est pris en charge les fichiers de formes compressés au format .zip et .rar. Par défaut, la source de données est créée en fonction de la projection présente sur le fichier .prj. obs : la lecture et l'écriture des fichiers .rar sont disponibles uniquement pour le système d'exploitation Linux. Seule la classe vectorio.compress.Rar (moteur à compresser) a cette restriction. Les autres classes sont disponibles pour n'importe quel système d'exploitation.

  • Préparation des données
  • Lire toutes les données du fichier .shp
  • Lecture et itération sur chaque fonctionnalité du fichier .shp.
  • Création d'un nouveau shapefile (Sont créés les fichiers .shp, .shx, .dbf, .prj)
Lire et écrire Shapefile compressé

Par défaut, l'algorithme recherchera de manière récursive les fichiers .shp, .shx, .dbf, .prj à l'intérieur du fichier compressé. L'algorithme recherchera le premier fichier de chaque extension, si le fichier compressé contient 2 (ou plus) fichiers .shp, ou 2 (ou plus) fichiers .prj, sera obtenu le premier fichier .shp et le premier fichier .prj .

Lire et écrire KML

Reprojeter un vecteur

La reprojection spatiale fonctionne avec le même type de géographie qui implémente l'interface IVectorIO. Si le srid d'entrée (in_srid) est omis, le srid de la géométrie sera utilisé.

Par défaut, le wkt est en référence spatiale WGS84.

Par défaut, le geojson est en référence spatiale WGS84.

Basculement rapide entre les données géographiques

Pour l'exécution de l'interrupteur rapide doit être utilisé le VecteurComposite présent sur l'emballage vectorio.vector.

Passage rapide de geojson à wkt
  • Préparation des données
  • Lecture de toute la géométrie de geojson en tant que wkt
  • Itération sur toutes les géométries comme wkt
  • Création d'un fichier wkt
Passage rapide de wkt à shapefile en tant que zip
  • Lecture de toute la géométrie de wkt
  • Itérer sur toutes les géométries
  • Création d'un shapefile sous forme de zip
Rechercher une zone UTM à partir de la géométrie

La méthode zone_from_biggest_geom obtenir une seule zone qui a la plus grande géométrie. Par exemple, si une grande géométrie se trouve dans la zone UTM 25N ou 26N, cette méthode calcule la surface (pour le polygone) ou la longueur (pour la ligne) de la géométrie et renvoie la zone UTM où la surface ou la longueur est la plus grande.

Signature: zone_from_biggest_geom(self, inp_ds : DataSource, wkt_prj_for_metrics : str, in_wkt_prj=PRJ_WGS84)

  • wkt_prj_for_metrics: obligatoire
  • inp_ds: obligatoire
  • in_wkt_prj: Optionnel. Cas utilisé, la source de données d'entrée n'a pas de CRS.
  • Les points sont ignorés dans ce calcul. Cette méthode n'est pas recommandée pour les géométries composées uniquement de points (MultiPoint, GeometryCollection of points e etc.), car la zone UTM renvoyée est basée sur des métriques de surface/longueur. Bientôt, sera implémentée une méthode pour obtenir la zone UTM uniquement pour les géométries ponctuelles.
Zone UTM pour les géométries avec des erreurs de topologie

Si votre géométrie comporte des erreurs de topologie, vous devez utiliser la méthode UTMZone().zones_from_iteration(wkt) passer un WKT comme argument. Cette méthode renverra un ensemble de zones utm que le wkt croise.

Lever une exception pour les avertissements de Gdal

Pour utiliser l'exception des avertissements gdal doit utiliser le décorateur gdal_warning_as_exception présente sur vectorio.gdal paquet. Ce décorateur lancera l'erreur lorsque le Est valable() méthode de géométrie() méthode sera utilisée.

Exceptions possibles
  • GDALSelfIntersectionGéométrie: exception levée lorsqu'un polygone contient une auto-intersection.
  • GDALMauvaisPolygoneFermé: exception levée lorsqu'un polygone ne se ferme pas correctement.
  • GDALExceptionInconnue: exception levée lorsqu'une erreur inconnue se produit.

Obs : Toutes les exceptions sont disponibles sur forfait vectorio.exceptions

Reprojeter une source de données directement

Pour effectuer une reprojection spatiale, utilisez la classe VectorIOReprojected mentionnée ci-dessus. Cependant, au même moment il faudra reprojeter directement une source de données, pour cela utiliser la classe DataSourceReprojected.

Signature: DataSourceReprojected(inp_ds : DataSource, in_srid : int=None, out_srid : int=None, in_wkt_prj=None, out_wkt_prj=None, use_wkt_prj : bool=False)

  • inp_ds: obligatoire. Source de données qui sera reprojetée.
  • in_srid: optionnel. Cas utilisé, la source de données d'entrée n'a pas de CRS. (Utilisé uniquement lorsque le drapeau use_wkt_prj est Faux).
  • out_srid: Obligatoire si le use_wkt_prj est Faux.
  • in_wkt_prj: optionnel. Projection WKT sur modèle OGC. Cas utilisé, la source de données d'entrée n'a pas de CRS. (Utilisé uniquement lorsque le drapeau use_wkt_prj est Vrai).
  • out_wkt_prj: Obligatoire si le use_wkt_prj est Vrai. Projection WKT sur modèle OGC.
  • use_wkt_prj: Booléen pour définir si la source de données sera reprojetée avec le SRID ou avec la Projection WKT.

Compter les sommets

Toutes les classes vectorielles ont la méthode pour compter tous les sommets. Ci-dessous, sont illustrés par un WKT, cependant, cette méthode fonctionne également avec Shapefile, GeoJson.

Signature: vertices_by_feature(ds: DataSource) -> dict


Vous rencontrerez généralement des données vectorielles sous forme de fichiers de formes (extension .shp ). Cependant, vous rencontrerez également des vecteurs sous forme d'ensembles de données de points qui incluent les coordonnées des observations et des informations sur ces points (souvent .csv ). À la base, tous les fichiers vectoriels sont simplement des collections de points. Pour les fichiers basés sur des points, cette déclaration est évidente. Les lignes et les polygones sont des collections de points, donc encore une fois, nous voyons que les classes d'objets qui composent les vecteurs sont « simplement » des collections de points.

Chargeons un fichier de formes. Concrètement, nous allons charger les polygones qui composent les battements policiers de la ville de Chicago. 3 Pour lire dans le shapefile qui comprend ces données spatiales, nous utiliserons la fonction readOGR() du package rgdal . La fonction a besoin de deux arguments : dsn , qui est généralement le dossier qui contient le fichier de formes, et layer , qui est généralement le nom du fichier de formes (en omettant l'extension .shp). Dans notre cas, le nom du dossier est « ChicagoPoliceBeats », et le nom du fichier de formes est « chiBeats.shp ». Chargeons le shapefile !

Noter: rgdal et sp font du bon travail pour le travail spatial dans R, mais si vous le trouvez un peu lent - ou si vous voulez être un peu plus à la pointe - je vous suggère de consulter le package sf. C'est plus récent, donc c'est

Maintenant, vérifions quelques éléments : la classe de l'objet que nous venons de lire dans R, ses dimensions et son tracé.

Cool, non ? Alors, qu'est-ce que cela signifie que ce « bloc de données de polygones spatiaux » a les dimensions 277 et 4 ? Le premier nombre est généralement le nombre de formes distinctes (caractéristiques) dans le fichier de formes. Dans notre cas, nous avons 277 polygones. Et le deuxième numéro ? Le deuxième nombre nous donne le nombre de champs, essentiellement le nombre de variables (alias caractéristiques) pour lesquelles nous avons des informations pour chaque polygone.


Tiers de confiance

Il appartient entièrement à l'attaquant de créer un système où il peut prouver que les données sont sécurisées. Une fois qu'ils ont les données en main, il est trop tard. C'est pourquoi ils doivent transmettre les données directement à un tiers en qui eux-mêmes et l'entreprise ont confiance. Il est également crucial que le tiers fournisse les informations nécessaires pour prouver à la fois qu'il dispose des données (par exemple, les noms et les tailles des fichiers) et que personne n'y a accédé (par exemple, un nombre de téléchargements).

Voici comment une attaque pourrait se dérouler :

  1. L'attaquant convainc un employé de télécharger les données sur un compte tiers via l'ingénierie sociale. L'attaquant contrôle ce compte et peut potentiellement télécharger les données.
  2. L'attaquant montre la preuve du téléchargement à l'entreprise et fait la demande.
  3. Lorsque la rançon est payée, l'attaquant remet le compte tiers à l'entreprise, qui verrouille ensuite l'attaquant hors du compte.
  4. L'entreprise vérifie auprès du tiers que les données n'ont pas été consultées, puis supprime le compte et toutes les données.

Notez à quel point les vecteurs d'attaque sont extrêmement limités. L'attaque ne fonctionne que si une personne en qui l'entreprise a confiance télécharge les données. De plus, l'attaquant est très limité dans l'infiltration de l'entreprise : s'il se donne par inadvertance un accès inauditable aux données, cela rendrait le pouf ci-dessus dénué de sens.


Comment traiter les données météorologiques GRIB2 pour les applications aéronautiques (Shapefile)

Écrit par
Elliott Wobler
Ingénieur de solutions techniques

Ce didacticiel explique comment utiliser les données de prévision globales de Spire Weather au format GRIB2 à l'aide de Python.

Ce didacticiel utilise l'entrée Shapefile.

Aperçu

Ce didacticiel explique comment utiliser les données de prévision mondiales de Spire Weather au format GRIB2 à l'aide de Python.

Si vous n'avez jamais travaillé avec des données GRIB2 auparavant, il est recommandé de commencer par le didacticiel de démarrage, car celui-ci abordera des sujets légèrement plus avancés.

Plus précisément, ce didacticiel montre comment récupérer des turbulences en air clair à différents niveaux de vol dans les limites d'un polygone complexe (par exemple, une frontière nationale).

À la fin de ce tutoriel, vous saurez comment :

  1. Charger des fichiers contenant des messages GRIB2 dans un programme Python
  2. Inspectez les données GRIB2 pour comprendre quelles variables météorologiques sont incluses
  3. Filtrez les données GRIB2 sur les variables météorologiques d'intérêt à différents niveaux de vol verticaux
  4. Recadrez les données GRIB2 dans la zone définie par une entrée Esri Shapefile
  5. Convertissez les données GRIB2 transformées en un fichier de sortie CSV pour une analyse et une visualisation plus poussées

Téléchargement des données

La première étape consiste à télécharger les données de l'API de fichier de Spire Weather.

Ce didacticiel s'attend à ce que les messages GRIB2 contiennent des données NWP du paquet de données Spire's Aviation.

Pour plus d'informations sur l'utilisation des points de terminaison de l'API de fichier Spire Weather, veuillez consulter la documentation de l'API et la FAQ.

La FAQ contient également des exemples détaillés expliquant comment télécharger plusieurs fichiers de prévision à la fois à l'aide de cURL.

Pour les besoins de ce didacticiel, il est supposé que les données GRIB2 ont déjà été téléchargées avec succès, sinon vous pouvez obtenir un échantillon de prévisions météorologiques mondiales ici.

Comprendre les noms de fichiers

Les fichiers téléchargés à partir des produits API de Spire Weather partagent tous la même convention de dénomination de fichier.

En regardant simplement le nom du fichier, nous pouvons déterminer :

  • la date et l'heure d'émission de la prévision
  • la date et l'heure de validité des données prévisionnelles
  • la résolution horizontale des données de prévision globales
  • les variables de données météorologiques incluses dans le fichier (voir la liste complète des ensembles de données commerciales Spire Weather)

Pour plus d'informations sur ce qui précède, veuillez vous référer à notre FAQ.

Configuration requise pour Python

Les packages Python suivants sont requis pour ce didacticiel.

Bien que l'utilisation de conda ne soit pas strictement requise, c'est la méthode officiellement recommandée pour installer PyNIO (voir le lien ci-dessous).

Une fois qu'un environnement conda a été créé et activé, les commandes suivantes peuvent être exécutées directement :

Inspection des données

Après avoir téléchargé les données et configuré un environnement Python avec les packages requis, l'étape suivante consiste à inspecter le contenu des données afin de déterminer les variables météorologiques auxquelles il est possible d'accéder. L'inspection des données peut être effectuée uniquement avec PyNIO, mais dans ce didacticiel, nous utiliserons plutôt PyNIO comme moteur pour xarray et chargerons les données dans un ensemble de données xarray pour la transformation. Notez qu'une importation explicite de PyNIO n'est pas requise, tant qu'il est correctement installé dans l'environnement Python actif.

Tout d'abord, importez le package xarray :

Ensuite, ouvrez les données GRIB2 avec xarray en utilisant PyNIO comme moteur (notez que les données GRIB2 doivent provenir de Spire’s Aviation paquet de données):

ds = xr.open_dataset("path_to_aviation_file.grib2", engine="pynio")

Enfin, pour chacune des variables, imprimez la clé de recherche, le nom lisible par l'homme et les unités de mesure :

La sortie de ce qui précède devrait ressembler à ceci, donnant un aperçu clair des champs de données disponibles :

Étant donné que ce didacticiel couvre l'utilisation des données GRIB2 à différents niveaux verticaux, nous devons également inspecter le champ level_type :

La sortie de ce qui précède devrait ressembler à ceci :

Notez que plusieurs des variables incluses ont la même valeur level_type d'altitude spécifique au-dessus du niveau moyen de la mer (m) .

Cela signifie que la même méthodologie exacte peut être utilisée pour traiter ces variables à différents niveaux de vol, bien que pour les besoins de ce didacticiel, nous nous en tiendrons uniquement à la turbulence en air clair : CAT_P0_L102_GLL0 .

Recadrez les données GRIB2 dans la zone définie par une entrée Esri Shapefile
.

Pour mieux comprendre les valeurs de niveau de vol disponibles, nous pouvons commencer par imprimer des métadonnées pour la variable de turbulence en air clair :

La sortie de ce qui précède devrait ressembler à ceci :

À la fin de la première ligne, la liste contient non seulement les noms d'index lat_0 et lon_0 auxquels nous nous attendrions après le didacticiel de base, mais aussi un nouveau appelé lv_AMSL0 .
Nous savons d'après notre déclaration d'impression précédente que lv_AMSL0 fait référence à l'altitude au-dessus du niveau moyen de la mer, et la sixième ligne de la dernière sortie d'impression fournit un aperçu des valeurs lv_AMSL0. A partir de ces informations, nous pouvons deviner qu'il existe un troisième index pour ces données appelé lv_AMSL0 qui correspond au niveau vertical des données.
Nous pouvons imprimer des métadonnées sur un index avec la même méthode que celle utilisée pour inspecter les variables :

La sortie de ce qui précède devrait ressembler à ceci :

La première ligne de la sortie nous indique qu'il y a 36 valeurs possibles, exprimées sous la forme (lv_AMSL0: 36) , et les lignes 2 à 5 répertorient ces valeurs exactes (en mètres).

Cela signifie que pour chaque paire unique de coordonnées lat_0 et lon_0, il existe 36 valeurs de turbulence en air clair distinctes à différentes altitudes qui peuvent être récupérées en spécifiant l'une des valeurs lv_AMSL0.

Afin d'améliorer la convivialité de notre script, c'est une bonne idée de créer un dictionnaire qui mappe les niveaux de vol à leurs valeurs associées en mètres :

Maintenant, nous pouvons implémenter notre script pour attendre une valeur d'entrée de niveau de vol de FL350 plutôt que d'exiger que l'utilisateur du script connaisse la valeur exacte en mètres.

Traitement des données

Maintenant que nous savons quelles variables météorologiques et quels niveaux verticaux sont inclus dans les données GRIB2, nous pouvons commencer à traiter notre jeu de données xarray.

Filtrage des données xarray vers une variable spécifique

Avec xarray , le filtrage du contenu de l'ensemble de données sur une seule variable d'intérêt est très simple :

La sélection de plusieurs variables est également possible en utilisant un tableau de chaînes comme entrée : ds = ds.get([var1, var2])

Il est recommandé d'effectuer cette étape avant de convertir le jeu de données xarray en un DataFrame pandas (plutôt que de filtrer le DataFrame ultérieurement), car cela minimise la taille des données converties et réduit donc le temps d'exécution global.

Conversion des données xarray en pandas.DataFrame

Pour convertir l'ensemble de données xarray filtré en un DataFrame pandas, exécutez simplement la commande suivante :

Filtrage des pandas.DataFrame à un niveau de vol spécifique

En utilisant les informations recueillies auprès du Comprendre les valeurs de niveau de vol section, nous pouvons maintenant filtrer le DataFrame à un niveau de vol spécifique.

La première chose que nous faisons est de copier les valeurs d'index de niveau vertical (en mètres) dans une nouvelle colonne DataFrame :

Ensuite, nous créons l'expression que nous utiliserons pour filtrer le DataFrame.

En utilisant la carte du dictionnaire créée ci-dessus, nous pouvons filtrer les données de turbulence en air clair à un niveau de vol de FL350 avec l'expression suivante :

Enfin, nous appliquons le filtre au DataFrame :

Notre DataFrame a maintenant été filtré pour n'inclure que les données de turbulence en air clair à un niveau de vol de FL350.

Chargement d'un fichier de forme Esri avec le package Python GDAL

Bien que le package que nous avons installé avec conda s'appelle gdal , nous l'importons dans Python sous le nom osgeo . Il s'agit d'une abréviation de l'Open Source Geospatial Foundation, par laquelle GDAL/OGR (un logiciel libre et open source) est concédé sous licence.

Le format Shapefile, développé à l'origine par Esri, est une norme courante dans le monde des systèmes d'information géographique. Il définit la géométrie et les attributs des entités référencées géographiquement dans trois fichiers ou plus avec des extensions de fichier spécifiques. Les trois extensions de fichier obligatoires sont .shp , .shx et .dbf — et tous ces fichiers de composants associés doivent être stockés dans le même répertoire de fichiers. De nombreux fichiers de formes préexistants peuvent être téléchargés gratuitement en ligne (par exemple, les frontières nationales, les zones économiques exclusives, etc.) et des formes personnalisées peuvent également être créées dans une variété d'outils logiciels gratuits. Dans cet exemple, nous utilisons le pays de l'Italie comme polygone complexe, mais cela pourrait tout aussi bien être la zone entourant un aéroport ou une autre petite région.

Lors de l'ouverture d'un Shapefile avec GDAL, il suffit de pointer sur le fichier avec l'extension .shp. Cependant, il est nécessaire que les autres fichiers de composants existent dans le même répertoire. Si nous ouvrons un fichier appelé italy.shp , il devrait au moins y avoir des fichiers nommés italy.shx et italy.dbf dans le même répertoire.

Une fois que nous avons défini le pilote et le chemin vers notre polygone Shapefile, nous pouvons le charger dans notre script comme ceci :

Obtenir la zone de délimitation qui contient une zone Shapefile

Finalement, nous recadrerons les données dans la zone précise définie par le Shapefile, mais il s'agit d'un processus coûteux en calcul, il est donc préférable de limiter d'abord la taille des données. En théorie, nous pourrions sauter complètement l'étape du recadrage dans une simple boîte, mais en pratique, cela vaut la peine de le faire pour réduire le temps d'exécution global.

GDAL facilite le calcul du cadre de délimitation grossier qui contient une zone de fichier de formes complexe :

Les valeurs de coordonnées sont ensuite accessibles individuellement à partir du tableau résultant :

L'ordre des coordonnées géospatiales est une source courante de confusion, alors prenez soin de noter que la fonction GetExtent de GDAL renvoie un tableau où les valeurs de longitude précèdent les valeurs de latitude.

Recadrage de pandas.DataFrame dans une zone de délimitation géospatiale

Maintenant que les données filtrées sont converties en un DataFrame pandas et que nous avons les limites contenant notre zone d'intérêt, nous pouvons recadrer les données dans une simple boîte.

La première étape de ce processus consiste à décompresser les valeurs de latitude et de longitude de l'index DataFrame’s, accessible via les noms d'index lat_0 et lon_0 :

Bien que les valeurs de latitude soient déjà dans la plage standard de -90 dégresser à +90 degrés, les valeurs de longitude sont dans la plage de 0 à +360.

Pour faciliter l'utilisation des données, nous convertissons les valeurs de longitude dans la plage standard de -180 degrés à +180 degrés:

Avec les données de latitude et de longitude maintenant dans les plages de valeurs souhaitées, nous pouvons les stocker en tant que nouvelles colonnes dans notre DataFrame existant :

Ensuite, nous utilisons les valeurs de la zone de délimitation de la section précédente (les composants du tableau bbox) pour construire les expressions de filtre DataFrame :

Enfin, nous appliquons les filtres à notre DataFrame existant :

Le DataFrame résultant a été rogné jusqu'aux limites de la zone qui contient la zone complexe du Shapefile.

Recadrage du pandas.DataFrame aux limites précises d'une zone Shapefile

Afin de recadrer le DataFrame aux limites précises de la zone complexe du Shapefile, nous devrons vérifier chaque paire de coordonnées dans nos données. Semblable à la section précédente où nous avons remappé chaque valeur de longitude, nous effectuerons cette action avec une expression de carte.

Pour transmettre chaque paire de coordonnées à la fonction de carte, nous créons une nouvelle colonne DataFrame appelée point où chaque valeur est un tuple contenant à la fois la latitude et la longitude :

Nous pouvons ensuite passer chaque valeur de tuple de paire de coordonnées dans la fonction map, ainsi que la zone Shapefile précédemment chargée, et les traiter dans une fonction appelée check_point_in_area que nous définirons ci-dessous. La fonction check_point_in_area retournera True ou False pour indiquer si la paire de coordonnées fournie est à l'intérieur de la zone ou non. En conséquence, nous allons nous retrouver avec une nouvelle colonne DataFrame de valeurs booléennes appelée inArea :

Une fois que la colonne inArea est remplie, nous effectuons un filtre simple pour supprimer les lignes où la valeur inArea est False . Cela supprime efficacement les données de tous les emplacements de points qui ne se trouvent pas dans la zone Shapefile’s :

Bien sûr, le succès de ce qui précède dépend de la logique à l'intérieur de la fonction check_point_in_area, que nous n'avons pas encore implémentée. Étant donné que la zone Shapefile a été chargée avec GDAL, nous pouvons utiliser une méthode de géométrie GDAL appelée Contains pour vérifier rapidement si la zone contient un point spécifique. Pour ce faire, la paire de coordonnées doit d'abord être convertie en une géométrie wkbPoint dans GDAL :

Tous les fichiers de formes sont composés d'une ou plusieurs fonctionnalités, donc pour nous assurer que notre code est robuste, nous devrons vérifier chaque fonctionnalité de composant individuellement. Les entités peuvent être récupérées par leur index numérique avec area.GetFeature(i) , et le nombre total d'entités peut être récupéré avec area.GetFeatureCount() . En utilisant ces deux fonctions avec la méthode Contains mentionnée ci-dessus, nous pouvons parcourir chaque entité et vérifier si elle contient la géométrie du point :

En rassemblant les pièces, voici ce que nous obtenons pour notre dernière fonction check_point_in_area :

Comme vous pouvez le voir, la seule variable renvoyée par la fonction check_point_in_area est une valeur booléenne indiquant si le point spécifié est dans la zone ou non. Ces valeurs booléennes remplissent ensuite la nouvelle colonne inArea DataFrame. Cela nous permet d'appliquer le filtre d'en haut et d'obtenir les données recadrées avec précision que nous voulons :

Analyser l'heure de prévision à partir du nom de fichier

Chaque fichier individuel contient des données de prévisions météorologiques mondiales pour le même moment.

En utilisant nos connaissances de la Comprendre les noms de fichiers section de ce didacticiel, et en supposant que l'argument du nom de fichier est dans le format attendu, nous pouvons écrire une fonction pour analyser l'heure de prévision valide à partir du nom de fichier :

Ensuite, nous pouvons créer une nouvelle colonne DataFrame pour le temps qui stocke cette valeur sous forme de chaîne dans chaque ligne :

Bien qu'il puisse sembler inutile de stocker la même valeur d'horodatage exacte dans chaque ligne, il s'agit d'une étape importante si nous voulons éventuellement concaténer notre DataFrame avec des données de prévision pour d'autres moments (démontré ci-dessous dans Traitement de plusieurs fichiers de données).

Enregistrement des données dans un fichier de sortie CSV

Effectuez un filtre final sur notre DataFrame pour sélectionner uniquement les colonnes que nous voulons dans notre sortie, où la variable est une chaîne comme "CAT_P0_L102_GLL0":

Enregistrez le DataFrame traité dans un fichier CSV de sortie :

La définition du paramètre index=False garantit que les colonnes d'index DataFrame ne sont pas incluses dans la sortie. De cette façon, nous excluons les valeurs lat_0 , lon_0 et lv_AMSL0 car nous avons déjà des colonnes pour la latitude et la longitude remappées (et toutes les données restantes sont à la même altitude après notre filtre de niveau de vol).

Veuillez noter que la conversion de GRIB2 en CSV peut entraîner fichiers de très grande taille, especially if the data is not significantly cropped or filtered.

Processing Multiple Data Files

It is often desirable to process multiple data files at once, in order to combine the results into a single unified CSV output file.

For example, let’s say that we have just used the Spire Weather API to download a full forecast’s worth of GRIB2 data into a local directory called forecast_data/ . We can then read those filenames into a list and sort them alphabetically for good measure:

From here, we can iterate through the filenames and pass each one into a function that performs the steps outlined in the Processing the Data section of this tutorial.

Once all of our final DataFrames are ready, we can use pandas to concatenate them together like so (where final_dataframes is a list of DataFrames):

We end up with a combined DataFrame called output_df which we can save to an output CSV file like we did before:

Complete Code

Below is an operational Python script which uses the techniques described in this tutorial and also includes explanatory in-line comments.

The script takes four arguments:

The weather data variable of interest

The local directory where the Aviation Weather GRIB2 data is stored

The local directory where the Shapefile components are stored

The flight level value ( FL100 to FL450 ) for filtering the data to a specific altitude

For example, the script can be run like this:

Here is the complete code:

Final Notes

Using the CSV data output from our final script, we can now easily visualize the processed data in a free tool such as kepler.gl. We can also set thresholds for alerts, generate statistics, or fuse with other datasets.

Spire Weather also offers pre-created visualizations for some data variables through the Web Map Service (WMS) API which you can read more about here.

For additional code samples, check out Spire Weather’s public GitHub repository.

Would you like to book a consultation?

Learn more about our Weather APIs and how Spire Weather can help you enable the data advantage.


How to count point shp inside polygon shp using Python? - Systèmes d'information géographique

Until now you have made a large number of exercises using GIS. These exercises have given you an insight in the vast possibilities of GIS. Your problem now (after having finished this manual) is to apply these techniques to your situation, with your data. You probably do not have your data in similar files as the ones you used for the exercises in this manual. There are several sources for these types of files most likely your country will have a national agency/institute for this type of information. Probably the Geography department of your national University will be able to point you to the right people. There are many public access sources on the internet, as www.esri.com, General Bathymetric Chart of the Oceans (http://www.ngdc.noaa.gov/mgg/gebco/gebco.html), etc. A look in search engines on the web will give results for specific data.

However, sometimes you will not be able to find the data you need, or you do not have the resources to acquire them. So the only solution is to make the files yourself. If you develop your own data sets in a GIS, you have to realize a very important issue, which is Quality of Data . Specialized GIS institutions can make high quality digitized Themes/layers as they use digitising tables and special software. With limited resources it will be difficult for you to reach their standards. However, with ArcView you can make Themes to analyse data of your own study area. In this chapter some guidelines are provided to create your own Themes/layers.

As you can imagine, there are three principle aspects when you are creating your own Themes/layers:

Georeferencing of the data. This is of utmost importance to avoid spatial distortion in areas and distances of your data as discussed in the Chapter Projection, on page 40.

Digitising your maps/images. The result of this will be a Theme (or several Themes) with points, and/or polygons, and/or lines.

Registering the images. It is important that ArcView knows the exact place in the mapped area where your data must be allocated to. When you register an image, map, data, you give ArcView georeferencing information of the image involved.

Once you have registered a Theme, map, or something similar, field data can be added to the Theme (see chapter Putting field data into a GIS, page 36).

18.1. Registering and digitizing an image

Making your own GIS starts with maps, aerial photographs, etc, of the area you are interested in, after which data is added. An example of a map is presented in Figure 18.1, where the results of a biodiversity study in the floodplains of Pais Pesca are presented. In the figure you see sampling stations (red circle with cross), water bodies (dark blue), villages (brown), rice fields (green), roads (dark brown/red), etc. What is more important, the latitudes and longitudes are clearly indicated. With this information this figure is georeferenced, and can be registered into ArcView.

Registration of images in ArcView can be done with an Avenue script, written by Mark Cederholm [42] . This extension will be used in the next exercise.

FIGURE 18.1
Pais Pesca floodplain biodiversity study

18.1.1. Registering

Only ‘.tiff’ and ‘.bmp’ images can be registered into ArcView with Cederholm’s extension. So when you scan an image, or make the image yourself, you have to make sure you save it as a ‘.tif’ or ‘.bmp’ file. You can find the file containing the above figure in the folder ‘RW_05_Geo’ with the file name: ‘Pais_Pesca_Biodiv_floodplain.tif’. First you have to copy this folder onto your hard disk, before you start working with this file (from the hard disk).

Before you can register this image, you will have to install Cederholm’s extension. Copy the file ‘Georeference Tiff image.avx’ from the folder ‘Extensions/Geo_ref’ on the CD to the folder of which the path ends like: ESRIAV_GIS30ARCVIEWEXT32 on your hard disk.

1. Open ArcView open a new View.

2. Set the working directory.

3. Set the projection to Category: ‘Projections of the World’, and Type: ‘Geographic’.

4. Activate via F ile/Extension. the ‘Georeference Tiff Image’ extension (Figure 18.2).

FIGURE 18.2
Georeference Tiff Image extension

After this activation, you will see that on the button bar a new green button has appeared, the colour of the button is rather diffuse (Figure 18.3).

FIGURE 18.3
The new Georeference Tiff image extension button

5. Add the image to the View. The image can be found in the file ‘Pais_pesca_Biodiv_floodplain.tif’, which is located in the folder ‘RW_05_Geo’. Remember that when you add this file to the View to change the data source type to ‘Image Data Source’ in the Add Theme window. Do you notice that the button of the Georeference Image extension has become brightly green? (Another point to remember: If you copy a file from a CD, its properties will be copied as well. In the case of a file on a CD, these are automatically set as read only. Remove this read only property as described on page 5.

The Georeference Tiff Image extension is straight forward. Two points on an image are selected, located diagonally and as far from each other as reasonably possible. Enter the geographic location of the two points and the map will be registered in ArcView. However, remember, ArcView works with Decimal Degrees (DD) and the map from the file contains the more commonly used Degrees, Minutes, and Seconds (DMS) system. Remember how to convert DMS into DD [43] (see Chapter: Points and degrees, page 35).

In the map we encircled the two points you need to use for this exercise (Figure 18.1), with the following coordinates:

TABLE 18.1
Georeferencing points of figure 18.1

6. Zoom in on Point 1, and press the Georeferencing Image button in the toolbar.

7. Your pointer is now changed into a little cross (the georeference pointer). Place this cross on the crossing of the longitude and latitude lines in the red circle, and select this point by clicking your mouse.

8. The message ‘Point collected’ appears, press OK.

9. Zoom to Full Extent, zoom in on Point 2, press the Georeferencing image button again, place the cross (georeference pointer) on the crossing of the longitude and latitude lines in the red circle, and select this point by clicking your mouse.

10. The message ‘Use current point to complete rectangle?’ appears, press Yes . 11. The Georeference Points window pops up, where you are requested to ‘Enter Map Coordinates’(Figure 18.4). Use the DD-coordinates from the table above and press OK .

FIGURE 18.4
Georeference Points window

12. After you have pressed OK , you will see that the View goes back to Full Extent. If you move the pointer over the View you will see that the coordinates (in the top right corner) are changing.

The image is now registered in this project. To check the registering process, open the Theme Properties, and you will see values for the Left, Right, Bottom and Top sides of the Theme (Figure 18.5).

FIGURE 18.5
Theme properties of the registered image

13. Save the project with a name of your own choice (menu bar: File/Save Project) You can check if the image was registered correctly by Adding the Theme ‘Flood_district.shp’ from the folder 㢰_calc_with_grids’ on your CD (remember that the Data Source Type of this Theme is Feature Data Source). After you have added this Theme, Zoom out, and you will see that the image fits perfectly into the ‘flood_district.shp’ Theme (Figure 18.6).

FIGURE 18.6
Two Themes, to check the registration process

14. Close the Project, exit ArcView.

With the exercise above you have learned how to register an image on which coordinates were clearly indicated. However, it is possible that coordinates are not clearly indicated on an image. You can still reference such an image as long as you have the exact coordinates of two locations on the map.

The longitude and latitude of the villages in the ‘Pais_pesca_Biodiv_floodplain.tif’ image were recorded with a GPS and are presented below, with the corresponding XY coordinates. The next exercise is to register the image with the GPS readings of the villages Vanda and Dura.

TABLE 18.2
Coordinates of villages in figure 18.1

1. Start ArcView, open a new View.

2. Set your Working directory.

3. Set the projection (Category: Projections of the World, Type: Geographic).

4. Activate the Georeference Tiff Image-extension ( File/Extensions. ).

5. Add the Image (Image data source) ‘Pais_pesca_Biodiv_floodplain.tif’ located in the folder ‘RW_05_Geo’ to the View.

6. Click the Georeference Image-button.

7. Place the georeference pointer in the centre of the village Vanda and click your mouse, a pop-up window shows again, with remark: ‘point collected’, press OK .

8. Place the georeference pointer in the centre of the village Dura and click your mouse. You will be asked: ‘Use current point to complete rectangle’, press Yes .

9. Enter the Map Coordinates of Vanda and Dura and press OK .

10. Add the Theme ‘Flood_district.shp’ from the folder 㢰_calc_with_grids’ on your CD again and check if you registered the image correctly.

11. Take care: the XY coordinates are not attached to the Theme. If you open the open the Image Theme in a new project, you will have to repeat the above described procedures. However, if you save the project, the coordinates will be added automatically when the project is opened.

12. In order to attach the coordinates to the image file you need to convert the image into a World tiff-file. This is easily done with a small script ‘make a world file.avx’ also written and provided by Mark Cederholm.

13. Copy the file ‘Make a world file.avx’ from the folder /Extensions/Geo_ref/ on the CD to the folder of which the end of the path looks like: ‘ESRIAV_GIS30ARCVIEWEXT32’on your hard disk.

14. In the F ile/Extensions menu activate the Make a World file extension.

15. On the button bar, the Make World File button appears.

16. Make the to be converted Image active (remember: pressing the legend of the Theme once, so that it looks like is a bit elevated). In this case the ‘Pais_pesca_Biodiv_floodplain.tif’.

17. Click the Make World File button.

18. A message ‘Done’ will pop up, press OK .

The extension has created a new small file with the extension ‘.Tifw’ where information on the coordinates is stored. As long as the original ‘.tif’ file and this ‘.Tifw’ file are located in the same folder, opening the ‘.tif’ file in a new View will result in an Image which is registered.

19. Close the Project and exit ArcView.

You can now proceed to digitize the data on the referenced map as explained below.

18.1.2. On screen digitizing in ArcView [44]

From your registered Image you can now abstract any data you need by digitising in ArcView. The newly created themes will automatically encompass the geo reference and can be further used to analyse your data in ArcView. Digitising is in principle nothing more then accurately tracing lines, polygons, points, etc. on your screen and saving them as a new theme.

The next exercise will show you how to do this by digitising the water bodies in the floodplain of Pais Pesca

1. Open ArcView and a new View.

2. Set your working directory.

3. Set the projection (Category: Projections of the World, Type: Geographic).

4. Add the image ‘Pais_pesca_Biodiv_floodplain.tif’ from the folder you used to work from in the previous exercise.

5. Check if the image is registered in the upper right corner you should see values around 90 and 23. You could also check the Theme Properties table, menu: T heme/ P roperties. , (see Figure 17.5).

If the image is not registered, go back to the previous exercise (where you used the ‘Make a world file’ extension), and try again to make the file with the ‘.Tifw’ extension. Make sure you save this file in the same folder as where the Image file is located that you use in this digitizing exercise.

6. Go via the menu bar to V iew/ N ew Theme. . In the dialogue box select ‘Polygon’ as the Feature type from the drop-down menu (Figure 18.7) and press OK .

FIGURE 18.7
The New Theme dialogue box

7. In the next dialogue box specify where to save the new Theme. Give the Theme a name, for the exercise you do now you can use the name: ‘Water_bodies.shp’. press OK .

A new Theme is now added to the View, which at the moment contains no features. The check box for the ‘Water_bodies.shp’ Theme in the View’s Table Of Content has a dashed line around it, indicating that the Theme can be edited (Figure 18.8).

FIGURE 18.8
Table Of Content indicating which Theme can be edited

8. From the button bar select the Drawing Tool button. The Drawing Tool button is the tool button in the second row, far right (Figure 18.9).

FIGURE 18.9
The Drawing Tool button

9. Move the pointer to the Drawing Tool button, click and hold the mouse button to bring up a drop down list of the Drawing Tool buttons. Select the Draw Polygon button.

10. Zoom in on the water body you want to digitize (for this exercise select Jugini). After you have zoomed in on the waterbody, and made sure you have pressed the Draw Polygon button, you need to trace the outer borders of the waterbody with the pointer. Place the pointer on the position where you wish to start and click once. Move the pointer over the outer border, clicking every point where you change direction (Figure 18.10). After you have followed the outline of the waterbody, and have returned to your starting point, make a double-click, indicating you have finished the polygon (Figure 18.11).

FIGURE 18.10
Following the outline of a waterbody

FIGURE 18.11
The completed polygon. See the eight squares around it, indicating it is active

11. If you have made a mistake you can always delete the polygon you made and start over: After you completed the polygon, you will see there are eight squares around it. This means the polygon is active, and you can delete it. If these squares are not around the polygon, you first have to activate it. Select the polygon by using the Pointer tool . First press the Pointer tool button in the button bar, then move the pointer to the polygon you want to delete, and click it (you will see eight squares appear around the polygon, meaning the polygon is active). When you now press the delete button on your keyboard, you will delete the polygon.

The first digitising is always a little bid rough and the shape you made never overlaps nicely with the original shape. So you have to continue editing.

12. To be able to edit the shape of the new polygon, you need to see the outline only: Go to the Legend Editor of the new Theme (either by double clicking on the legend, or through T heme/Edit L egend. ), double click on the Symbol, the fill palette comes up, click on the colour button, select in the Color dropdown box ‘Foreground’, select the square with the cross, indicating you do not want any colour as foreground. After this, select in the Color dropdown box ‘Outline’, select a colour for the outline (in this example red is used). Close the Fill Palette, press Apply in the Legend Editor and close the Legend Editor (Figure 18.12).

FIGURE 18.12
The polygon, after changing the Legend to show only the outline in red

FIGURE 18.13
The polygon, after pressing the Vertex Edit tool button

13. Click on the Vertex Edit tool button and you will see the outline of your Theme with small squares indicating the places were you clicked your mouse (Figure 18.13).

14. Zoom in on a part of the outline of the Polygon and press the Vertex Edit tool button again. You can now adjust the shape of the polygon (as explained below) until its outline follows the shape of the original image.

15. Move the pointer onto one of the squares, you will see that once you are on the square the pointer changes from an arrow into a cross. As soon as the pointer is changed into a cross, you can move the square on the outline by clicking your mouse keeping the mouse button pressed, moving the square to the spot where you want it and release the mouse button (‘drag and drop’). By doing this you can change the shape of the polygon.

If you need more points where the direction of the outline changes, you can add a square to the outline. Move the pointer to the outline, and you will see that it changes from an arrow to a circle with a cross inside. Once you see this, you can click your mouse. By doing so, you have added a square which you can move around as described above.

Once you are satisfied with the outline of the digitized water body, use the same procedure to digitize the other water bodies of the Image, each time starting with using the Draw Polygon tool button.

16. After you have digitized all polygons satisfactory, the editing session should be ended: Go via the menu bar to T heme/Stop E diting . A message will pop up asking you ‘Save edits to Water_bodies.shp’ Press Yes to save your edits.

17. Open the Attributes table of the ‘Water_bodies.shp’ Theme. You will see only two columns. You might want to add a field to add the names of the different waterbodies. In Chapter ‘How to edit your Theme table and do calculations?’ on page 31 is explained how to add a field to an attribute table. One problem rests though, and that is how to know which record belongs to which polygon. As explained in the same chapter, you can edit the table by adding a field for the names of the different water bodies. Now add per record a number (Figure 18.14).

FIGURE 18.14
The Attributes table of the New Theme, with the added field ‘WaterBodyNumber’

Now you have to determine which record belongs to which polygon: Press the Select tool button on the button bar.

18. Select one record in the attribute table by moving the pointer to a record, click your mouse. You will notice that the record you selected becomes yellow.

FIGURE 18.15
The Attributes table with a record selected

19. Close the attribute table of ‘Water_bodies.shp’. Now you will see in the View that one of the polygons is coloured yellow. This is the polygon that you selected in the attributes table. So now you know the number of the polygon, and you know which polygon in the View is connected to that record. The name of the waterbody is shown in the ‘Pais_pesca_Biodiv_floodplain.tiff’ Image.

20. Determine the Names of all waterbodies, and add them to the Attributes of Water_bodies.shp table.

After making this Theme, containing polygons representing waterbodies, another Theme to be made could be a Theme showing the positions of the roads. What would be the feature type of this Theme? [45]

1. Open the View of your current project.

2. Go via the menu bar to V iew/ N ew Theme. . In the dialogue box select ‘Line’ as the Feature type from the drop-down menu and press OK .

3. Save the New Theme with the filename: ‘road.shp’ in the dialogue box that pops up. Make sure you save it in the appropriate folder and on the appropriate drive.

The ‘road.shp’ Theme is added to the map legend, with the dashed line round its checkbox indicating that the Theme is editable. Go to the Legend Editor (menu bar: T heme/Edit L egend. ) and change the Symbol in such a way that it will be clearly visible (change the colour into purple, and width to 2).

4. From the Drawing tool button drop-down list select the Draw Line tool button.

5. Choose a road in the Image you want to digitize, and follow the same procedure as you have followed digitizing a polygon Click along the length of the road, until you have reach the end of the road, end the line by double-clicking. Select the line tool again and digitize the next part of the roads, etc.

After you have digitized all roads that you wanted to digitize, save your work: Via the menu bar go to: T heme/Stop E diting . Click Yes to save edits.

The line that you worked on last is highlighted (with yellow). Click the Clear Selected Features tool button to deselect it.

After a polygon Theme, and a line Theme, you will make in the following exercise a point Theme. This type of Theme is the ideal way to store data, if you have made specific measurements on a specific location. For this exercise you will digitize the locations of the sampling stations on the ‘Pais_pesca_Biodiv_floodplain.tif’ Image.

6. Open the View of your current project.

7. Go via the menu bar to View/New Theme. . In the dialogue box select ‘Point’ as the Feature type from the drop-down menu and press OK .

8. Save the New Theme with the filename: ‘Sampling_stations.shp’ in the dialogue box that pops up. Make sure you save it in the appropriate folder and on the appropriate drive.

The ‘Sampling_stations.shp’ Theme is added to the map legend, with the dashed line round its checkbox indicating that the Theme is editable.

9. From the Drawing tool button drop-down list select the Draw Point tool button.

10. Zoom in on the sampling station you want to digitize (sampling stations are indicated by a red circle with a red cross). After you have zoomed in on the sampling station, and made sure you have pressed the Draw Polygon button, move the pointer to the centre of the sampling station, and click it with your mouse. After having selected one sampling station, go to the next and select if following this procedure.

11. When you have added all sampling stations you can finish editing the new Theme go via the menu bar to T heme/Stop E diting . Press Yes to save the edits.

After these exercises you know how to register an image of which at least two locations are georeferenced, and you can use this image to make three types of Themes (polygon, line, and point Themes). Previously (Putting field data into a GIS, page 36) you have learned how to add data to existing Themes, which you practised in these exercises as well.

18.1.3. Real world exercise, Rio Uneiuxi, Brazil

If you want to get some more experience registering, then try to register the map of the Rio Uneixu [46] in Brazil.

1. Start ArcView, Open a new View

2. Set the working directory

3. Set the projection (Category: Projections of the World, Type: Geographic) and

4. Add the image ‘Rio_Uneiuxi.tif’ located in the folder ‘RW_06_Brasil_Geo’.

You will see the scan of a map of an area in the north of Brazil. Unfortunately, the latitudes and longitudes of the lines you see in the image, which are normally written on the map, are not known. The only data that came with this map is that one point in the red squares (Figure 18.16) has the values: X: -65.06941 and Y:-0.72384, while the other point the following data contains: Latitude: 0䓛ָ” (S), Longitude: 65䓏ֽ” (W) (The letters S and W mean south and west).

FIGURE 18.16
Image of Rio Uneiux with georeferenced points

Recalculate the obvious DMS values of the second point into decimal degrees. Pay special attention to the fact that the points are south of the Equator, and west of the Greenwich Meridian [47] .

Determine which values belong to which point, keep in mind the fact that the points are south of the Equator, and west of the Greenwich Meridian [48] .

5. With the (calculated) values of both points, register the image, using these points, in ArcView.


Voir la vidéo: QGIS: counting points in polygons