Archivo de la categoría: Monográficos

Leer archivos Excel con Python

Entrada sobre la importación de Excel con Python, un aporte que sirve para mi documentación y que es posible que sea de ayuda para muchos que se estén iniciando en el uso de Python y Pandas, aunque en este caso para la lectura del Excel usaremos tanto Pandas como la librería xlrd.

Lectura de Excel con Pandas

Lo más sencillo para importarnos en Python un Excel y crearnos un data frame de Pandas es:

import pandas as pd
archivo = 'C:/Users/Documents/ejemplo.xlsx'

df = pd.read_excel(archivo, sheetname='Hoja1')

df.describe()

La función read_excel será suficiente en el 80% de las ocasiones que realicemos esta tarea. Como es habitual en la ayuda tenéis perfectamente descritas sus posibilidades.

Lectura de Excel con xlrd

Es posible que necesitemos realizar tareas más complejas a la hora de leer archivos Excel y podemos usar xlrd. Vemos algunas de las posibilidades:

import xlrd 
 
archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 
print(hoja.nrows) 
print(hoja.ncols) 
print(hoja.cell_value(0, 0))

open_workbook nos abre el Excel para trabajar con él. Seleccionamos hojas por índice (empezando por el 0) y con la hoja seleccionada podemos ver el número de filas (nrows) o columnas (ncols). Seleccionar una celda lo hacemos con cell_value mediante índices (empezando por el 0). Otras posibilidades:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_name('Hoja1') 
for i in range(0,hoja.nrows):
    print(hoja.cell_value(i,1))

Si por ejemplo deseamos saber las cabeceras, los nombres de las columnas:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 
nombres = hoja.row(0)  
print(nombres)

Y mediante xlrd podemos crear data frames de pandas con lo que es posible realizar lecturas de rangos:

archivo = 'C:/Users/rvaquerizo/Documents/ejemplo.xlsx'
  
wb = xlrd.open_workbook(archivo) 

hoja = wb.sheet_by_index(0) 

# Creamos listas
filas = []
for fila in range(1,hoja.nrows):
    columnas = []
    for columna in range(0,2):
        columnas.append(hoja.cell_value(fila,columna))
    filas.append(columnas)

import pandas as pd
df = pd.DataFrame(filas)
df.head()

Hay alguna librería que lo hace de forma más elegante pero la importación de rangos de Excel con xlrd, una vez te familiarizas, me parece bastante sencilla. Espero que sea de utilidad

Los pilares de mi simulación de la extensión del COVID19

No debería publicar esta simulación de la extensión del CODVID10 o coronavirus porque puede disparar alarmas, provocar insultos, levantar ampollas,… el caso es que yo llevo 7 días de aislamiento más que el resto de España porque sólo había que ver los datos de Italia para saber lo que iba a pasar y no avisé a nadie para no disparar alarmas, provocar insultos, levantar ampollas… Y AL FINAL YO TENÍA RAZÓN. Así que os voy a exponer el motivo por el cual estoy muy asustado, bueno, hoy quiero mostraros el inicio de una simulación mala y sin fundamento que estoy realizando sobre la extensión en España del COVID19. Para hacerla vamos a emplear la siguiente información:

Y allá voy a comentaros que estoy montando. Se trata de poner a los 47 millones de españoles en una tabla, situarlos en unas coordenadas y dadas 5 personas iniciales ver como se propaga el virus municipio a municipio y, en 98 días, determinar cuantas personas pueden estar contagiadas, cuantas enfermas, cuantas sanas o cuantas desgraciadamente muertas. Esto no es que tenga lagunas, es que está inventado, pero no os creáis que las cifras oficiales son más fiables. Evidentemente, lo voy a hacer con R y dplyr. No lo subo a git porque el equipo que uso tiene un usuario de github que no es el adecuado, pero ya sabéis que el código está a vuestra disposición.

Creación de la tabla de personas edad

library(tidyverse)
library(pxR)
library(sqldf)

#detach("package:dplyr", unload=TRUE)

censo = "C:\\temp\\personales\\covid19\\0001.px"
datos <-  read.px(censo)
datos <- datos$DATA[[1]]
names(datos) = c("rango_edad", 'seccion', 'sexo', 'habitantes' )
datos <- data.frame(lapply(datos, as.character), stringsAsFactors=FALSE)

muestra <- datos %>% mutate(habitantes=round(as.numeric(habitantes)/10,0)) %>% 
  filter(seccion != "TOTAL" & sexo == "Ambos Sexos" & rango_edad != "Total") %>% select(-sexo) %>% 
  mutate(rango_edad = case_when( 
    rango_edad %in% c('0-4', '5-9', '10-14', '15-19', '20-24') ~ '<25',
    rango_edad %in% c('80-84', '85-89', '90-94', '95-99', '100 y más') ~ '80 >',
    TRUE ~rango_edad  )) 

muestra <- muestra %>% group_by(seccion,rango_edad) %>% summarise(habitantes=sum(habitantes))


espania <- muestra %>% group_by(seccion,rango_edad) %>% expand(count = seq(1:habitantes)) %>% as_tibble()

Nota: si no funciona la creación de la muestra hacéis detach de dplyr

Leemos los datos del censo que nos hemos descargado del INE, es un fichero px pero con el paquete pxR podemos manejarlo. Los datos que tenemos están a nivel de sección censal, rangos de edad, sexo y disponemos del número de habitantes. Con esta tabla de frecuencias generamos con expand una tabla que repite un registro tantas veces como digamos en una variable, es decir, repetirá la edad, el sexo, la sección censal tantas veces como habitantes tenga. Manejo una muestra del 10% porque el tema tiene un tiempo importante de procesamiento. Con esto también hago un llamamiento por si Amazon, Microsoft o Google pueden poner buenas máquinas en manos de los gestores de información (mal llamados científicos de datos ahora) de forma altruista. En fin, tenemos a todos los españoles, ahora vamos a ubicarlos con la cartografía por sección censal del INE.

Creación de la tabla de centroides municipal

library(maptools)
library(sf)
ub_shp = "C:\\temp\\mapas\\Seccion_censal\\SECC_CPV_E_20111101_01_R_INE.shp"
seccion_censal = readShapeSpatial(ub_shp)

mapa <- map_data(seccion_censal)

centroides <- mapa %>% group_by(OBJECTID = as.numeric(region)) %>% 
  summarise(centro_long=mean(long), centro_lat=mean(lat))

ggplot(data = centroides, aes(x = centro_long, y = centro_lat, group = 1)) +
  geom_polygon() 

secciones <- seccion_censal@data %>% mutate(seccion=as.character(CUSEC), municipio=as.character(CUMUN)) %>%
  select(OBJECTID,seccion,municipio)

municipios <- left_join(secciones,centroides) %>% group_by(municipio) %>% 
  summarise(centro_long=mean(long), centro_lat=mean(lat)) %>%
  select(municipio, centro_long, centro_lat)  

#Matriz de distancias
distancias <- sqldf(" select a.municipio, b.municipio as municipio2, 
                    a.centro_long, a.centro_lat, b.centro_long as centro_long2, b.centro_lat as centro_lat2
                    from municipios a , municipios b where a.municipio != b.municipio")

distancias <- distancias %>% mutate(distancia=sqrt((centro_long - centro_long2)**2 + (centro_lat-centro_lat2)**2))

Os habéis desgarcado el shapefile con las secciones censales de España y con ella calculamos el centroide de cada municipio, también he calculado una matriz de distancias porque, como veréis más adelante, la distancia de desplazamiento puede ser interesante para determinar como se mueve y como se expande el virus. En este punto está mi otra de mis reclamaciones, las compañías de telefonía podían ofrecer datos de movilidad para ayudarnos y controlar el movimiento de personas.

En fin, si cruzáis ambas tablas empieza la simulación (de mierda):

#Proceso
indices <- sample( 1:nrow( espania ), nrow(espania)/2)
espania2 <- espania[indices,]
espania2 <- espania2 %>% left_join(secciones) %>% 
  mutate(id_persona=row_number(),
         dia=0,contagiado=0, evolucion_dias=0)

sanos = espania2
contagiados = espania2[0,]
enfermos = espania2[0,]
curados = espania2[0,]
muertos = espania2[0,]

#Dia 1
#5 contagiados
dia <- sample_n(filter(espania2,seccion %in% c('2807908161','0810205003')) , 5)
contagiados <- inner_join(dia, select(sanos,id_persona)) %>% 
  mutate(contagiado=1)
lista_contagiados = unique(contagiados$id_persona)
sanos <- sanos %>% filter(id_persona %notin% lista_contagiados)

max_distancia =max(distancias$distancia,na.rm = T)

Tenemos una tabla con la población española por edad y ubicación, son 5 personas al azar de Igualada y Madrid las que empiezan todo… Veré si me atrevo a seguir contando porque lo que sigue me lo he inventado completamente.

Gráficos de calendarios con series temporales

Cuando se realizan gráficos de series temporales se emplean gráficos de líneas donde el eje X contiene la fecha y el eje Y contiene el valor a representar. Hoy quiero traer al blog otra forma de representar series temporales, los gráficos de calendario y su realización con R. Para ilustrar el ejemplo vamos a emplear las cotizaciones históricas del índice bursatil IBEX35:

require(quantmod)
require(ggplot2)
require(reshape2)
require(dplyr)
library(lubridate)

# Obtenemos las cotizaciones del IBEX 35 desde 2010
getSymbols('^IBEX', from = '2010-01-01')

# data frame de trabajo
df<-data.frame(date=index(IBEX),IBEX)

Mediante quantmod extraemos las cotizaciones del IBEX desde 2010 y creamos un data frame de trabajo que llamamos df. Vamos a realizar dos tipos de gráficos, un mapa de calor por años, meses, semanas y días y un calendario de un año puntual.

Calendario como mapa de calor por

Este es un gráfico basado en un trabajo anterior (¡de 2012!) y es una forma imaginativa de representar el cierre del IBEX 35 desde 2010 en una sola imagen. El primer paso será crear las variables a representar en el mapa de calor, el mes, el día de la semana y la semana dentro del mes:

df <- df %>% mutate(año=year(date),
                    mes=factor(month(date),levels=(1:12),
                               labels = c("ENE","FEB","MAR","ABR","MAY","JUN","JUL",
                                              "AGO","SEP","OCT","NOV","DIC"),ordered = T),
                    dia=factor(wday(date)-1,levels=rev(1:7),
                          labels=rev(c("L","M","X","J","V","S","D"))),
                    semanames=ceiling(day(date) / 7))

Ahora sólo queda representar el gráfico mediante ggplot2 donde los paneles de facet_grid serán los años en eje X y los meses en eje Y:

# Realizamos el calendario
calendario1<- ggplot(df, aes(semanames, dia, fill = IBEX.Adjusted)) + 
  geom_tile(colour = "white") + facet_grid(año~mes) + 
  scale_fill_gradient(low = "red", high = "darkgreen", na.value = "black") +  
  labs(title="Cierre histórico del IBEX", x ="Semana del mes", y = "")
calendario1

Un gráfico que me gusta bastante y una original forma de representar series temporales muy largas, no he usado paletas de colores pero imagino que los resultados mejorarán, podéis aportar esas mejoras en los comentarios.

Calendario con openair y calendarPlot

Si deseamos representar un calendario de un año concreto tenemos la función calendarPlot de openair (que me ha costado instalar en Ubuntu) que no puede ser más sencilla:

library(openair)
calendarPlot(df, pollutant = "IBEX.Adjusted", year = 2019, cols = "Greens")

Este último calendario no lo he usado pero la sintaxis es muy sencilla y el resultado queda bastante bien. Ahora vosotros mismos podéis juzgar si hay o no hay rally de fin de año.

Me rindo, es necesario trabajar en Agile

Imagen de previsualización de YouTube

“Agile sounds good” y representa todo eso que critico. Tenía compuesta y preparada una canción que versiona el “Me cago en el amor” de Tonino Carotone, “Me cago en el Agile” se llamaba. ¿Por qué este cambio de opinión tan radical?  Porque no se trabaja de forma horizontal, se trabaja de forma vertical y cada uno hace la guerra por su cuenta. Me voy a mi terreno Agile Analytics

Echamos y vendemos humo. Empecemos: arquitectos, ingenieros, científicos de datos,  analistas de negocio, diseñadores,…  y al final nos olvidamos que estamos ahí para ganar € o hacer ganar a nuestra organización €, €€€€€€ no estamos para conectar nodos edge, ingestar data lakes con millones de registros que nadie usa, diseñar algoritmos, crear contenedores, APIs, visualizaciones espectaculares, mapas interactivos, inteligencias artificiales, etecé, etecé. ¿Objetivo final del proyecto? ¿Cuántos € retorna? Bueno pues ha sido necesario de crear una parafernalia con pizarras y posit para que no nos desviemos de esos objetivos y que en último término sean capaces de medir cuantos € supone cada proyecto. Si todos los implicados en un proyecto trabajaran de forma conjunta, ¿serían necesarias estas figuras? No, pero como eso no pasa, no es que sea necesario, es que se torna imprescindible trabajar en Agile.

Me gustaría plantear una visión distinta de las personas y los roles dentro de un proyecto analítico:

  • Jefe de proyecto (project manager sounds better). Esto ni Agile, ni cascadas, hace falta alguien que se las lleve y que reparta, no queda más remedio. Es muy importante su rol en la documentación y entregables de resultados.
  • Científico, ingeniero o arquitecto de soluciones. Al final son los que trabajan tutelados por un jefe de proyecto.
  • Sponsor de negocio. Este es el que plantea el caso de negocio,  el que pone en marcha toda la maquinaria para plantear un proyecto.
  • Usuario de negocio. Persona que en último término va a usar las herramientas que estén desarrollando y el que tiene un problema que la analítica puede resolver. La línea entre sponsor y usuario es fina.
  • Product Owner (la prefiero en english). Esta figura es la clave porque es la verdadera figura horizontal a las anteriores, rol multidisciplinar tiene que entender el lenguaje de todas las personas implicadas y velar por la correcta administración del proyecto, el que se traga la correcta documentación e imputación de las tareas, que si el resto lo hiciera bien no sería necesario. Junto con el jefe de proyecto cuida del cumplimiento de los objetivos a corto plazo.
  • Scrum máster. Es la persona encargada de que todo lo anterior se lleve a cabo y para ello dispone de una serie de ceremonias y de ritos ridículos pero que sin ellos todas las figuras anteriormente citadas harían lo que les viene en gana. Tiene una tarea de interlocución a más alto nivel en la organización muy importante ya que en el caso de producirse retrasos, problemas o impedimentos en el desarrollo de los proyectos debe ser capaz de eliminar esos problemas. También mide junto con negocio y dirección cual es el impacto económico de los proyectos llevados a cabo.

Esta es mi visión de un proyecto en Agile Analytics y alguno llevo encima (con cierto éxito).  No entro en las ceremonias, que tienen su pena y su gloria también pero obligan a que todos los elementos participantes de un proyecto se reúnan, se escriban actas y se gestionen tareas. Si alguien quiere hablaré sobre ellas.

Ahora bien, hace 20 años yo trabajaba sentado con Ana (de Negocio), tenía 2 pcs porque con uno ejecutaba Business Object y si no me salía la cucaracha con el otro realizaba los análisis con Access/Excel. Vale que vendíamos politonos, paquetes de SMS y móviles de concha pero no era necesaria toda esta parafernalia y éramos 2. Podemos opinar sobre la complejidad de los problemas de negocio actuales, la información disponible o el entorno competitivo son argumentos para trabajar de este modo, pero en mi opinión el problema está en la lejanía y la especialización.  Al final no salían las cosas, así que algunos listos han vestido de metodología (que poco me gusta esta palabra) el sentido común y han dado una solución a este problema y con un halo marketiniano muy molón que llamamos Agile. A todo esto la mayor beneficiada ha sido 3M que en 2018 estaba en máximos de cotización, creo que ahora está cayendo un poco quizá haya vuelto el sentido común, no creo, hasta yo me he rendido.

Obteniendo los parámetros de mi modelo GAM

Vimos como los modelos GAM iban más allá del GLM porque en el momento de obtener los parámetros asociados al modelo de un factor nos proponían, en vez de una función lineal una función de suavizado no paramétrica para aquellos factores susceptibles de transformar en variables numéricas ordinales con un sentido determinado. Se trabajó con un modelo de riesgo con una sola variable como era la edad y al sumarizar el modelo no era posible obtener los parámetros en la salida. En último término nuestra intención con este tipo de modelos es obtener esos parámetros para transformarlos en relatividades. Qué sentido tiene obtener un buen modelo para Negocio si su resultado no se puede expresar en términos de incrementos o descuentos, en términos de relatividades.

La entrada del blog que ahora os propongo nos permite extraer los parámetros de cualquier modelo GLM o GAM a partir de la función predict y una de las opciones más olvidadas por todos nosotros:

predict(modelo, newdata = datos,  type = "terms")

con type = “terms” lo que obtenemos en el momento de realizar la predicción son los parámetros del modelo que aplicamos, no es el resultado de la predicción.

Obteniendo las relatividades de nuestro modelo GAM

Partimos del ejemplo que estamos manejando en la serie de entradas:

library(dplyr)

varib <- c(edad = 2L, sexo = 1L, zona = 1L, clase_moto = 1L, antveh = 2L,
           bonus = 1L, exposicion = 8L, nsin = 4L, impsin = 8L)

varib.classes <- c("integer", rep("factor", 3), "integer",
                   "factor", "numeric", rep("integer", 2))

con <- url("https://staff.math.su.se/esbj/GLMbook/mccase.txt")
moto <- read.fwf(con, widths = varib, header = FALSE,
                 col.names = names(varib),
                 colClasses = varib.classes,
                 na.strings = NULL, comment.char = "")


library(mgcv)

moto$edad_numero <- as.numeric(moto$edad)

gam.1 <- gam(nsin ~ s(edad_numero,bs="cr",k=3) + zona, data=filter(moto,exposicion>0), 
             offset = log(exposicion), family = poisson(link='log'))
summary(gam.1)

Ejecutad este código y obtendréis un modelo GAM con la zona por la que circula el riesgo y una función de suavizado de la edad del asegurado. A la hora de sumarizar el modelo para la edad, la variable suavizada, no vemos parámetros solo una función, si queremos obtener parámetros solo aparece la zona, ¿cómo puedo obtener las relatividades que me arroja este modelo? Empleando predict como se indicó con anterioridad:

terminos <- data.frame(exp(predict(gam.1, newdata = moto, type = "terms")))
names(terminos) <- c("rela_zona","rela_edad") 
terminos <- cbind.data.frame(terminos,select(moto,zona,edad))

Se crea el data.frame terminos que tiene el exponencial del parámetro asociado a ese registro para los factores participantes en el modelo. Cabe señalar que predict no respeta el orden de las variables en el modelo, primero pone las variables que no están suavizadas y después las suavizadas. Después de obtener los parámetros registro a registro lo que hacemos es añadir al data frame los factores de los que deseamos obtener las relatividades y como os podéis imaginar la tabla de relatividades finalmente es el resultado de seleccionar los distintos elementos:

rela_zona <- distinct(select(terminos,zona,rela_zona))
rela_edad <- distinct(select(terminos,edad,rela_edad))

Ya sabéis, no subestiméis a predict...
 

Geometría básica con R. Triángulos, circunferencias, estrellas, distancias, ángulos,…

Trabajar con triángulos y R es bien sencillo con el paquete learnGeom. La entrada viene a cuento por una duda en lista de correo de ayuda en R que no pude ayudar a resolver por no disponer de un equipo informático en ese momento. Es un paquete que nos permite visualizar los aspectos básicos de la geometría que todos tenemos olvidada. Un ejemplo de uso sería:

#install.packages("LearnGeom")
library(LearnGeom)

x_min <- 0; x_max <- 100
y_min <- 0; y_max <- 100

CoordinatePlane(x_min, x_max, y_min, y_max)

A <-c(50,50)
B <- c(70,70)
C <- c(70,50)

triangulo <- CreatePolygon(A, B, C)
Draw(triangulo, "grey")
PolygonAngles(triangulo)

 

Fijamos un plano, en este caso de 0 a 100 en ambos ejes y sobre ese plano pintamos un polígono indicando los vértices y como resultado obtenemos un triángulo rectángulo, podemos ver los ángulos que forman los vértices también y hay otras funciones interesantes como distancias entre puntos que nos sirven para recordar a Pitágoras;

DistancePoints(A,B)
sqrt(20^2+20^2)

Llegué a este paquete por lo sencillo que resultaba obtener los ángulos entre los puntos y poderlos graficar:

angle <- Angle(A, B, C, label = TRUE)
angle <- Angle(A, C, B, label = TRUE)
angle <- Angle(B, A, C, label = TRUE)

Trazar circunferencias con dirección es otra de las posibilidades con las que estoy trabajando:

CoordinatePlane(x_min, x_max, y_min, y_max)
Draw(triangulo, "transparent")
direction <- "anticlock"
inicio = 0 
fin = 45
Arc2 <- CreateArcAngles(A, 20, inicio, fin, direction)
Draw(Arc2, "red")

Por algún motivo que desconozco mi cabeza sólo puede trabajar con la dirección contraria a las agujas del reloj, es curioso. Por último, por si alguien tiene que hacer ese tipo de estructuras geométricas podemos trazar estrellas con R fijando el inicio y el ángulo de rotación:

CoordinatePlane(x_min, x_max, y_min, y_max)
Star(A, 180, 10, color= "blue")

Esta función te lleva a otra más interesante (Scissor). Ya sabéis learnGeom un paquete de R para trabajar aspectos básicos (o no tan básicos) de la geometría. Yo esto intentando hacer un proceso que haga la vuelta perfecta para todos los circuitos automovilísticos del mundo con R. Seguramente abandone el proyecto, pero siempre es bueno compartir algún conocimiento adquirido.

Mapa de códigos postales con R. Aunque el mapa es lo de menos

Entrada para facilitar la realización de mapas de códigos postales de España con R. Todo parte del trabajo de Íñigo Flores al que ya mencionamos en otra entrada. Íñigo descargó de Cartociudad y recopiló los objetos shape file para realizar estos gráficos y los subió a su repositorio, están desactualizados pero puede ser suficiente para la realización de mapas de códigos postales. Íñigo subió en formato .zip todos los archivos necesarios provincia a provincia como lo tenía Cartociudad. Podemos clonarnos el repositorio o leer directamente de github, en cualquier caso necesitamos una función en R que nos permita leer archivos comprimidos en formato zip y cuando lea el zip seleccionar que expresamente lea el archivo shp que contiene el spatial data.

Función para la lectura de archivos comprimidos zip con R

library(maptools)

leer.zip <- function(archivozip) {
  zipdir <- tempfile()
  dir.create(zipdir)

  unzip(archivozip, exdir=zipdir)
  
  archivo <- list.files(zipdir)
   
  archivo <- archivo[grepl("shp",archivo)>0]
  archivo <- paste(zipdir, archivo, sep="/")
  readShapeSpatial(archivo)
} 

Esta función leer.zip permite leer archivos zip, guardarlos en un directorio temporal y posteriormente sólo lee aquel archivo extraído que en su nombre contenga el texto “shp”. Función interesante que modificada ligeramente os permitirá descomprimir cualquier archivo y leer el elemento que deseáis, además de un buen ejemplo de uso de unzip. En este punto, como comentamos antes, podemos leer directamente de github con R.

Leer archivo zip de github con R

url <- 'https://github.com/inigoflores/ds-codigos-postales/raw/master/archive/42605-NAVARRA.zip'

tf = tempfile(tmpdir=tempdir(), fileext=".zip")
download.file(url, tf)
navarra <- leer.zip(tf)

Creamos un temporal para descargarnos el zip pero es necesario especificar la extensión. Descargamos de la url correspondiente el archivo con los elementos comprimidos y el objeto navarra será el resultado de la lectura del shapefile con los códigos postales de Navarra. La otra forma es clonar el repositorio y acceder directamente al directorio:

navarra <- leer.zip('C:\\temp\\personales\\wordpress\\ds-codigos-postales-master\\archive\\42605-NAVARRA.zip')

Otro de los motivos de esta entrada es mostraros como podemos realizar mapas de modo rápido con la librería tmap.

Ejemplo de mapa quick & dirty con R

library(tmap)

navarra <- leer.zip('C:\\temp\\personales\\wordpress\\ds-codigos-postales-master\\archive\\42605-NAVARRA.zip')
navarra@data$dato <- rpois(nrow(navarra@data),2)
qtm(shp = navarra, fill = "dato", fill.palette = "Blues")

La función qtm se traduce como -Quick thematic plot- y quick es muy quick. El mejor balance entre rápido y sencillo que hay (bajo mi punto de vista). En el ejemplo se pinta un dato aleatorio pero podéis hacer una left join con vuestros datos (que me conozco a algunos). Y si queremos crear un objeto con cada uno de los elementos que preparó Íñigo podemos hacer.

Lectura de archivos y creación de data frame mediante un bucle que lee otro data frame

trabajo <- 'C:/temp/personales/wordpress/ds-codigos-postales-master/archive/'
provincias <- list.files(trabajo)
provincias <- data.frame(archivo=provincias)
provincias$nombre <- substr(provincias$archivo,regexpr("-",provincias$archivo)+1,length(provincias$archivo))
provincias$nombre <- gsub('.zip','',provincias$nombre)

for (i in 1:nrow(provincias) ){
  instruccion <- paste0(provincias$nombre[i],' <- leer.zip("',trabajo,provincias$archivo[i],'")')
  eval(parse(text=instruccion))
}

Código rudimentario que crea un data frame a partir de los archivos de un directorio de trabajo, los archivos son los .zip que nos clonamos de github y con ellos vamos a crear 52 data frame para cada una de las provincias. El nombre de los archivos es XXXX-provincia.zip por eso tenemos que usar algunas funciones de texto para obtener el nombre de la provincia como regexpr que nos permite encontrar la primera posición en la se encuentra un patrón dentro de un texto, por otro lado gsub nos sirve para sustituir un patrón de texto por otro. Así leemos desde el – y posteriormente tenemos que eliminar el .zip para tener el nombre de cada provincia. Y por último un clásico en mis programas de R herencia de los tiempos en los que trabajaba con macros en SAS, tengo que recorrer ese data frame con los elementos del directorio y el nombre del objeto será una columna del data frame y el archivo a leer otra columna, para evaluar un texto el mítico eval ( parse ( text = nunca me falla, habrá formas más elegantes pero esta son dos líneas. Siempre hay que poner talento en la construcción de la instrucción y acordarse de cerrar paréntesis y demás. Ejecutando eso tendríamos un objeto para cada provincia, si queremos toda España.

Creación de un mapa de España por códigos postales

trabajo <- 'C:/temp/personales/wordpress/ds-codigos-postales-master/archive/'
provincias <- list.files(trabajo)
provincias <- data.frame(archivo=provincias)
provincias$nombre <- substr(provincias$archivo,regexpr("-",provincias$archivo)+1,length(provincias$archivo))
provincias$nombre <- gsub('.zip','',provincias$nombre)

for (i in 1:nrow(provincias) ){
  instruccion <- paste0('borrar <- leer.zip("',trabajo,provincias$archivo[i],'")')
  eval(parse(text=instruccion))
  if (i == 1) {espania <- borrar}
  espania <- rbind(espania,borrar)
  remove(borrar)
}

plot(espania)

Otro bucle con la marca de la casa pero que funciona perfectamente, leemos uno a uno cada zip con las provincias y con rbind podemos unir los objetos spatial para poder pintar el mapa de España y cuidado que esto si genera un objeto de casi 120 MB. Podéis manejar los objetos spatial data y así reducir su tamaño, así como idea por si pongo el código en el repositorio.

El caso es que ya sabéis como hacer un mapa de España de códigos postales con R, incluso si sois avezados podéis guardar el objeto final resultante y utilizarlo con QGIS u otra herramienta que uséis para hacer mapas. Además esta entrada es todo un compendio de malas prácticas en programación con R que funcionan a las mil maravillas, desde leer archivos zip con R, seleccionar el que deseamos a funciones de texto para extraer con condiciones, ejemplo de gráfico de mapa rápido con tmap y un bucle que lee un data frame y genera objetos con él.

Guardar objeto de R como shapefile (shp)

 writeSpatialShape(espania, "C:/temp/personales/wordpress/espania.shp")

Por último podemos guardar el objeto resultante de R para usarlo directamente con QGIS, se generan todos los archivos necesarios, el shp, el dbf y el otro.

Inteligencia Arficial frente a un juego de niños. La partícula tonta de Nicolás

Pablo Picasso decía que en aprender a pintar como los pintores del renacimiento tardó unos años pero pintar como los niños le llevó toda la vida y en ocasiones creo que hacemos las cosas difíciles porque nos creemos que hacemos cosas difíciles y entonces llega un niño de nueve años y dice “Papá un punto que primero vaya a la izquierda y luego a la derecha no es tan difícil”.
Os pongo en antecedentes, el pasado 7 de mayo fui al AWS Summit de Madrid porque Sergio Caballero iba a contar uno de los casos de uso. Los de AWS no se deben ni imaginar de las maravillas que ha hecho Sergio en el Ayuntamiento de Alcobendas porque sólo dejaron que hablara 10 minutos, muy torpes ellos, su trabajo es mejor escaparate que el planteado por Mai-Lan Tomsen, un error en el planteamiento de la jornada. El caso es que había una “competición” de vehículos que circulaban por un circuito guiados por complicados algoritmos de inteligencia artificial. Vimos algún “bucanero serio” de alguno de los participantes, ya sabemos reinforcement learning, pero reinforcement reinforcement. Otros participantes más o menos honrosos, en fin, distraído. Viendo la competición me entraron ganas de participar y al llegar a casa me siento a preparar un algoritmo que recorriera el circuito del Jarama de Madrid, no un circuito cualquiera un circuito donde yo he visto ganar carreras a Jorge Martínez Aspar.

Portátil y R, empiezo mi trabajo con imager, busco en la Wikipedia el circuito, lo cargo, genero un data frame, selecciono puntos y comienzo a diseñar mi propia estrategia de reinforcement learning combinadas con técnicas de machine learning, algo como “SVM direccionables” se acerca por detrás mi hijo y me suelta “Papá un punto que primero vaya a la izquierda y luego a la derecha no es tan difícil”. Bueno, pues en 20 minutos sale esto:

De momento no funciona pero no me digáis que no es genial la idea, lo que hace con pocas líneas de código y una consulta en sql. En el repositorio de analisisydecision tenéis el código en R que realiza esta maravilla, he llamado al código partícula tonta y tiene aspectos interesantes en cuanto al uso de la librería imager de R para el tratamiento de imágenes y como transformo una imagen en un data frame de coordenadas y por supuesto la genial idea de Nicolás.

Por cierto, al ver el resultado Nicolás dijo que no sólo derecha e izquierda, también era necesario un arriba y abajo. Tengo abandonado el proyecto, como muchos, pero la anécdota me ayudó en mi trabajo.

Lectura de archivos csv con Python y Pandas

A continuación os planteo un acercamiento básico a la lectura de archivos csv con Python y algunos trucos para facilitar la vida cuando realizamos importaciones basados en la experiencia como son leer los primeros registros del csv o realizar una lectura de observaciones aleatoria por si el archivo es muy voluminoso. Para realizar las importaciones vamos a emplear Pandas y la función read_csv con sus infititas opciones:

 pd.read_csv(filepath_or_buffer, sep=', ', delimiter=None, header='infer', names=None, index_col=None, usecols=None, squeeze=False, prefix=None, mangle_dupe_cols=True, dtype=None, engine=None, converters=None, true_values=None, false_values=None, skipinitialspace=False, skiprows=None, nrows=None, na_values=None, keep_default_na=True, na_filter=True, verbose=False, skip_blank_lines=True, parse_dates=False, infer_datetime_format=False, keep_date_col=False, date_parser=None, dayfirst=False, iterator=False, chunksize=None, compression='infer', thousands=None, decimal=b'.', lineterminator=None, quotechar='"', quoting=0, escapechar=None, comment=None, encoding=None, dialect=None, tupleize_cols=None, error_bad_lines=True, warn_bad_lines=True, skipfooter=0, doublequote=True, delim_whitespace=False, low_memory=True, memory_map=False, float_precision=None)

Para trabajar la entrada vamos a necesitar dos archivos de texto:

Como costumbre poner la ubicación del archivo y después la lectura:

path = 'C:/temp/'

import pandas as pd
df = pd.read_csv (path + 'index.csv')
df.head()

En este caso la vida es maravillosa y ha salido todo a la primera pero sabemos que eso no pasa siempre, ejecutáis:

df = pd.read_csv (path + 'bank-additional-full.csv')
df.head()

El separador es distinto:

df = pd.read_csv (path + 'bank-additional-full.csv', sep = ';')
df.head()

La vida sigue sin ser muy complicada porque el archivo de ejemplo tiene pocos registros, pero imaginad que leéis unas docenas de GB por ello previamente es mejor ejecutar:

df = pd.read_csv (path + 'bank-additional-full.csv', nrows= 200)
df.shape

con nrows = 200 leemos las primeras 200 líneas y podemos comprobar si lo estamos leyendo correctamente y podemos ahorrarnos disgustos, tiempo y trabajo. E incluso estaría bien no leer las docenas de GB porque no tenemos suficiente memoria o porque no necesitamos leer entero el archivo podemos leer por trozos:

meses = ['may', 'jul']
df = pd.DataFrame()
for trozo in pd.read_csv(path + 'bank-additional-full.csv', sep=';',
                             chunksize=1000):
    df = pd.concat([df,trozo[trozo['month'].isin(meses)]])

df.month.value_counts()

Con chunksize estamos leyendo el archivo csv en trozos (chunks) de 1000 en 1000 y nos quedamos sólo con aquellos que cumplan un determinado requisito, en este caso que el campo month sea may o jul. E incluso podéis leer el csv extrayendo una muestra aleatoria mientras leéis el fichero por partes y no sobre pasar la memoria:

df2 = pd.DataFrame()
for trozo in pd.read_csv(path + 'bank-additional-full.csv', sep=';',
                             chunksize=1000):
    df2 = pd.concat([df2,trozo.sample(frac=0.25)])
df2.shape

Este último truco puede servir para leer csv extremadamente grandes y realizar los primeros análisis aproximativos a nuestro problema porque como dice un buen amigo “si en 200.000 registros no encuentras una señal no hace falta que cargues millones”.

Longitud de las frases del Quijote con #rstats

Siempre he querido hacer cosas con Rstats y el Quijote y ayer se me ocurrió medir la longitud de las frases del Quijote y crear un histograma que describa esta longitud. Aunque confieso que no me lo he leído, me he quedado en el capítulo 7 u 8 (no recuerdo) el caso es que me pareció hipnótico con sus ritmos, es musical. Además tengo muchas ganas de meter mano al proyecto Gutemberg porque esos ritmos, esa musicalidad, el uso de palabras esdrújulas,… me llama la atención.
Bueno, al lío, todo el código está subido al repositorio por si lo queréis, pero hay algunas funciones y algunas ideas que me parecen interesantes.

library(dplyr)
library(ggplot2)
library(plotly)

#Leemos el fichero desde proyecto Gutemberg
ubicacion <- "https://www.gutenberg.org/cache/epub/2000/pg2000.txt"
quijote <- read.table (ubicacion,sep="\r", encoding="UTF-8")
quijote <- data.frame(quijote)
names(quijote) <- 'linea'

#Transformaciones e identificar el inicio del libro.
quijote <- quijote %>%
  mutate(linea = toupper(linea),
         inicio = grepl("EN UN LUGAR DE LA MANCHA",linea)>0)

Leemos directamente un txt desde Gutemberg y prefiero transformarlo en data frame para usar dplyr. Todas las palabras las pongo en mayúsculas e identifico donde empieza el Quijote, para evitar prólogos y demás. Ya tengo unos datos con los que poder trabajar:

#Marcamos lo que vamos a leer
desde <- which(quijote$inicio)
hasta <- nrow(quijote)

#Texto de trabajo
texto <- quijote[desde:hasta,1]

#El texto lo transformamos en una lista separada por espacios
texto_split = strsplit(texto, split=" ")

#Deshacemos esa lista y tenemos el data.frame
texto_col = as.character(unlist(texto_split))
texto_col = data.frame(texto_col)
names(texto_col) = 'palabra'

En este caso los datos los quiero de tal forma que disponga de un data frame con una sola variable que sea cada palabra del Quijote. Ahora voy a medir las frases identificando donde hay puntos en esas palabras:

#Identificamos donde tenemos puntos y un autonumérico del registro
texto_col <- texto_col %>% filter(!is.na(palabra)) %>%
  mutate(punto = ifelse(grepl('.',palabra,fixed=T),"FIN","NO"),
         posicion = row_number())

¿Qué se me ha ocurrido? Trabajar con autonuméricos, tengo identificados los puntos, ahora tengo que fijar una posición inicial y una posición final:

#Si unimos las posiciones con puntos con lag podemos calcular la longitud
pos_puntos1 <- filter(texto_col,punto=="FIN") %>% 
  select(posicion) %>% mutate(id = row_number())

pos_puntos2 <- pos_puntos1 %>% mutate(id = id + 1) %>%
  rename(posicion_final = posicion)

pos_puntos <- left_join(pos_puntos1,pos_puntos2) %>%
  mutate(longitud = ifelse(is.na(posicion_final), posicion, posicion - posicion_final))

Como no soy un tipo muy brillante opto por una opción sencilla de cruzar una tabla consigo misma, como me ponen los productos cartesianos “con talento”. La idea es seleccionar solo los registros que marcan el final de la frase, un autonumérico me marca cual es cada frase, ahora si hago una left join por el id de la frase y el id + 1 de la frase creo una especie de lag. La longitud de la frase será donde está el punto menos donde estaba el final de la anterior frase. Creo que me he explicado de pena, pero si veis el data frame final lo entenderéis mejor. Ahora ya pinto un histograma:

#GRaficamos la longitud
plot_ly(data = pos_puntos, x = ~longitud, type = "histogram") %>%
  layout(title = "Longitud de las frases del Quijote",
         xaxis = list(title = "Longitud"), yaxis = list(title = ""))

Y queda una gamma perfecta, yo diría que hasta bonita. Ahora quedaría identificar los parámetros de esta gamma y compararlos con otros libros, e incluso comparar lenguas. Pero esas tareas se las dejo a los “buenos”.