Geopandas.ipynb

Kernel: Python 3 (ipykernel)

Geopandas y operaciones GIS

@author: Roberto mendoza

import pandas as pd
from pandas import Series, DataFrame
import numpy as np
import matplotlib.pyplot as plt 

import geopandas as gpd  # manejo de datos georefereciados
from geopandas import GeoSeries # series de datos georerenciados
from shapely.geometry import Point, LineString, Polygon, MultiLineString # objetos geométricos
from shapely.ops import nearest_points  # operaciones entre objetods geometricos
import contextily as cx  # Fondo Goole maps, fondo de mapa 
from pyproj import CRS, Geod # proyecciones a sistemas planares


import matplotlib.patches as mpatches
import haversine as hs # distancia de grat-cricle entre puntos
from  geopy import distance # distancia entre puntos 
from tqdm import tqdm # contador de tiempo en un loop
from matplotlib.lines import Line2D


import warnings
warnings.filterwarnings('ignore') # eliminar warning messages 

no module name, solo instalelo usando !pip install

1.0 Polygonos

# Load .gdb dataset.
# La particuladridad es que almacena basese de datos por capas. Podemos acceder a cada una de ellas con layer
# El layer 5 contiene el shapefile de polygono de cada distrito 

dist_mita = gpd.read_file(r'..\data\Mita\mita.gdb', layer=5
                         )
dist_mita

dist_mita.crs

# sistema de coordenadas proyectadas en un plano cartesiano (unidad de medida en metros) bajo el método Transverse Mercator
# compatible para países de Argentina, Brazil, Chile, Colombia, Ecuadro y Peru

<Derived Projected CRS: EPSG:32718>
Name: WGS 84 / UTM zone 18S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 78°W and 72°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Brazil. Chile. Colombia. Ecuador. Peru.
- bounds: (-78.0, -80.0, -72.0, 0.0)
Coordinate Operation:
- name: UTM zone 18S
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Sistema de coordendas proyectada usada en la Mita

link

# Conversión de sistemas
# del sistema EPSG:32718 hacia el sistema de coordenadas estándar EPSG:4326 (WGS84)

# El WGS 84 es un sistema geodésico de coordenadas geográficas usado mundialmente,
# que permite localizar cualquier punto de la Tierra por medio de tres unidades dadas.


dist_mita.to_crs(epsg=4326,inplace=True)
dist_mita

dist_mita.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# Proyección cilindrica equidistante 
    
mita_proy_cil = dist_mita.to_crs(epsg=4087)
mita_proy_cil

mita_proy_cil.crs

<Derived Projected CRS: EPSG:4087>
Name: WGS 84 / World Equidistant Cylindrical
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: World Equidistant Cylindrical
- method: Equidistant Cylindrical
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# Se retorna al sistema de coordenadas iniciales

mita_cartesiano = dist_mita.to_crs(epsg=32718)
mita_cartesiano

mita_cartesiano.plot()

<AxesSubplot:>

2.0 Point

# Coordenad de las capitales de los distritos dentro de la MITA minera

centroids = gpd.read_file(r'..\data\Mita\mita.gdb', layer=0
                         )
centroids 

centroids.crs

<Derived Projected CRS: PROJCS["mita_equi",GEOGCS["WGS 84",DATUM["WGS_1984 ...>
Name: mita_equi
Axis Info [cartesian]:
- [east]: Easting (metre)
- [north]: Northing (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unnamed
- method: Equidistant Cylindrical
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# De columnas de latitud longitud al objeto geometry (Point)

centroids['Point'] = gpd.points_from_xy(centroids.LON, centroids.LAT, crs="EPSG:4326")
centroids

centroids.plot()

<AxesSubplot:>

# Del objeto Geometry (proyección cilindrico equidistante) a columnas latitud y longitud por separado 

centroids["longitude"] = centroids.geometry.map(lambda p: p.x)
centroids["latitude"] = centroids.geometry.map(lambda p: p.y)

centroids

3.0 Linestring

# Load a shapefile in Python 

mita_boundary = gpd.read_file(r'..\data\Mita\MitaBoundary.shp')

mita_boundary 

mita_boundary.plot()

<AxesSubplot:>

mita_boundary.crs

<Derived Projected CRS: EPSG:32718>
Name: WGS 84 / UTM zone 18S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 78°W and 72°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Brazil. Chile. Colombia. Ecuador. Peru.
- bounds: (-78.0, -80.0, -72.0, 0.0)
Coordinate Operation:
- name: UTM zone 18S
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

# load shapefile a nivel distrital

distritos = gpd.read_file(r'..\data\geopandas\LIMITE_DISTRITO\LIMITE_DIST.shp')

distritos

# centroide de polygonos. En este caso del poligono a nivel distrital 
distritos['centroid'] = distritos.geometry.centroid

distritos

3.0 Merge Geospatial information

# select columns and then convert to geodataframe. 

gdf = gpd.GeoDataFrame(
    centroids[['Point','UBIGEO']], geometry= centroids.Point)

gdf.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 305 entries, 0 to 304
Data columns (total 3 columns):
 #   Column    Non-Null Count  Dtype   
---  ------    --------------  -----   
 0   Point     305 non-null    geometry
 1   UBIGEO    305 non-null    float64 
 2   geometry  305 non-null    geometry
dtypes: float64(1), geometry(2)
memory usage: 7.3 KB

3.1 Join - Contains

# Merge: el poligono del distrito que contenga (contains) el centroide de los distritos ubicadas en la MITA minera

merge = gpd.sjoin(distritos, gdf , how="left", op="contains")
merge[merge.Point != None ] 

3.2 Join - whitin

# el centroide de los distritos ubicadas en la MITA minera está al interior del polygono del distrito (whitin) 

merge = gpd.sjoin( gdf, distritos , how="left", op="within")
merge[merge.Point != None ] 

3.3 Join - Intersects

# Merge: el poligono del distrito que intersecte (intersects) el centroide de los distritos ubicadas en la MITA minera

data_geo= gpd.sjoin(distritos, gdf , how="inner", op="intersects").reset_index()

f, ax = plt.subplots(figsize=(8,8))

data_geo['geometry'].plot(color='none', edgecolor='black', zorder=0.5, ax = ax)

data_geo['centroid'].plot(color = 'r', markersize=10, ax = ax)

<AxesSubplot:>

f, ax = plt.subplots(figsize=(8,8))

data_geo['geometry'].plot(color='none', edgecolor='black', zorder=0.5, ax = ax)

data_geo['Point'].plot(color = 'r', markersize=10, ax = ax)

mita_boundary.plot(color = 'b', markersize=10, ax = ax) 

<AxesSubplot:>

# Se proyecta al sistema convencional de coordenadas EPSG:4326
# D esta manera sea compatible con los otros mapas 

mita_boundary.to_crs(epsg=4326,inplace=True)

f, ax = plt.subplots(figsize=(15,15))

# capa de limite distrital

data_geo['geometry'].plot(color='none', linewidth=0.8,edgecolor='black', zorder=0.5, ax = ax)

# capital de cada distrito
data_geo['Point'].plot(color = '#DC143C', markersize=14, ax = ax, marker = "o", label="Distric's capital")

# Una parte de los limites de la MITA 
mita_boundary.plot(color = '#834333' ,linewidth=4, ax = ax, label='Mita Boundary target') 

# añade un mapa de fondo 

cx.add_basemap(ax, crs="EPSG:4326", attribution = False, zoom = 11)

# sin valores en los ejes 

plt.xticks([])
plt.yticks([])

# Añadir texto

f.text(0.65,0.7,'Mita Boundaries',color = 'black', size = 20,
        bbox=dict(facecolor='none', edgecolor='none', pad=5.0)) # bbox: añade una caja al texto, pad : largo y ancho de la caja

# Costumize legend

plt.legend(loc='upper left',
           title = "",frameon=True,
            bbox_to_anchor=(0, 0.15), prop={'size': 15})

<matplotlib.legend.Legend at 0x257a18bd720>

f, ax = plt.subplots(figsize=(15,15))

data_geo['geometry'].plot(color='none', edgecolor='gray', zorder=0.5, ax = ax)

data_geo['Point'].plot(color = '#DC143C', markersize=10, ax = ax, marker = "o", label="Distric's capital")

mita_boundary.plot(color = '#FFBF00' ,linewidth=2, ax = ax, label='Mita Boundary target') 

cx.add_basemap(ax, crs="EPSG:4326", source=cx.providers.NASAGIBS.ViirsEarthAtNight2012, attribution = False, zoom = 8)

plt.xticks([])
plt.yticks([])

# Añadir texto

f.text(0.65,0.7,'Mita Boundaries',color = 'black', size = 15,
        bbox=dict(facecolor='none', edgecolor='none', pad=5.0)) # bbox: añade una caja al texto, pad : largo y ancho de la caja

# Costumize legend

plt.legend(loc='upper left',
           title = "",frameon=True,
            bbox_to_anchor=(0, 0.15), prop={'size': 15})

<matplotlib.legend.Legend at 0x257c70c9e70>

4.0 Distances

data_geo

# selección de un point

data_geo.loc[0,"Point"]

# Plot de uno de los boundaries 

mita_boundary.loc[0,"geometry"] 

print( data_geo.loc[0,"Point"].distance(mita_boundary.loc[0,"geometry"] ) )

data_geo.loc[0,"Point"].distance(mita_boundary.loc[1,"geometry"] )

# Se buene observar que la primera capital está cercana al boundary de index : 0
# distance permite hallar la minima distancia geodésica
# Las unidades están en grados

0.3780424516996454

0.7991220422560418

# Coordenadas de Mina Potosi y Huancavelica

#potosi = (-19.571598207409757, -65.75495755044265)
#huanca = (-12.808993982792076, -74.97525207494975)

potosi = (-65.75495755044265, -19.571598207409757)
huanca = (-74.97525207494975, -12.808993982792076)

point = (data_geo.loc[0,"Point"].x, data_geo.loc[0,"Point"].y)

(-73.185, -13.95583333)

Referencia

Haversine

# Haversine distance asume que la tierra es una esfera. A esta medida tambien se llama Great-Circle

print("Distancia en kilometros: ",hs.haversine( huanca , point, unit = 'km'))

print("Distancia en metros: ", hs.haversine( huanca , point, unit = 'm'))

print("Distancia en millas: ", hs.haversine( huanca , point, unit = 'mi'))

Distancia en kilometros:  202.10753203071803
Distancia en metros:  202107.53203071802
Distancia en millas:  125.58379809010545

# Usando geopy 

print("Distancia geodésica en kilometros : ", geopy.distance.geodesic(huanca , point).km )

print("Distancia great-circle en kilometros : ", geopy.distance.great_circle(huanca , point).km )

Distancia geodésica en kilometros :  202.86294742689847
Distancia great-circle en kilometros :  202.1075383753187

Referencia de Geopy

Geo-py

# longitud 
print(mita_boundary.loc[1,"geometry"].length)

# area

data_geo.loc[0,"geometry"].area

# la unidad de medidas esta en grados (redianes)

6.399277433462688

0.011908123743330926

# Se declara el sistema geodésico de coordenadas: ellipsoide WGS84 ( EPSG 4326 )

geod = Geod(ellps="WGS84")

# minima distancia geodesica en metros 
# nearest_point halla los extremos de la linea mas corta entre boundary y point 
# LineString usas esos extremos para crear una linea
# luego se calcula la distancia de esa linea 

distance = geod.geometry_length(LineString(nearest_points(mita_boundary.loc[1,"geometry"], data_geo.loc[0,"Point"] )))

distance  # medida en metros 

87884.59488352262

5.0 Proyecciones a un plano cartesiano (Distancias Euclideanas)

data_geo.to_crs(epsg=32718, inplace = True) # proyeccion que permite calcular distancias precisas en el Perú 

data_geo

# Loop para proyectar cada columna del tipo de objeto "geometry"

for column in ['geometry','centroid','Point']:

    data_geo = data_geo.set_geometry(column).to_crs(epsg=32718) 
    
data_geo

# Proyeccion de la mita  
mita_boundary.to_crs(epsg=32718, inplace=True)
mita_boundary

print("Distancia minima geodésica en metros al Boundary : ", data_geo.loc[0,"Point"].distance(mita_boundary.loc[1,"geometry"] ) )

Distancia minima geodésica en metros al Boundary :  87900.81443711776

# longitud 

print(mita_boundary.loc[1,"geometry"].length/1000)  # distancia en kilometros del boundary

# area

data_geo.loc[0,"geometry"].area/1000**2  # 142.4 KM^2

698.6938653544977

142.422770084007

# Distancia euclideana 

distance = LineString(nearest_points(mita_boundary.loc[1,"geometry"], data_geo.loc[0,"Point"] )).length
distance 

# 87884.59488352262 cercano a la distancia geodésica , diferencia de 16 metros 

87900.81443711776

7.0 Buffers

# Buffer de linestring 

f, ax = plt.subplots(1, figsize=(8,8))

base1 = mita_boundary.loc[1,"geometry"].buffer(5000) # Linestring
gpd.GeoSeries(base1).plot(ax=ax) # Se convierte a geoseriies para graficarlo

base2 = mita_boundary.loc[0,"geometry"].buffer(3000) 
gpd.GeoSeries(base2).plot(ax=ax, color = 'red')

mita_boundary.plot(color='black', ax=ax)

plt.xticks([])
plt.yticks([])

([], [])

# Buffer de un polígono 

f, ax = plt.subplots(1, figsize=(8,8))

base = data_geo.loc[0,"geometry"].buffer(1000) # poligono del distrito con bordes extendido en 3 km
gpd.GeoSeries(base).plot(ax=ax) # Se convierte a geoseriies para graficarlo
gpd.GeoSeries(data_geo.loc[0,"geometry"]).plot(color="darkgray",ax=ax, linewidth= 1) # poligono del distrito original

plt.xticks([])
plt.yticks([])

([], [])

#gpd.tools.geocode("PUCP Peru, Avenida Universitaria, 15032, San Miguel, Lima")

References

Geopandas

Contexltly