GitHub Repository: robertopucp/1eco35_2022_2
Path: blob/main/Lab12/Geopandas.ipynb
²⁷¹⁴ views

Kernel: Python 3 (ipykernel)

Geopandas y operaciones ARCGIS

@author: Roberto mendoza

In [4]:

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 googlemaps #google request  

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 pyreadstat  # import spss files
import swifter # parallel procesing
import unidecode 

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

no module name, solo instalelo usando !pip install

1.0 Polygonos

In [5]:

# 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

Out[5]:

In [6]:

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

Out[6]:

<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

In [7]:

# 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

Out[7]:

In [8]:

dist_mita.crs

Out[8]:

<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

In [9]:

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

Out[9]:

In [10]:

mita_proy_cil.crs

Out[10]:

<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

In [11]:

# Se retorna al sistema de coordenadas iniciales

mita_cartesiano = dist_mita.to_crs(epsg=32718)
mita_cartesiano

Out[11]:

In [12]:

mita_cartesiano.plot()

Out[12]:

<AxesSubplot:>

2.0 Point

In [13]:

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

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

Out[13]:

In [14]:

centroids.crs

Out[14]:

<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

In [15]:

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

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

Out[15]:

In [16]:

centroids.plot()

Out[16]:

<AxesSubplot:>

In [17]:

# 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)

In [18]:

centroids

Out[18]:

3.0 Linestring

In [19]:

# Load a shapefile in Python 

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

In [20]:

mita_boundary

Out[20]:

In [21]:

mita_boundary.plot()

Out[21]:

<AxesSubplot:>

In [22]:

mita_boundary.crs

Out[22]:

<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

In [23]:

# load shapefile a nivel distrital

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

distritos

Out[23]:

In [24]:

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

distritos

Out[24]:

3.0 Merge Geospatial information

In [25]:

# select columns and then convert to geodataframe. 

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

In [26]:

centroids

Out[26]:

In [27]:

gdf

Out[27]:

In [28]:

gdf.info()

Out[28]:

<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

In [29]:

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

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

Out[29]:

3.2 Join - whitin

In [30]:

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

merge = gpd.sjoin( gdf, distritos , how="inner", op="within")
merge

Out[30]:

3.3 Join - Intersects

In [31]:

# 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()
data_geo

Out[31]:

In [32]:

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)

Out[32]:

<AxesSubplot:>

In [33]:

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)

Out[33]:

<AxesSubplot:>

In [34]:

mita_boundary.crs

Out[34]:

<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

In [35]:

# 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)

In [33]:

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 

# contextuly 
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})

Out[33]:

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [33], in <cell line: 16>()
mita_boundary.plot(color = '#834333' ,linewidth=4, ax = ax, label='Mita Boundary target') 
# añade un mapa de fondo 

# contextuly 
---> 16 cx.add_basemap(ax, crs="EPSG:4326", attribution = False, zoom = 11)
# sin valores en los ejes 
plt.xticks([])
File ~\anaconda3\envs\env_geopandas\lib\site-packages\contextily\plotting.py:121, in add_basemap(ax, zoom, source, interpolation, attribution, attribution_size, reset_extent, crs, resampling, **extra_imshow_args)
   left, right, bottom, top = _reproj_bb(
       left, right, bottom, top, crs, {"init": "epsg:3857"}
   )
# Download image
--> 121 image, extent = bounds2img(
   left, bottom, right, top, zoom=zoom, source=source, ll=False
)
# Warping
if crs is not None:
File ~\anaconda3\envs\env_geopandas\lib\site-packages\contextily\tile.py:222, in bounds2img(w, s, e, n, zoom, source, ll, wait, max_retries)
x, y, z = t.x, t.y, t.z
tile_url = provider.build_url(x=x, y=y, z=z)
--> 222 image = _fetch_tile(tile_url, wait, max_retries)
tiles.append(t)
arrays.append(image)
File ~\anaconda3\envs\env_geopandas\lib\site-packages\joblib\memory.py:594, in MemorizedFunc.__call__(self, *args, **kwargs)
def __call__(self, *args, **kwargs):
--> 594     return self._cached_call(args, kwargs)[0]
File ~\anaconda3\envs\env_geopandas\lib\site-packages\joblib\memory.py:537, in MemorizedFunc._cached_call(self, args, kwargs, shelving)
       must_call = True
if must_call:
--> 537     out, metadata = self.call(*args, **kwargs)
   if self.mmap_mode is not None:
       # Memmap the output at the first call to be consistent with
       # later calls
       if self._verbose:
File ~\anaconda3\envs\env_geopandas\lib\site-packages\joblib\memory.py:779, in MemorizedFunc.call(self, *args, **kwargs)
if self._verbose > 0:
   print(format_call(self.func, args, kwargs))
--> 779 output = self.func(*args, **kwargs)
self.store_backend.dump_item(
   [func_id, args_id], output, verbose=self._verbose)
duration = time.time() - start_time
File ~\anaconda3\envs\env_geopandas\lib\site-packages\contextily\tile.py:252, in _fetch_tile(tile_url, wait, max_retries)
@memory.cache
def _fetch_tile(tile_url, wait, max_retries):
--> 252     request = _retryer(tile_url, wait, max_retries)
   with io.BytesIO(request.content) as image_stream:
       image = Image.open(image_stream).convert("RGBA")
File ~\anaconda3\envs\env_geopandas\lib\site-packages\contextily\tile.py:395, in _retryer(tile_url, wait, max_retries)
"""
Retry a url many times in attempt to get a tile

   (...)
request object containing the web response.
"""
try:
--> 395     request = requests.get(tile_url, headers={"user-agent": USER_AGENT})
   request.raise_for_status()
except requests.HTTPError:
File ~\anaconda3\envs\env_geopandas\lib\site-packages\requests\api.py:75, in get(url, params, **kwargs)
def get(url, params=None, **kwargs):
   r"""Sends a GET request.

   :param url: URL for the new :class:`Request` object.
   (...)
   :rtype: requests.Response
   """
---> 75     return request('get', url, params=params, **kwargs)
File ~\anaconda3\envs\env_geopandas\lib\site-packages\requests\api.py:61, in request(method, url, **kwargs)
# By using the 'with' statement we are sure the session is closed, thus we
# avoid leaving sockets open which can trigger a ResourceWarning in some
# cases, and look like a memory leak in others.
with sessions.Session() as session:
---> 61     return session.request(method=method, url=url, **kwargs)
File ~\anaconda3\envs\env_geopandas\lib\site-packages\requests\sessions.py:529, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
send_kwargs = {
   'timeout': timeout,
   'allow_redirects': allow_redirects,
}
send_kwargs.update(settings)
--> 529 resp = self.send(prep, **send_kwargs)
return resp
File ~\anaconda3\envs\env_geopandas\lib\site-packages\requests\sessions.py:645, in Session.send(self, request, **kwargs)
start = preferred_clock()
# Send the request
--> 645 r = adapter.send(request, **kwargs)
# Total elapsed time of the request (approximately)
elapsed = preferred_clock() - start
File ~\anaconda3\envs\env_geopandas\lib\site-packages\requests\adapters.py:440, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
try:
   if not chunked:
--> 440         resp = conn.urlopen(
           method=request.method,
           url=url,
           body=request.body,
           headers=request.headers,
           redirect=False,
           assert_same_host=False,
           preload_content=False,
           decode_content=False,
           retries=self.max_retries,
           timeout=timeout
       )
   # Send the request.
   else:
       if hasattr(conn, 'proxy_pool'):
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\connectionpool.py:703, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, **response_kw)
   self._prepare_proxy(conn)
# Make the request on the httplib connection object.
--> 703 httplib_response = self._make_request(
   conn,
   method,
   url,
   timeout=timeout_obj,
   body=body,
   headers=headers,
   chunked=chunked,
)
# If we're going to release the connection in ``finally:``, then
# the response doesn't need to know about the connection. Otherwise
# it will also try to release it and we'll have a double-release
# mess.
response_conn = conn if not release_conn else None
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\connectionpool.py:386, in HTTPConnectionPool._make_request(self, conn, method, url, timeout, chunked, **httplib_request_kw)
# Trigger any extra validation we need to do.
try:
--> 386     self._validate_conn(conn)
except (SocketTimeout, BaseSSLError) as e:
   # Py2 raises this as a BaseSSLError, Py3 raises it as socket timeout.
   self._raise_timeout(err=e, url=url, timeout_value=conn.timeout)
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\connectionpool.py:1040, in HTTPSConnectionPool._validate_conn(self, conn)
# Force connect early to allow us to validate the connection.
if not getattr(conn, "sock", None):  # AppEngine might not have  `.sock`
-> 1040     conn.connect()
if not conn.is_verified:
   warnings.warn(
       (
           "Unverified HTTPS request is being made to host '%s'. "
   (...)
       InsecureRequestWarning,
   )
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\connection.py:358, in HTTPSConnection.connect(self)
def connect(self):
   # Add certificate verification
--> 358     self.sock = conn = self._new_conn()
   hostname = self.host
   tls_in_tls = False
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\connection.py:174, in HTTPConnection._new_conn(self)
   extra_kw["socket_options"] = self.socket_options
try:
--> 174     conn = connection.create_connection(
       (self._dns_host, self.port), self.timeout, **extra_kw
   )
except SocketTimeout:
   raise ConnectTimeoutError(
       self,
       "Connection to %s timed out. (connect timeout=%s)"
       % (self.host, self.timeout),
   )
File ~\anaconda3\envs\env_geopandas\lib\site-packages\urllib3\util\connection.py:85, in create_connection(address, timeout, source_address, socket_options)
   if source_address:
       sock.bind(source_address)
---> 85     sock.connect(sa)
   return sock
except socket.error as e:
KeyboardInterrupt: 

In [ ]:

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})

4.0 Distances

In [36]:

data_geo

Out[36]:

In [37]:

# selección de un point

data_geo.loc[0,"Point"]

Out[37]:

In [38]:

mita_boundary

Out[38]:

In [39]:

# Plot de uno de los boundaries 

mita_boundary.loc[0,"geometry"]

Out[39]:

In [40]:

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

Out[40]:

0.3780424516996454

0.7991220422560418

In [41]:

# 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)

In [42]:

point

Out[42]:

(-73.185, -13.95583333)

Referencia

Haversine

In [43]:

# 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'))

Out[43]:

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

In [44]:

# Usando geopy (from geopy import distance as )

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

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

Out[44]:

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

Referencia de Geopy

Geo-py

In [45]:

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

# area

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

# la unidad de medidas esta en grados (redianes)

Out[45]:

6.399277433462688

0.011908123743330926

In [46]:

# 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"] )))

In [47]:

distance  # medida en metros

Out[47]:

87884.59488352262

5.0 Proyecciones a un plano cartesiano (Distancias Euclideanas)

In [48]:

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

In [49]:

data_geo

Out[49]:

In [50]:

# 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

Out[50]:

In [51]:

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

Out[51]:

In [52]:

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

Out[52]:

Distancia minima euclideana en metros al Boundary :  87900.81102564473

In [53]:

# 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

Out[53]:

698.6938671485821

142.42277155133257

In [54]:

# 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

Out[54]:

87900.81102564473

7.0 Buffers

In [55]:

# Buffer de linestring 

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

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

base2 = mita_boundary.loc[0,"geometry"].buffer(3000)  # 3 kilometros a la redonda
gpd.GeoSeries(base2).plot(ax=ax, color = 'red')

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

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

Out[55]:

([], [])

In [79]:

# Buffer de un polígono 

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

base = data_geo.loc[0,"geometry"].buffer(2000) # poligono del distrito con bordes extendido en 2 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([])

Out[79]:

([], [])

8.0 Geocoding

1.0 Free geocoding and Google API

In [57]:

data2017, meta = pyreadstat.read_sav(r"../data/crime_data/2017/611-Modulo1312/sample_Denuncia_de_Delitos_2017.sav", encoding="latin1")
data2017

Out[57]:

In [58]:

# NOMBREDI : nombre del distrito
# IH205_DD: codigo de 2 digitos del departamento (callao : 07)
# H206: TIPO DE VIA DONDE OCURRIÓ EL DELITO
# IH207_B: NOMBRE DE LA DE VÍA DONDE OCURRIÓ EL DELITO
# IH207_A: Nº CUADRA DONDE OCURRIÓ EL DELITO

data2017[['NOMBREDI','H206','IH207_B','IH207_A']]

Out[58]:

In [59]:

data2017.IH207_A.unique()

Out[59]:

array(['SN', '1', '2', '6', 'NULL', '7', '8', '4', '13', '3', '21', '5',
       '11', '20', '10', '55', '9', '33', '14', '28', '49', '41', '38',
       '40', '31', '16', '15', '39', '3,', '2,', '5,', '1,', '33,', '6,',
       '31,', '4,', '20,', '8,', '18,', '16,', '40,', '7,', '9,', '19,',
       '30,'], dtype=object)

In [60]:

# filtro callaro y selecciono las primeras 20 observaciones 

callao = data2017.iloc[0:20]
callao[['NOMBREDI','H206','IH207_B','IH207_A']]

Out[60]:

In [61]:

#COMISARÃA BELLAVISTA : erro por no reconocer tilde. Se aplica unicode para corregirlo 

callao['IH207_B'] = callao['IH207_B'].apply( lambda x : unidecode.unidecode( x ) ) # unicode para corregir tildes no detectados
callao

Out[61]:

In [63]:

# Free geocoding 

coordinates = gpd.tools.geocode("PUCP Peru, Avenida Universitaria, 15032, San Miguel, Lima")
coordinates
# Point(x,y): x: longitud, y:latitud

Out[63]:

In [64]:

coordinates['address'][0]

Out[64]:

'PUCP Peru, Avenida Universitaria, 15032, Avenida Universitaria, San Miguel, Lima, Perú'

In [65]:

# Google API

gmaps = googlemaps.Client( key = 'AIzaSyA-Cj5npW0l2P7bZKP6WmQSq4uXSMvfHCM' ) 
# se crea el objeto gmaps con el token de permiso
    
result_api = gmaps.geocode( "PUCP, Avenida Universitaria, 15032, San Miguel, Lima" , region = 'PE' ) # region = PE 
result_api

# Google devuelve las solicitudes en formato Json

Out[65]:

[{'address_components': [{'long_name': '1801',
    'short_name': '1801',
    'types': ['street_number']},
   {'long_name': 'Avenida Universitaria',
    'short_name': 'Av. Universitaria',
    'types': ['route']},
   {'long_name': 'Fund Pando',
    'short_name': 'Fund Pando',
    'types': ['political', 'sublocality', 'sublocality_level_1']},
   {'long_name': 'San Miguel',
    'short_name': 'San Miguel',
    'types': ['locality', 'political']},
   {'long_name': 'Provincia de Lima',
    'short_name': 'Provincia de Lima',
    'types': ['administrative_area_level_2', 'political']},
   {'long_name': 'Provincia de Lima',
    'short_name': 'Provincia de Lima',
    'types': ['administrative_area_level_1', 'political']},
   {'long_name': 'Peru',
    'short_name': 'PE',
    'types': ['country', 'political']},
   {'long_name': '15088', 'short_name': '15088', 'types': ['postal_code']}],
  'formatted_address': 'Av. Universitaria 1801, San Miguel 15088, Peru',
  'geometry': {'location': {'lat': -12.0689502, 'lng': -77.0780608},
   'location_type': 'ROOFTOP',
   'viewport': {'northeast': {'lat': -12.0676218697085,
     'lng': -77.07669776970849},
    'southwest': {'lat': -12.0703198302915, 'lng': -77.0793957302915}}},
  'partial_match': True,
  'place_id': 'ChIJpUAI1BLJBZERLobll7e_oNc',
  'plus_code': {'compound_code': 'WWJC+CQ San Miguel',
   'global_code': '57V4WWJC+CQ'},
  'types': ['establishment', 'point_of_interest', 'university']}]

In [66]:

result_api[0]

Out[66]:

{'address_components': [{'long_name': '1801',
   'short_name': '1801',
   'types': ['street_number']},
  {'long_name': 'Avenida Universitaria',
   'short_name': 'Av. Universitaria',
   'types': ['route']},
  {'long_name': 'Fund Pando',
   'short_name': 'Fund Pando',
   'types': ['political', 'sublocality', 'sublocality_level_1']},
  {'long_name': 'San Miguel',
   'short_name': 'San Miguel',
   'types': ['locality', 'political']},
  {'long_name': 'Provincia de Lima',
   'short_name': 'Provincia de Lima',
   'types': ['administrative_area_level_2', 'political']},
  {'long_name': 'Provincia de Lima',
   'short_name': 'Provincia de Lima',
   'types': ['administrative_area_level_1', 'political']},
  {'long_name': 'Peru', 'short_name': 'PE', 'types': ['country', 'political']},
  {'long_name': '15088', 'short_name': '15088', 'types': ['postal_code']}],
 'formatted_address': 'Av. Universitaria 1801, San Miguel 15088, Peru',
 'geometry': {'location': {'lat': -12.0689502, 'lng': -77.0780608},
  'location_type': 'ROOFTOP',
  'viewport': {'northeast': {'lat': -12.0676218697085,
    'lng': -77.07669776970849},
   'southwest': {'lat': -12.0703198302915, 'lng': -77.0793957302915}}},
 'partial_match': True,
 'place_id': 'ChIJpUAI1BLJBZERLobll7e_oNc',
 'plus_code': {'compound_code': 'WWJC+CQ San Miguel',
  'global_code': '57V4WWJC+CQ'},
 'types': ['establishment', 'point_of_interest', 'university']}

In [67]:

result_api[0].keys()

Out[67]:

dict_keys(['address_components', 'formatted_address', 'geometry', 'partial_match', 'place_id', 'plus_code', 'types'])

In [68]:

result_api[0]['geometry']['location']

Out[68]:

{'lat': -12.0689502, 'lng': -77.0780608}

In [69]:

lat = result_api[0]['geometry']['location']['lat']
lon = result_api[0]['geometry']['location']['lng']  
(lon,lat)

Out[69]:

(-77.0780608, -12.0689502)

In [70]:

result_api[0]['types']

Out[70]:

['establishment', 'point_of_interest', 'university']

In [71]:

# componentes de la dirección 

result_api[0]['address_components'][4]

Out[71]:

{'long_name': 'Provincia de Lima',
 'short_name': 'Provincia de Lima',
 'types': ['administrative_area_level_2', 'political']}

In [72]:

callao[['NOMBREDI','H206','IH207_B','IH207_A']]

Out[72]:

In [73]:

# crear el tipo de calle 

# condiciones

filt1 = callao['H206'] == 1
filt2 = callao['H206'] == 2
filt3 = callao['H206'] == 3
filt4 = callao['H206'] == 4
filt5 = callao['H206'] == 5
filt6 = callao['H206'] == 6
filt7 = callao['H206'] == 7

# se filtra, se crea la vafriable tipo_calle y se asigno valores 
callao.loc[filt1, 'tipo_calle'] = "Avenida"
callao.loc[filt2, 'tipo_calle'] = "Jiron"
callao.loc[filt3, 'tipo_calle'] = "Calle"
callao.loc[filt4, 'tipo_calle'] = "Pasaje"
callao.loc[filt5, 'tipo_calle'] = "Carretera"
callao.loc[filt6, 'tipo_calle'] = "Prolongación"
callao.loc[filt7, 'tipo_calle'] = " "

callao.replace({'IH207_B': "SN", 'IH207_A': "SN"}, " ", inplace = True)

In [74]:

callao[['NOMBREDI','H206','tipo_calle','IH207_B','IH207_A']]

Out[74]:

In [75]:

# function to geocode

def georef_geocode( row ):
       
    address = f"{row['tipo_calle']} {row['IH207_B']} {row['IH207_A']}, {row['NOMBREDI']}, Callao, Perú"
    
    geo = gpd.tools.geocode(address)
    
    # Information
    
    try:
        lat = geo['geometry'].y[0]
        lon = geo['geometry'].x[0]  
        
    except Exception as e:
        
        lat = np.nan
        lon = np.nan
    
    return ( lat, lon )

# function to Geocoding Google API


def georef_google_api( row ):
       
    address = f"{row['tipo_calle']} {row['IH207_B']} {row['IH207_A']}, {row['NOMBREDI']}, Callao, Perú"
    
    # Set Geolocation
    
    gmaps = googlemaps.Client( key = 'AIzaSyCJoBnbPDYdUMLngeRKsLLMAsUFdFF28do' )
    
    result_api = gmaps.geocode( address , region = 'PE' )
    
    # Information
    try:
        lat = result_api[0]['geometry']['location']['lat']
        lon = result_api[0]['geometry']['location']['lng']   
        
    except Exception :
        
        
        lat = np.nan
        lon = np.nan
    
    return ( lat, lon )

In [76]:

callao['coordinates_geo'] = callao.swifter.apply( lambda x: georef_geocode( x )  , axis = 1 )  # apply function by each row
callao['coordinates_google'] = callao.swifter.apply( lambda x: georef_google_api( x )  , axis = 1 )  # apply function by each row

Out[76]:

Pandas Apply:   0%|          | 0/20 [00:00<?, ?it/s]

Pandas Apply:   0%|          | 0/20 [00:00<?, ?it/s]

In [77]:

callao[['lat_geo', 'lng_geo']] = pd.DataFrame(callao.coordinates_geo.tolist() ,  index = callao.index )
callao[['lat_google', 'lng_google']] = pd.DataFrame(callao.coordinates_google.tolist() ,  index = callao.index )

In [78]:

callao[['tipo_calle','IH207_B','IH207_A','lat_geo', 'lng_geo','lat_google', 'lng_google']]

Out[78]:

References

Geopandas

Contexltly

Geopandas y operaciones ARCGIS

@author: Roberto mendoza

no module name, solo instalelo usando !pip install

1.0 Polygonos

Sistema de coordendas proyectada usada en la Mita

2.0 Point

3.0 Linestring

3.0 Merge Geospatial information

3.1 Join - Contains

3.2 Join - whitin

3.3 Join - Intersects

4.0 Distances

Referencia

Referencia de Geopy

5.0 Proyecciones a un plano cartesiano (Distancias Euclideanas)

7.0 Buffers

8.0 Geocoding

1.0 Free geocoding and Google API

References

Product

Resources

Company