Principes de base de Python x SIG (1)

Cette fois, j'utiliserai Automatiser les processus SIG fourni par l'Université d'Helsinki comme matériel d'étude. Il se compose de 7 semaines, et vous pouvez en apprendre plus sur la visualisation en utilisant GeoPnadas et matplotlib. De plus, des devoirs (questions pratiques) pour les étudiants sont joints chaque semaine, et cet article fournit principalement des exemples de réponses. Partie 1: Semaine1-Semaine2 Partie 2: Semaine3-Semaine4 Partie 3: Semaine5-Semaine6 (Semaine7 n'a pas de problèmes) Partie 4 (en construction): Affectation finale

Ce qui suit fait partie du contenu du matériel didactique hebdomadaire.

Week 1

1-1 Création d'une fonction pour créer des points, des lignes et des polygones

En tant que condition du problème, certaines conditions sont requises pour que la fonction fonctionne correctement.

from shapely.geometry import Point, LineString, Polygon

#point
def create_point_geom(x, y):
    point = Point(x, y)
    return point

#Prend une liste de lignes comme argument
def create_line_geom(points):
    assert type(points)=="list", "Input should be a list!"
    assert len(points)>=2, "LineString object requires at least two Points!"
    line = LineString([points[0], points[1]])
    return line

#Prend une liste de polygones comme argument
def create_poly_geom(coords):
    assert type(coords) is list, "Input should be a list!"
    assert len(coords)>=3, "Polygon object requires at least three Points!"
    for i in coords:
        assert type(i) is tuple, "All list values should be coordinate tuples!"
    poly = Polygon(coords)
    return poly

1-2 Acquisition du centre de gravité, de la surface du polygone et de la distance

L'instruction assert n'est pas essentielle, certaines parties sont donc omises.

#L'argument du centre de gravité est un point, une ligne ou un polygone
def get_centroid(gem):
    assert type(gem) == 'shapely',  "Input should be a Shapely geometry!"
    assert gem.geom_type in ['Point', 'LineString', 'Polygon'], "Input should be a Shapely geometry!"
    centroid = gem.centroid
    return centroid
#zone
def get_area(poly):
    return poly.area
#distance
def get_length(geom):
    if geom.geom_type == 'LineString':
        return geom.length
    elif geom.geom_type == 'Polygon':
        return geom.exterior.length

1-3 Acquisition et affichage de données textuelles

Il s'agit d'un simple exercice de traitement de données utilisant des données mobiles d'Helsinki. from_x, from_y indiquent le point de départ, to_x, to_y indiquent le point d'arrivée, et le but est d'obtenir la distance totale parcourue à partir de ces données.

import pandas as pd

#Lecture des données
data = pd.read_table("data/travelTimes_2015_Helsinki.txt", sep=";",)
data = data[['from_x','from_y', 'to_x', 'to_y', 'total_route_time',]]
#Création de données de coordonnées
orig_points = []
dest_points = []
from shapely.geometry import Point
for index, row in data.iterrows():
    orig = Point(row['from_x'], row['from_y'])
    dest = Point(row['to_x'], row['to_y'])
    orig_points.append(orig)
    dest_points.append(dest)
#Créer une ligne mobile
from shapely.geometry import LineString
lines = []
for orig, dest in zip(orig_points, dest_points):
    line = LineString([orig, dest])
    lines.append(line)
#Obtenez la distance totale parcourue
total_length = 0
for line in lines:
    total_length += line.length

Week2 Apprenez une bibliothèque appelée Geopandas, une extension de Pandas. Comme les pandas, il est très facile à utiliser et la manipulation de l'espace est facile.

2-1 Conversion des données de coordonnées en GeoDataFrame

import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

#Coordonner les données à utiliser
longitudes = [29.99671173095703, 31.58196258544922, 27.738052368164062, 26.50013542175293, 26.652359008789062, 25.921663284301758, 22.90027618408203, 23.257217407226562,
           23.335693359375, 22.87444305419922, 23.08465003967285, 22.565473556518555, 21.452774047851562, 21.66388702392578, 21.065969467163086, 21.67659568786621,
           21.496871948242188, 22.339998245239258, 22.288192749023438, 24.539581298828125, 25.444232940673828, 25.303749084472656, 24.669166564941406, 24.689163208007812,
           24.174999237060547, 23.68471908569336, 24.000761032104492, 23.57332992553711, 23.76513671875, 23.430830001831055, 23.6597900390625, 20.580928802490234, 21.320831298828125,
           22.398330688476562, 23.97638702392578, 24.934917449951172, 25.7611083984375, 25.95930290222168, 26.476804733276367, 27.91069221496582, 29.1027774810791, 29.29846954345703,
           28.4355525970459, 28.817358016967773, 28.459857940673828, 30.028610229492188, 29.075136184692383, 30.13492774963379, 29.818885803222656, 29.640830993652344, 30.57735824584961,
           29.99671173095703]
latitudes = [63.748023986816406, 62.90789794921875, 60.511383056640625, 60.44499588012695, 60.646385192871094, 60.243743896484375, 59.806800842285156, 59.91944122314453,
           60.02395248413086, 60.14555358886719, 60.3452033996582, 60.211936950683594, 60.56249237060547, 61.54027557373047, 62.59798049926758, 63.02013397216797,
           63.20353698730469, 63.27652359008789, 63.525691986083984, 64.79915618896484, 64.9533920288086, 65.51513671875, 65.65470886230469, 65.89610290527344, 65.79151916503906,
           66.26332092285156, 66.80228424072266, 67.1570053100586, 67.4168701171875, 67.47978210449219, 67.94589233398438, 69.060302734375, 69.32611083984375, 68.71110534667969,
           68.83248901367188, 68.580810546875, 68.98916625976562, 69.68568420410156, 69.9363784790039, 70.08860778808594, 69.70597076416016, 69.48533630371094, 68.90263366699219,
           68.84700012207031, 68.53485107421875, 67.69471740722656, 66.90360260009766, 65.70887756347656, 65.6533203125, 64.92096710205078, 64.22373962402344, 63.748023986816406]
#Polygonisation
coordpairs = list(zip(longitudes, latitudes))
poly = Polygon(coordpairs)
#Créer un GeoDataFrame
geo = gpd.GeoDataFrame(index=[0], columns=['geometry'])
geo['geometry'] = poly
#Illustré
import matplotlib.pyplot as plt
geo.plot()
#Enregistrer le fichier shp
fp = 'polygon.shp'
geo.to_file(fp)

2-2 Tracé des données de coordonnées stockées dans csv

Le CRS haineux est là. La fonction Apply fonctionne bien avec les pandas et peut être traitée avec un code court.

#Bibliothèque
import pandas as pd
from shapely.geometry import Point, LineString, Polygon
from pyproj import CRS
import matplotlib.pyplot as plt
#lire csv
data = pd.read_csv('data/some_posts.csv')
#Création de données ponctuelles
data = pd.read_csv('data/some_posts.csv')
make_point = lambda row:Point(row['lat'],row['lon'])
data['geometry'] = data.apply(make_point, axis=1)

#La conversion en géométrie GeoDataFrame et système de coordonnées doit être spécifiée.
geo = gpd.GeoDataFrame(data, geometry='geometry',crs=CRS.from_epsg(4326).to_wkt())

#terrain
geo.plot()

2-3 Calcul de la distance parcourue

Calculez la distance parcourue par chaque utilisateur à partir des données de coordonnées affichées sur SNS. La conversion du système de coordonnées apparaîtra, mais notez les points suivants. CRS est une nuisance partout. .. ..

-Parce que les données GeoDataFrame ne sont pas converties uniquement en définissant le système de coordonnées (ex. Data.crs = CRS.from_espg (4276)), la conversion du système de coordonnées (ex. Data = data.to_crs (epsg = 4276)) est nécessaire. -Si le système de coordonnées de GeoDataFrame n'est pas défini, définissez le système de coordonnées avant de convertir le système de coordonnées.

import geopandas as gpd
import pandas as pd
from pyproj import CRS
from shapely.geometry import Point, LineString, Polygon
#Lecture des données
data = gpd.read_file('Kruger_posts.shp')
#Conversion du système de coordonnées
data = data.to_crs(epsg=32735)
#Créer une ligne mobile en classant par identifiant(S'il n'y a qu'un seul article, il n'a pas été déplacé,Entrez une valeur Aucun)
grouped = data.groupby('userid')
movements = gpd.GeoDataFrame(columns=['userid', 'geometry'])
for key, group in grouped:
    group = group.sort_values('timestamp')
    if len(group['geometry'])>=2:
        line = (LineString(list(group['geometry'])))
    else:
        line=None
    movements.at[count, 'userid'] = key
    movements.at[count, 'geometry'] = line
movements.crs = CRS.from_epsg(32735)
#Calcul de la distance parcourue
def cal_distance(x):
    if x['geometry'] is None:
        return None
    else:
        return x['geometry'].length
movements['distance'] = movements.apply(cal_distance, axis=1)
#Distance parcourue moyenne,Valeur maximum,valeur minimum
print(mean(movements['distance'])) #138871.14194459998
print(max(movements['distance'].dropna())) #8457917.497356484
print(min(movements['distance'])) #0.0