6
votes

Tracer une carte à l'aide de géopandas et de matplotlib

J'ai un petit csv qui a 6 coordonnées de Birmingham en Angleterre. J'ai lu le csv avec des pandas, puis je l'ai transformé en GeoPandas DataFrame en changeant mes colonnes de latitude et de longitude avec Shapely Points. J'essaie maintenant de tracer mon GeoDataframe et tout ce que je peux voir, ce sont les points. Comment puis-je également représenter la carte de Birmingham? Une bonne source de documentation sur GeoPandas serait également très appréciée.

from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

df = pd.read_csv('SiteLocation.csv')
df['Coordinates'] = list(zip(df.LONG, df.LAT))
df['Coordinates'] = df['Coordinates'].apply(Point)
# Building the GeoDataframe 
geo_df = gpd.GeoDataFrame(df, geometry='Coordinates')
geo_df.plot()  


2 commentaires

@Kanak pouvez-vous ajouter cela comme réponse?


@joris Ajouté en tant que réponse.


3 Réponses :


0
votes

Essayez df.unary_union . La fonction agrégera les points en une seule géométrie. Jupyter Notebook peut le tracer


1 commentaires

Désolé mon pote - ce n'était pas ce que je cherchais. Mais merci . Vous regardez d'autres questions ouvertes



10
votes

La documentation GeoPandas contient un exemple sur la façon d'ajouter un arrière-plan à une carte ( https://geopandas.readthedocs.io/en/latest/gallery/plotting_basemap_background.html ), qui est expliqué plus en détail ci-dessous.


Vous devrez traiter tuiles , qui sont des images (png) servies via un serveur Web, avec une URL comme

http: //.../Z/X/Y.png , où Z est le niveau de zoom, et X et Y identifient la tuile

Et la documentation de geopandas montre comment définir des tuiles comme arrière-plans pour vos tracés, en récupérant les bons et en faisant tout le travail autrement difficile de synchronisation spatiale, etc ...

p >


Installation

En supposant que GeoPandas est déjà installé, vous avez besoin du contextily en plus. Si vous êtes sous Windows, vous pouvez jeter un œil à Comment installer Contextily?

Cas d'utilisation

Créez un script python et définissez les fonction d'aide contextuelle

import matplotlib.pyplot as plt
from shapely.geometry import Point
import geopandas as gpd
import pandas as pd

# Let's define our raw data, whose epsg is 4326
df = pd.DataFrame({
    'LAT'  :[-22.266415, -20.684157],
    'LONG' :[166.452764, 164.956089],
})
df['coords'] = list(zip(df.LONG, df.LAT))

# ... turn them into geodataframe, and convert our
# epsg into 3857, since web map tiles are typically
# provided as such.
geo_df = gpd.GeoDataFrame(
    df, crs  ={'init': 'epsg:4326'},
    geometry = df['coords'].apply(Point)
).to_crs(epsg=3857)

# ... and make the plot
ax = geo_df.plot(
    figsize= (5, 5),
    alpha  = 1
)
add_basemap(ax, zoom=10)
ax.set_axis_off()
plt.title('Kaledonia : From Hienghène to Nouméa')
plt.show()

et jouez

import contextily as ctx

def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

entrez la description de l'image ici


Remarque: vous peut jouer avec le zoom pour trouver la bonne résolution de la carte.

Par exemple / Ie :

 enter image description here

... et de telles résolutions appellent implicitement à changer les limites x / y.


8 commentaires

Une question @ Herc01?


Merci d'avoir répondu. Que voulez-vous dire par la dernière remarque sur le niveau de zoom? Normalement, vous ne devriez pas avoir besoin de modifier manuellement les limites x / y si vous changez le niveau de zoom.


@joris, ce que je voulais dire, c'est que l'utilisation de la valeur définie pour le paramètre zoom équivaut en quelque sorte à changer la résolution (image) de la carte. Êtes-vous d'accord avec cela?


Oui, c'est exactement ce que signifie le "zoom". Le nom peut être un peu déroutant car la carte elle-même n'est pas agrandie (en regardant une zone plus petite), mais un zoom plus bas / plus haut zoom pour la zone entière et résultant ainsi en une résolution inférieure / supérieure .


@joris Êtes-vous sous Windows? il semble que vous ignorez à quel point il est difficile de l'installer en contexte. Cette difficulté est interdite pour la plupart des utilisateurs de Windows. D'où la section qui expliquait comment (enfin réussir) une installation contextuelle. Faire pip3 installer contextuellement sous Windows ne fonctionne tout simplement pas.


désolé, j'avais l'intention d'ajouter un commentaire dans l'édition pourquoi je l'ai supprimé, mais j'ai oublié cela: je sais que cela peut certainement être difficile à installer, mais je pense qu'ici cela détourne de la réponse réelle (et ce n'était pas non plus la seule possible façon d'installer, par exemple pour quelqu'un utilisant conda qui aurait été la mauvaise façon). Souhaitez-vous poser une nouvelle question "Comment installer en contexte?", Et y ajouter votre réponse sur cette partie. Et puis nous pouvons établir un lien vers cette question à partir d'ici (et garder cette réponse concentrée sur le comment).


@Kanak; oui j'ai une question sur l'epsg; où dois-je trouver. J'ai trouvé une carte et je voudrais faire de même et qu'entendez-vous par convertir puisque notre epsg est 3857.?


@ Herc01 Il y a l'epsg de vos points, et celui de la carte sur laquelle vous voulez projeter vos points, et les deux doivent être identiques. Je ne sais pas à quel point vous êtes familier avec les trucs géographiques, mais c'est comme si vous aviez un plan ax / y dont les unités étaient exprimées en cm d'une part, et quelques points dont les coordonnées étaient en inch d'autre part. L'epsg de la carte utilisée dans l'exemple est 3857. Si vous souhaitez utiliser la même carte, il vous suffit de déterminer quel est l'epsg de vos points. Il n'y a pas de moyen systématique de le faire AFAIK. Avez-vous des paires de lat-long à me montrer?



0
votes

Je veux juste ajouter le cas d'utilisation concernant le zoom dans lequel le fond de carte est mis à jour en fonction des nouvelles coordonnées xlim et ylim . Une solution que j'ai trouvée est:

  • Commencez par définir les rappels sur ax qui peuvent détecter xlim_changed et ylim_changed
  • Une fois que les deux ont été détectés comme modifiés, obtenez la nouvelle plot_area appelant ax.get_xlim () et ax.get_ylim () li >
  • Effacez ensuite l ' axe et retracez le fond de carte et toutes les autres données

Exemple de carte du monde montrant les capitales. Vous remarquez que lorsque vous effectuez un zoom avant, la résolution de la carte est en cours de mise à jour.

affine==2.3.0
attrs==19.3.0
autopep8==1.4.4
Cartopy==0.17.0
certifi==2019.11.28
chardet==3.0.4
Click==7.0
click-plugins==1.1.1
cligj==0.5.0
contextily==1.0rc2
cycler==0.10.0
descartes==1.1.0
Fiona==1.8.11
geographiclib==1.50
geopandas==0.6.2
geopy==1.20.0
idna==2.8
joblib==0.14.0
kiwisolver==1.1.0
matplotlib==3.1.2
mercantile==1.1.2
more-itertools==8.0.0
munch==2.5.0
numpy==1.17.4
packaging==19.2
pandas==0.25.3
Pillow==6.2.1
pluggy==0.13.1
py==1.8.0
pycodestyle==2.5.0
pyparsing==2.4.5
pyproj==2.4.1
pyshp==2.1.0
pytest==5.3.1
python-dateutil==2.8.1
pytz==2019.3
rasterio==1.1.1
requests==2.22.0
Rtree==0.9.1
Shapely==1.6.4.post2
six==1.13.0
snuggs==1.4.7
urllib3==1.25.7
wcwidth==0.1.7

Fonctionne sous Linux Python3.8 avec les installations pip suivantes

import geopandas as gpd
import matplotlib.pyplot as plt
import contextily as ctx


figsize = (12, 10)
osm_url = 'http://tile.stamen.com/terrain/{z}/{x}/{y}.png'
EPSG_OSM = 3857
EPSG_WGS84 = 4326

class MapTools:
    def __init__(self):
        self.cities = gpd.read_file(
            gpd.datasets.get_path('naturalearth_cities'))
        self.cities.crs = EPSG_WGS84
        self.cities = self.convert_to_osm(self.cities)

        self.fig, self.ax = plt.subplots(nrows=1, ncols=1, figsize=figsize)
        self.callbacks_connect()

        # get extent of the map for all cities
        self.cities.plot(ax=self.ax)
        self.plot_area = self.ax.axis()

    def convert_to_osm(self, df):
        return df.to_crs(epsg=EPSG_OSM)

    def callbacks_connect(self):
        self.zoomcallx = self.ax.callbacks.connect(
            'xlim_changed', self.on_limx_change)
        self.zoomcally = self.ax.callbacks.connect(
            'ylim_changed', self.on_limy_change)

        self.x_called = False
        self.y_called = False

    def callbacks_disconnect(self):
        self.ax.callbacks.disconnect(self.zoomcallx)
        self.ax.callbacks.disconnect(self.zoomcally)

    def on_limx_change(self, _):
        self.x_called = True
        if self.y_called:
            self.on_lim_change()

    def on_limy_change(self, _):
        self.y_called = True
        if self.x_called:
            self.on_lim_change()

    def on_lim_change(self):
        xlim = self.ax.get_xlim()
        ylim = self.ax.get_ylim()
        self.plot_area = (*xlim, *ylim)
        self.blit_map()

    def add_base_map_osm(self):
        if abs(self.plot_area[1] - self.plot_area[0]) < 100:
            zoom = 13

        else:
            zoom = 'auto'

        try:
            basemap, extent = ctx.bounds2img(
                self.plot_area[0], self.plot_area[2],
                self.plot_area[1], self.plot_area[3],
                zoom=zoom,
                url=osm_url,)
            self.ax.imshow(basemap, extent=extent, interpolation='bilinear')

        except Exception as e:
            print(f'unable to load map: {e}')

    def blit_map(self):
        self.ax.cla()
        self.callbacks_disconnect()
        cities = self.cities.cx[
            self.plot_area[0]:self.plot_area[1],
            self.plot_area[2]:self.plot_area[3]]
        cities.plot(ax=self.ax, color='red', markersize=3)

        print('*'*80)
        print(self.plot_area)
        print(f'{len(cities)} cities in plot area')

        self.add_base_map_osm()
        self.callbacks_connect()

    @staticmethod
    def show():
        plt.show()


def main():
    map_tools = MapTools()
    map_tools.show()

if __name__ == '__main__':
    main()


0 commentaires