CoCalc -- Geopandas

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

Kernel: Python 3 (ipykernel)

¿Geopandas?

Geopandas enables the use and manipulation of geospatial data in Python. It extends the common datatype used in pandas to allow for the many and unique geometric operations: GeoSeries and GeoDataFrame. Geopandas is also built on top of shapely for its geometric operation; its underlying datatype allows Geopandas to run blazingly fast and is appropriate for many machine learning pipelines that require large geospatial datasets.

Instalation

Following steps:

1. create new environment

conda creare -name new_env

2. Activate new environment

conda activate new_env

3. Install package

conda install -c conda-forge geopandas

Geopandas data types

Shapefile (.shp): vector data that stores geospatial information such as topology, shape geometry, coordinate etc

GeoJSON (.geojson): geospatial information and relevant characteristics (object or string)

Finally, Geopandas creates an additional geometry column

Coordinate Reference Systems (crs)

The coordinate reference systems tell us how GeodataFrame o GeoSeries relate to places on Earth

relevant operations: setting a projection and re-projecting

setting projection: geopandas.read_file() most of time includes projection information. Check using "Geoseries.crs" If GeoDataframe does not include projection, hence use data["geometry"].set_crs(epsg=4326) for WGS-84.

WGS-84:

The Global Positioning System uses the World Geodetic System (WGS84) as its reference coordinate system. It's made up of a reference ellipsoid, a standard coordinate system, altitude data, and a geoid

Reprojecting:

In somes cases, due to distortions in bidimensional plane and your application of geospatial data, re-projection is required.

A commonly used default is Universal Transverse Mercator (UTM)

In [1]:

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

In [34]:

# objects from Geopandas 

import geopandas as gpd
from geopandas import GeoSeries
from shapely.geometry import Point, LineString

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

In [4]:

# chardet permite reconocer el formato del texto dentro del archiv .csv

base = open(r'../data/geopandas_data/Datos_panel_distritos.csv', 'rb').read()
det = chardet.detect(base)
charenc = det['encoding']
charenc

Out[4]:

'utf-8'

Covid Data

In [6]:

# Import csv file from panel covid information at distric level
# Panel data from march-2020 to first two wees of january-2022

cv_data = pd.read_csv(r'../data/geopandas_data/Datos_panel_distritos.csv', encoding = charenc)
cv_data.head( 5 )

Out[6]:

In [11]:

cv_data.dtypes

Out[11]:

Código              int64
Departamento       object
Provincia          object
Distrito           object
Mes                 int64
Año                 int64
Casos               int64
Muertes Covid       int64
Muertes totales     int64
dtype: object

In [7]:

# rename variable to merge later 
cv_data.rename({'Código':'UBIGEO'}, axis =1 , inplace = True)
cv_data.head( 5 )

Out[7]:

Public Health Centers Data

In [9]:

# import csv file of public health centers: MINSA, ESSALUD and Regional Goverment
# coordinates information (NORTE: longitud) (ESTE: latitud)

health_centers = pd.read_csv(r'../data/geopandas_data/centros_hospitalarios.csv')
health_centers.head( 5 )

# Latitud (Y) ;  Longitud (X)

Out[9]:

In [11]:

# We can see longitude cordinates

health_centers.NORTE

Out[11]:

    -81.310632
    -81.306508
    -81.303974
    -81.275139
    -81.273840
          ...    
        NaN
        NaN
        NaN
        NaN
        NaN
Name: NORTE, Length: 8464, dtype: float64

In [12]:

# DataFrame to GeoDataFrame

Public_health = gpd.GeoDataFrame( health_centers, geometry=gpd.points_from_xy( health_centers.NORTE, health_centers.ESTE ) )
Public_health.head( 5 )

Out[12]:

Types and classification public centers

We focus on Postas médicas and Hospitales

code2 (Type)

ESTABLECIMIENTO DE SALUD SIN INTERNAMIENTO: 1

SERVICIO MÉDICO DE APOYO: 2

ESTABLECIMIENTO DE SALUD CON INTERNAMIENTO: 3

code1 (Clasification)

ATENCION PRE HOSPITALARIA: 1

CENTROS DE ATENCION PARA DEPENDIENTES A SUSTANCIAS PSICOACTIVAS Y OTRAS DEPENDENCIAS: 2

CENTROS DE SALUD CON CAMAS DE INTERNAMIENTO,CENTROS DE SALUD CON CAMAS DE INTERNAMIENTO: 3

CENTROS DE SALUD O CENTROS MEDICOS: 4

CENTROS DE SALUD O CENTROS MEDICOS,CENTROS DE SALUD O CENTROS MEDICOS: 5

CENTROS DE VACUNACION: 6

CENTROS MEDICOS ESPECIALIZADOS: 7

CONSULTORIOS MEDICOS Y DE OTROS PROFESIONALES DE LA SALUD: 8

ESTABLECIMIENTOS DE RECUPERACION O REPOSO: 9

HEMODIALISIS: 10

HOSPITALES O CLINICAS DE ATENCION ESPECIALIZADA: 11

HOSPITALES O CLINICAS DE ATENCION GENERAL: 12

INSTITUTOS DE SALUD ESPECIALIZADOS: 13

MEDICINA FISICA,REHABILITACION: 14

PATOLOGIA CLINICA: 15

PATOLOGIA CLINICA,ANATOMIA PATOLOGICA: 16

PATOLOGIA CLINICA,DIAGNOSTICO POR IMAGENES: 17

POLICLINICOS: 18

PUESTOS DE SALUD O POSTAS DE SALUD: 19

PUESTOS DE SALUD O POSTAS DE SALUD,PUESTOS DE SALUD O POSTAS DE SALUD: 19

SERVICIO DE TRASLADO DE PACIENTES: 20

SERVICIO DE TRASLADO DE PACIENTES,ATENCION DOMICILIARIA: 21

SERVICIO DE TRASLADO DE PACIENTES,ATENCION DOMICILIARIA,ATENCION PRE HOSPITALARIA: 22

SERVICIO DE TRASLADO DE PACIENTES,ATENCION PRE HOSPITALARIA: 23

CENTROS DE SALUD CON CAMAS DE INTERNAMIENTO: 24

In [13]:

# code1 = 19 to postas médicas 
# code1 = 12 & 24 to Hospitals

Public_health1 = Public_health[Public_health.code1 == "19"]

Public_health2 = Public_health[(Public_health.code1 == "12") | (Public_health.code1 == "24")]
Public_health2.head( 5 )

Out[13]:

Shape File

In [14]:

# Upload shape file at district level

maps = gpd.read_file(r'../data/geopandas_data/shape_file/DISTRITOS.shp')

In [15]:

maps.head(10)

Out[15]:

In [16]:

# Polygon by each district 

maps.geometry

Out[16]:

     POLYGON ((-75.31797 -9.29529, -75.31710 -9.297...
     POLYGON ((-74.64136 -8.82302, -74.64036 -8.828...
     POLYGON ((-75.02253 -8.74193, -75.02267 -8.742...
     POLYGON ((-75.13864 -8.56712, -75.13956 -8.569...
     POLYGON ((-75.01589 -8.44637, -75.01585 -8.446...
                              ...                        
  POLYGON ((-76.08083 -9.13017, -76.08026 -9.130...
  POLYGON ((-75.88828 -9.00906, -75.88756 -9.010...
  POLYGON ((-75.91141 -8.88593, -75.91182 -8.886...
  POLYGON ((-70.13203 -15.33382, -70.12355 -15.3...
  POLYGON ((-77.37630 -10.89734, -77.37604 -10.8...
Name: geometry, Length: 1873, dtype: geometry

In [17]:

maps = maps[['IDDIST', 'geometry']]
maps.rename({'IDDIST':'UBIGEO'}, axis =1, inplace = True )
maps

Out[17]:

In [18]:

# UBIGEO from srting to int

maps['UBIGEO'] = maps['UBIGEO'].astype(str).astype(int)
maps

Out[18]:

In [20]:

# Plot at district level 

fig, ax = plt.subplots(figsize=(20, 20))

maps.plot(  ax=ax, color='blue',
            edgecolor='black' )
plt.show()

Out[20]:

Merge: Covid Data & Shapefile

In [24]:

# Merge using inner way to drop missing values 

dataset_cv = pd.merge(maps, cv_data, how="left", on="UBIGEO")
dataset_cv

# cv_data información de casos de covid por distrito y (años y mes)

Out[24]:

In [25]:

# Function to filter dataset_covid by year, month and administrative space(province)

def subdata( a , b, c ):
    
    if c == "All-Peru":
        
        resultado = dataset_cv.loc[(dataset_cv['Año']== a) & (dataset_cv['Mes']== b)]
    else:
        resultado = dataset_cv.loc[(dataset_cv['Año']== a) & (dataset_cv['Mes']== b) & (dataset_cv['Provincia']== c)]
    return resultado

In [30]:

#Second wave 
base = subdata( 2021 , 8, "All-Peru")  # all peru Aug-2020

#Second wave
base_lima = subdata( 2021 , 8, "LIMA")  # Lima Aug-2020

Plot covid infection at disctrict level

In [31]:

# Plot covid infection at disctrict level 

base.plot( column='Casos', cmap='Reds', 
          figsize=(20, 20),
          linestyle='-',
          edgecolor='gray',
           missing_kwds= dict(color = "#DADADB"), 
          legend = True )

Out[31]:

<AxesSubplot:>

In [32]:

# Covid & Lima Matropolitana 

base_lima.plot(column='Casos', cmap='Purples', figsize=(20, 20),
                      edgecolor='black', 
                      legend = True)

# save map
plt.savefig(r'Lima_covid_Heatmap.png')

Out[32]:

In [33]:

# Observamos el sitema de coordenadas (EPSG:4326)

base_lima.crs

Out[33]:

<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

Get Centroid

In [35]:

# Centroid 

base_lima['centroid'] = base_lima['geometry'].centroid

In [37]:

base_lima.head(10)

Out[37]:

In [38]:

props = dict(boxstyle='round', facecolor = 'linen', alpha =1)
props

Out[38]:

{'boxstyle': 'round', 'facecolor': 'linen', 'alpha': 1}

In [39]:

# check distric's name 
for point in base_lima.iterrows():
    print(point[1]['Distrito'])

Out[39]:

LURIN
ANCON
SAN ISIDRO
MAGDALENA DEL MAR
SAN BORJA
LINCE
SANTIAGO DE SURCO
PUEBLO LIBRE
SANTA ROSA
CARABAYLLO
JESUS MARIA
PUCUSANA
SANTA MARIA DEL MAR
CHORRILLOS
BARRANCO
VILLA MARIA DEL TRIUNFO
SAN JUAN DE MIRAFLORES
MIRAFLORES
SURQUILLO
PUNTA NEGRA
SAN LUIS
SAN MIGUEL
VILLA EL SALVADOR
PUNTA HERMOSA
SAN BARTOLO
PACHACAMAC
LA VICTORIA
LA MOLINA
BRE<d1>A
SANTA ANITA
CIENEGUILLA
LIMA
EL AGUSTINO
RIMAC
ATE
CHACLACAYO
INDEPENDENCIA
SAN MARTIN DE PORRES
LOS OLIVOS
COMAS
LURIGANCHO
SAN JUAN DE LURIGANCHO
PUENTE PIEDRA

Add District Names in Centroid

In [41]:

fig, ax = plt.subplots(figsize=(20, 20))

# Heatmap Lima and Covid infection

base_lima.plot(column='Casos', cmap='viridis', linestyle='--',
                      edgecolor='black', 
                      legend = True, ax=ax)

Out[41]:

<AxesSubplot:>

In [42]:

# Add text by district's name 

fig, ax = plt.subplots(figsize=(20, 20))

# Heatmap Lima and Covid infection

base_lima.plot(column='Casos', cmap='viridis', linestyle='--',
                      edgecolor='black', 
                      legend = True, ax=ax)

# add district name in litle box

props = dict(boxstyle='round', facecolor = 'linen', alpha =1)

for point in base_lima.iterrows():
        ax.text( point[1]['centroid'].x,
                 point[1]['centroid'].y,
                 point[1]['Distrito'],
                 horizontalalignment='center',
                 fontsize=10,
                 bbox=props )

#plt.savefig(r'Lima_covid_Heatmap_2.png')

Out[42]:

In [45]:

#Using legend by district 


# Customize legend 
leg_kwds={'title':'District Name',
          'loc': 'upper left',
          'bbox_to_anchor':(1, 1.03),
          'ncol':3}


base_lima.plot(column='Distrito', cmap='tab20c', figsize=(20, 20), legend = True, linestyle='-', 
               edgecolor='black',
               legend_kwds = leg_kwds)

Out[45]:

<AxesSubplot:>

In [46]:

# Postas Médicas and Hospitales públicos - LIMA

PH_lima_postas = Public_health1[Public_health1.Provincia == "LIMA"]

PH_lima_hospitales = Public_health2[Public_health2.Provincia == "LIMA"]
PH_lima_hospitales

Out[46]:

Covid y Postas

In [47]:

PH_lima_postas

Out[47]:

In [51]:

# Covid deaths and Postas médicas 

fig, ax = plt.subplots(figsize=(20, 20))

#Heatmap Lima and Covid infection

base_lima.plot( column='Muertes Covid', 
               cmap='viridis', 
               linestyle='--',
               edgecolor='black', 
               legend = True, 
               ax=ax )

PH_lima_postas.plot(ax = ax, color = 'red', markersize=50, edgecolor='black')

Out[51]:

<AxesSubplot:>

Hospitales Públicos y Postas Médicas

In [52]:

# Plot de Postas anh Hospital 

ax = base_lima.plot(figsize=(25, 25), color='none', edgecolor='black', zorder=3)

PH_lima_postas.plot(color='teal', label = "Postas médicas", markersize=25, ax=ax)
PH_lima_hospitales.plot(color='red', label = "Hospitales públicos", markersize=75, ax=ax)

plt.legend(loc='upper right')

Out[52]:

<matplotlib.legend.Legend at 0x2109022df60>

In [54]:

# Covid infection and Postas médicas 

fig, ax = plt.subplots(figsize=(20, 20))

#Heatmap Lima and Covid infection

base_lima.plot(column='Casos', cmap='viridis', linestyle='--',
                      edgecolor='black', 
                      legend = True, ax=ax)

PH_lima_postas.plot(color='red', label = "Postas médicas", markersize=25, ax=ax)
PH_lima_hospitales.plot(color='blue', label = "Hospitales públicos", markersize=75, ax=ax)

plt.legend(loc='upper right')

Out[54]:

<matplotlib.legend.Legend at 0x210906aceb0>

Relative Size of Markers from panda series

In [56]:

# drop raws with nan values 

PH_lima_hospitales = PH_lima_hospitales.dropna(subset=['CAMAS'])

In [58]:

# Covid infection and Public Hospital by number of beds 

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

markersize=PH_lima_hospitales['CAMAS']

#Heatmap Lima and Covid infection

base_lima.plot(column='Casos', cmap='viridis', linestyle='--',
                      edgecolor='black', 
                      legend = True, ax=ax)

PH_lima_hospitales.plot( ax=ax, color='red', label = "Hospitales públicos", markersize = markersize,
                       edgecolor='black')
plt.legend(loc='upper right')

Out[58]:

<matplotlib.legend.Legend at 0x21095ae0f40>

In [59]:

# Covid infection & postas médicas at the national level 

gf = base.plot(column='Casos', cmap='Oranges', 
          figsize=(30, 30),
          linestyle='-',
          edgecolor='black',
          legend = True)

Public_health1.plot(ax = gf, color = 'red', markersize=15)

Out[59]:

<AxesSubplot:>

Select Districts

In [60]:

# districs name as a index
base_lima = base_lima.set_index("Distrito")

In [61]:

base_lima.head(20)

Out[61]:

In [64]:

#Plot selectec districs 

ATE = base_lima.loc["ATE", "geometry"]
SJL = base_lima.loc["SAN JUAN DE LURIGANCHO", "geometry"]
SI = base_lima.loc["SAN ISIDRO", "geometry"]
LC = base_lima.loc["LINCE", "geometry"]

gpd.GeoSeries([ATE,SJL,SI,LC]).plot(figsize=(8, 8),linestyle='-',
                                    edgecolor='black')

Out[64]:

<AxesSubplot:>

Limits - Boarders

In [65]:

# Logic operations:

# Lince and San Isidro share limits 

LC.touches(SI)

Out[65]:

True

In [66]:

LC.touches(ATE)

Out[66]:

False

In [67]:

#Distance ATE y (SAN ISIDRO, LINCE)

ATE.distance(SI)

Out[67]:

0.021430816237973403

In [68]:

ATE.distance(LC)

Out[68]:

0.03353670040995021

In [69]:

# Add line that cross Lince and ATE

ate_point = base_lima['centroid']['ATE']
lc_point = base_lima['centroid']['LINCE']
sb_point = base_lima['centroid']['SAN BORJA']
lr_point = base_lima['centroid']['LURIN']
ac_point = base_lima['centroid']['ANCON']

In [70]:

# Add line that cross Lince and ATE

ate_lince = LineString(zip((ate_point.x , lc_point.x ), (ate_point.y, lc_point.y)))

fig, ax = plt.subplots(figsize=(20,20))
base_lima.plot(color='lightgray', ax=ax, linestyle='-',
                                    edgecolor='black')

gpd.GeoSeries([ate_lince]).plot(color='red', linestyle='--', ax = ax)

Out[70]:

<AxesSubplot:>

In [71]:

base_lima['area']=base_lima.area
base_lima

Out[71]:

In [72]:

#Distamce between San Juan de Lurigancho to other districs 

SJL_point = base_lima['centroid']["SAN JUAN DE LURIGANCHO"]
    
base_lima['distance'] = base_lima['centroid'].distance(SJL_point)

In [73]:

base_lima['distance'].sort_values(ascending=False)

Out[73]:

Distrito
PUCUSANA                   0.555547
SANTA MARIA DEL MAR        0.505632
SAN BARTOLO                0.490367
PUNTA NEGRA                0.443395
PUNTA HERMOSA              0.395445
LURIN                      0.334878
VILLA EL SALVADOR          0.274567
ANCON                      0.273336
PACHACAMAC                 0.268347
CHORRILLOS                 0.249065
SANTA ROSA                 0.237993
CIENEGUILLA                0.233781
VILLA MARIA DEL TRIUNFO    0.230486
SAN JUAN DE MIRAFLORES     0.211952
CHACLACAYO                 0.209580
BARRANCO                   0.204123
MIRAFLORES                 0.184133
SANTIAGO DE SURCO          0.181479
SAN MIGUEL                 0.176307
MAGDALENA DEL MAR          0.175867
SURQUILLO                  0.172533
LURIGANCHO                 0.170624
SAN ISIDRO                 0.166317
PUEBLO LIBRE               0.159857
LINCE                      0.153848
SAN BORJA                  0.153424
JESUS MARIA                0.152794
LA MOLINA                  0.149307
CARABAYLLO                 0.138590
BRE<d1>A                   0.138095
PUENTE PIEDRA              0.136796
LA VICTORIA                0.134441
SAN LUIS                   0.130681
LIMA                       0.130485
ATE                        0.127825
SAN MARTIN DE PORRES       0.125958
LOS OLIVOS                 0.106220
RIMAC                      0.098832
SANTA ANITA                0.097599
EL AGUSTINO                0.095466
INDEPENDENCIA              0.085326
COMAS                      0.069492
SAN JUAN DE LURIGANCHO     0.000000
Name: distance, dtype: float64

In [74]:


base_lima['distance'].sort_values(ascending=False).plot.bar(figsize=(12, 7), color = 'red', edgecolor='black')

Out[74]:

<AxesSubplot:xlabel='Distrito'>

In [75]:

# SJL and other districs 

d1 = LineString(zip((SJL_point.x , lc_point.x ), (SJL_point.y, lc_point.y)))
d2 = LineString(zip((SJL_point.x , ate_point.x ), (SJL_point.y, ate_point.y)))
d3 = LineString(zip((SJL_point.x , sb_point.x ), (SJL_point.y, sb_point.y)))
d4 = LineString(zip((SJL_point.x , lr_point.x ), (SJL_point.y, lr_point.y)))
d5 = LineString(zip((SJL_point.x , ac_point.x ), (SJL_point.y, ac_point.y)))

fig, ax = plt.subplots(figsize=(20,20))
base_lima.plot(color='white', ax=ax, linestyle='-',
                                    edgecolor='black')

for i in range(1,6):
    gpd.GeoSeries(globals()[f'd{i}']).plot(color='red', linestyle='--', ax = ax, linewidth = 2.5)

Out[75]:

ArcGis Operations in Geopandas

In [77]:

dist_shp = gpd.read_file(r'../data/geopandas_data/shape_file/DISTRITOS.shp')

In [78]:

dist_shp

Out[78]:

Get Centroids

Get the centroids of each geometry (district).

In [79]:

dist_shp.centroid

Out[79]:

      POINT (-75.46654 -9.63436)
      POINT (-74.86638 -9.00132)
      POINT (-75.04802 -8.82493)
      POINT (-75.27992 -8.95052)
      POINT (-75.04903 -8.64130)
                   ...             
   POINT (-76.03798 -9.20808)
   POINT (-76.00137 -9.08220)
   POINT (-75.98297 -8.98490)
  POINT (-70.12642 -15.40623)
  POINT (-77.44523 -11.00358)
Length: 1873, dtype: geometry

In [80]:

dist_shp.crs

Out[80]:

<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

Proyecting from epsg = 4326 to bidimensional epsg = 24891

We get warnings since we need to reproject this shapefile object to a CRS in meters. This CRS should be a plane area to the specific zone of the shapefile (more info). In the case of Peru, we use 24891 (more info).

In [81]:

projected_dist_shp = dist_shp.to_crs( epsg = 24891 ) # proyectar de un sistema a otro

projected_dist_shp[ 'centroids'] = projected_dist_shp.centroid

In [82]:

projected_dist_shp

Out[82]:

Then, we reproject again the centroids to the original CRS.

In [83]:

dist_shp[ 'Centroid_Latitude' ] = projected_dist_shp[ 'centroids' ].to_crs( epsg = 4326 ).y
dist_shp[ 'Centroid_Longitude' ] = projected_dist_shp[ 'centroids' ].to_crs( epsg = 4326 ).x

Intersection

In [84]:

health_centers = pd.read_csv( r"../data/geopandas_data/salud_places_peru.csv")

In [85]:

ht1 = health_centers[ health_centers.latitud != 0 ]

In [86]:

# Make ht1 a geopandas
ht1_geo = gpd.GeoDataFrame( ht1 , 
                           geometry= gpd.points_from_xy( 
                                      ht1.longitud, 
                                      ht1.latitud))

Check CRS

We found that geolocation is recognized by google maps. We set the CRS that correspond to google maps.

In [87]:

ht1_geo_crs = ht1_geo.set_crs( epsg = 4326 )

Find Ubigeo for each point

In [88]:

# Just take 10% of observations
ht2_geo = ht1_geo_crs.sample( frac = 0.10 )

In [89]:

# Just take 10% of observations
ht2_geo

Out[89]:

Aggregation

We have shapefiles at district levels. It is possible to aggregate districts at Department or Province level.

In [90]:

fig, ax = plt.subplots( figsize = ( 12, 15 ) )
dist_shp.plot( ax = ax )

Out[90]:

<AxesSubplot:>

In [91]:

prov_shp = dist_shp.dissolve( by = 'IDPROV' )

In [92]:

fig, ax = plt.subplots( figsize = ( 12, 15 ) )
prov_shp .plot( ax = ax )

Out[92]:

<AxesSubplot:>

In [93]:

dpt_shp = dist_shp.dissolve( by = 'IDDPTO' )

In [94]:

fig, ax = plt.subplots( figsize = ( 12, 15 ) )
dpt_shp.plot( ax = ax )

Out[94]:

<AxesSubplot:>

Checking CRS

In [96]:

# Border of Switzerland
ch_borders = gpd.read_file( r'../data/geopandas_data/eu_maps/CHE_adm0.shp')

# Border of Germany
de_borders = gpd.read_file( r'../data/geopandas_data/eu_maps/gadm36_DEU_0.shp')

# Switzerland lakes shape files
ch_lakes = gpd.GeoDataFrame.from_file( r'../data/geopandas_data/eu_maps/g2s15.shp')

# muni
Geo_muni = gpd.read_file( r"../data/geopandas_data/eu_maps/VG250_GEM.shp" )

In [97]:

ch_borders

Out[97]:

Plots with different CRS

In [98]:


fig = plt.figure(figsize=(40, 20), dpi=100)
ax = fig.add_subplot(111, frame_on = False )

# German Plot
Geo_muni.plot( ax = ax )

# switzerland borders
ch_borders.geometry.boundary.plot(color = 'black', ax=ax, edgecolor='0.2', linewidth=1.5) 

# # Swiss lakes
ch_lakes.plot(color='lightblue' , ax=ax, edgecolor='0.2', linewidth=0.4, legend=True)

ax.axis('off')

Out[98]:

(-46058.36469219931, 967356.691928138, -305024.0777151486, 6406513.607306512)

Check CRS

In [99]:

shp_files = [ ch_borders, de_borders, ch_lakes, Geo_muni ]

In [100]:

for shp in shp_files:
    print( shp.crs )

Out[100]:

epsg:4326
epsg:4326
epsg:21781
epsg:25832

Reproject CRS

In [101]:

Geo_muni_crs = Geo_muni.to_crs( 21781 )

In [102]:

ch_borders_crs = ch_borders.to_crs( 21781 )

In [103]:


fig = plt.figure(figsize=(40, 20), dpi=100)
ax = fig.add_subplot(111, frame_on = False )

# German Plot
Geo_muni_crs.plot( ax = ax )

# switzerland borders
ch_borders_crs.geometry.boundary.plot(color = 'black', ax=ax, edgecolor='0.2', linewidth=1.5) 

# # Swiss lakes
ch_lakes.plot(color='lightblue' , ax=ax, edgecolor='0.2', linewidth=0.4, legend=True)

ax.axis('off')

Out[103]:

(453126.0049465788, 1163640.4776774019, 11798.918320549732, 1157415.473468453)

References:

¿Geopandas?

Instalation

1. create new environment

2. Activate new environment

3. Install package

Geopandas data types

Coordinate Reference Systems (crs)

relevant operations: setting a projection and re-projecting

WGS-84:

Reprojecting:

Covid Data

Public Health Centers Data

Types and classification public centers

We focus on Postas médicas and Hospitales

code2 (Type)

code1 (Clasification)

Shape File

Merge: Covid Data & Shapefile

Plot covid infection at disctrict level

Get Centroid

Add District Names in Centroid

Covid y Postas

Hospitales Públicos y Postas Médicas

Relative Size of Markers from panda series

Select Districts

Limits - Boarders

ArcGis Operations in Geopandas

Get Centroids

Proyecting from epsg = 4326 to bidimensional epsg = 24891

Intersection

Aggregation

Checking CRS

Check CRS

Reproject CRS

References:

Geopandas

Coordinate Reference System

Types of projections

Universal Transverse Mercator (UTM)