Ciencia de Datos

Behind the scenes

Mi intención hoy es relatar cómo una cosa que parece sencilla en trabajo geoespacial siempre termina convirtiéndose en un proceso largo, pero que por fortuna, sabiendo programar, este proceso cotidiano no pasa de una hora. Esta narrativa es a manera de un “lo que callamos los cdds”, para que la gente que quiera entrar a ciencia de datos entienda que hay mucha más complejidad en el día a día, pero que fácilmente se sortea sabiendo conociendo los datos.

Y lo más interesante es cómo va creciendo el script para resolver particularidades de los datos.

Empieza así

Para un proyecto, un alumno necesitaba una tabla de las localidades de Veracruz con sus respectivos climas, algo así:

   NOMLOC            NOMMUN  TIPO_C         
 1 Barranquillas     Acajete Templado húmedo
 2 Colexta           Acajete Templado húmedo

lo cual sería el producto de una función sencilla como: ::: {.cell}

st_join(localidades, climas)

:::

Pero esto nos dío el error de que las capas no están en el mismo sistema de referencia de coordenadas (SRC)

Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared,  : 
  st_crs(x) == st_crs(y) is not TRUE

Revisamos la capa de climas, aparecía como que no tenía SRC:

Simple feature collection with 1695 features and 13 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 572600.4 ymin: 124.68 xmax: 4427355 ymax: 2767388
CRS:           NA

Entonces el archivo prj de la capa estaba mal configurado, en fin, como sabía que sí eran iguales solo le puse el mismo de la otra capa:

climas =  st_set_crs(climas, st_crs(localidades))
st_join(localidades, climas)

Pero había un error geométrico:

Error in `stopifnot()`:
ℹ In argument: `lengths(.predicate(x, y, ...)) > 0`.
Caused by error in `wk_handle.wk_wkb()`:
! Loop 3738 edge 61 has duplicate near loop 3747 edge 204

Entonces le pondremos más código, con la función para corregir las geometrías:

climas =  st_set_crs(climas, st_crs(localidades)) |>
  st_make_valid()
st_join(localidades, climas)

Y con eso, digamos, ya salía la cosa; sin embargo, hubo que hacer algunas cosas antes: INEGI proporciona dos tipos de capas de localidades: las que son POLYGON y las que son POINT, las primeras son las ciudades y zonas más urbanizadas y las segundas rancherías y villas más pequeñas, sin embargo, algunas existen tanto como punto y como polígono

Entonces lo primero que hubo que hacer fue quitar los puntos que estuvieran dentro de un polígono, recordando meter la capa de polígonos en un st_combine

st_filter(locs_p, st_combine(locs_a), .predicate = st_disjoint)

Y después de re-visitar el código

Cuando al día siguiente uno relee el script, se suelen ver cosas que fueron innecesarias: el primer intento fue con una capa de climas con SRC EPSG:4326, por lo cual reproyectamos las capas de localidades que estaban en una proyección cónica conforme de Lambert (LCC) que utiliza INEGI (y lamentablemente no tiene EPSG). Esto fue el origen de algunas geometrías inválidas.

facepalm