Bases de Datos Geoespaciales
Lo primero que tienes que hacer es extraer el contenido del zip con el censo 2010. Para esta práctica sólo vamos a utilizar las carpetas del DF y del Estado de México.
Una vez que hayas extraido las carpetas necesitas crear una base de datos en postgres con la extensión Postgis (aquí puedes ver un tutorial sobre cómo hacerlo: Crear una base de datos espacial.
Lo primero que necesitamos hacer es cargar los shapefiles que vamos a utilizar en la base de datos que acabamos de crear (asumiremos que la base de datos se llama practica_sig). Para esto, en linux abrimos una terminal en la carpeta donde está el shape de las manzanas del DF y ejecutamos el siguiente comando(*):
shp2pgsql -c -s 4326 -I -W LATIN1 df_manzanas.shp public.merge_manzanas | sudo -u postgres psql practica_sig
(*) Para una descripción completa del cargador 'shp2pgsql' puedes consultar aquí: manpage shp2pgsql. para un instructivo más detallado: tutorial shp2pgsql
Ahora vamos a agregar la capa de manzanas del Estado de México pero en lugar de crar una nueva tabla las vamos a pegar en la tabla que ya hicimos ('public.manzanas_zmvm'):
shp2pgsql -a -s 4326 -W LATIN1 mex_manzanas.shp public.merge_manzanas | sudo -u postgres psql practica_sig'
Nota que cambiamos el switch -c por el -a, esto es para subir el shape en modo append en lugar de crear una nueva tabla. Además quitamos el modificador -I para no crear un índice espacial (el índice ya está creado y volverlo a hacer arrojaría error).
Hay otras formas más amigables de subir capas a la base de datos, por ejemplo la interfaz gráfica de shp2pgsql, búscala en el menú de inicio de windows. En escencia funciona igual que la linea de comando, utiliza esta herramienta para subir las capas de AGEBS y de calles. Recuerda que siempre vamos a querer crear una sola tabla con el merge de los datos para el DF y el Estado de México.
Finalmente, agregamos otras dos capas que están incluidas en el directorio data de este repositorio:
shp2pgsql -c -s 32614 -I -W LATIN1 estaciones_metro_final.shp public.estaciones_metro | sudo -u postgres psql practica_sig'
shp2pgsql -c -s 32614 -I -W LATIN1 limite_metropolitano.shp public.limite_metropolitano | sudo -u postgres psql practica_sig'
Nota que estas dos capas están en una proyección diferente (32614 quiere decir UTM Zona 14 Norte elipsoide WGS84), más adelante regrersaremos a esto.
Pues hemos terminado de subir los datos para la práctica 1, ahora podemos proceder a hacer las consultas y a la diversión de Postgis. Ten en cuenta que para ejecutar las consultas y visualizarlas necesitas algun cliente para la base de datos PGAdmin III provee una interfaz amigable para administrar las bases de datos y hacer consultas pero no puedes visualizar capas geográficas. QGis, provee algunas funcionalidades de consulta y administración y además permite visulizar y analizar información geográfica.
Lo primero que vamos a hacer es proyectar las geometrías de las manzanas para que sean compatibles con las estaciones del metro y el límite metropolitano. Es mejor proyectar las manzanas para que todo está en coordenadas planas, de ese modo los cálculos geométricos son mucho más rápidos:
ALTER TABLE merge_manzanas
ALTER COLUMN geom
TYPE Geometry(Polygon, 32614)
USING ST_Transform(geom, 32614);
Repite la misma operación para todas las capas que estén en coordenadas geográficas (SRID: 4326)
Ahora vamos a cortar las manzanas con el polígono del límite metropolitano (lo que se llama un clip, pues) y meter el resultado en la tabla manzanas_zmv:
create table manzanas_zmv as(
select merge_manzanas.*
from merge_manzanas
inner join limite_metropolitano on
st_intersects(limite_metropolitano.geom,merge_manzanas.geom)
)
Nota como lo que hicimos fue en realidad un inner join pero como condición utilizamos una relación espacial: st_intersects
Ahora vamos a crear un índice espacial sobre la geometría:
create index manzanas_zmvm_gix on manzanas_zmvm using GIST(geom);
TAREA: investiga qué son y para qué sirven los índices espaciales.
Antes de continuar tenemos que checar la consistencia de los datos, por ejemplo, ver si los gid son únicos:
select count(*) from manzanas_zmvm group by gid order by count(*) desc;
¿Por qué no son únicos los gids?
Como la tabla tiene unos gid's repetidos, alteramos la columna para que sean únicos:
CREATE SEQUENCE "manzanas_zmvm_gid_seq";
update manzanas_zmvm set gid = nextval('"manzanas_zmvm_gid_seq"');
Creamos los constraints necesarios y agregamos un PK:
ALTER TABLE manzanas_zmvm ALTER COLUMN "gid" SET NOT NULL;
ALTER TABLE manzanas_zmvm ADD UNIQUE ("gid");
ALTER TABLE manzanas_zmvm ADD PRIMARY KEY ("gid");
Ahora podemos empezar a hacer algunas preguntas interesantes, por ejemplo:
¿Cuántas manzanas quedan a 500 metros de cada estación del metro?
select foo.* from
(with buf as (select st_buffer(estaciones_metro.geom,500.0) as geom , estaciones_metro.nombreesta as estacion from estaciones_metro)
select count(manzanas_zmvm.gid), buf.estacion from manzanas_zmvm join buf on
st_intersects(buf.geom,manzanas_zmvm.geom)
group by buf.estacion) as foo;
¿Cuánta gente vive a 500 m de una estación del metro?
Nota: La columna pob1 contiene la población de cada manzana
select foo.* from
(with buf as (select st_buffer(estaciones_metro.geom,500.0) as geom , estaciones_metro.nombreesta as estacion from estaciones_metro)
select sum(manzanas_zmvm.pob1), buf.estacion from manzanas_zmvm join buf on
st_intersects(buf.geom,manzanas_zmvm.geom)
group by buf.estacion) as foo;
PUNTO EXTRA: ¿Cuántas personas no viven a 500 metros de una estación de metro?
Hint: Tienes que sumar el resultado de la expresión de arriba y restarla de la población total. Como es el primer quiz, puedes hacerlo en dos querys
EJERCICIO # 1.- En los datos del Censo encontrarás shapes con las calles del DF y del Estado de México, así como de las AGEBS. Agrega estos shapes como capas en Postgis y repite los pasos 1 a 6 de este archivo para obtener un corte de las calles y de las AGEBS con la forma de la ZMVM en tablas indexadas espacialmente y con llave primaria (PK)
Lo que vamos a hacer en esta parte de la práctica es completar todo un flujo de trabajo para obtener mapas de cantidad y densidad de población de habitantes por colonia.
Lo primero que vamos a hacer es relacionar las tablas de manzanas (que es donde tenemos datos de población), con la de colonias. Haciendo una unión espacial podemos asignar un identificador de colonia a la capa de manzanas:
-- Agregamos un campo para el identificador de colonia
alter table manzanas_zmv add column id_colonia int;
-- Populamos el campo a partir la unión entre manzanas y colonias
update manzanas_zmv set id_colonia = sub.id_colonia
from
(select m.id, c.id_colonia
from manzanas_zmv as m
join colonias as c
on st_intersects(c.geom, st_pointonsurface(m.geom))) as sub
where sub.id = manzanas_zmv.id
Noten como al hacer un select into, creamos implícitamente la tabla manzanas_colonias.
EJERCICIO # 2: Crea un índice espacial sobre la geometría y agrega una llave primaria a la tabla que acabas de crear.
EJERCICIO # 3: De la misma forma en que creamos la tabla manzanas_colonias, crea una tabla que una la geometría de las calles con el id de las colonias (recuerda crear su propio índice espacial y llave primaria).
Fíjense que ahora tenemos una manera de unir las tablas de manzanas y de colonias, de hecho, sería relatívamente fácil agregar la población de las manzanas de acuerdo al identificador de colonia y así obtener la población en cada colonia. Lo que vamos a hacer es un poco distinto: para acabar con mapas más bonitos, vamos a conservar la geometría de las manzanas pero agregada por las colonias, es decir, vamos a crear multipolígonos para representar a las colonias:
select st_union(geom) as geom , id_colonia
into manzanas_union
from manzanas_zmv
group by id_colonia;
Aquí estamos usando una función agregada (noten la cláusula group by), para unir la geometría de las manzanas en un sólo multipolígono para cada colonia.
Ahora sí, vamos a pegarle a nuestra tabla manzanas_unionlos datos de población:
alter table manzanas_union add column poblacion int;
update manzanas_union set poblacion = sub.poblacion
from (select sum("POB1") as poblacion, id_colonia
from manzanas_zmv
group by id_colonia) as sub
where manzanas_union.id_colonia = sub.id_colonia
Como puedes ver, esta es una consulta bastante más complicada, pero en realidad lo que hacemos es relativamente fácil si lo lees de adentro hacia afuera: en el query más interno unimos las capas de manzanas y de colonias, en el query externo agregamos por colonia y luego unimos todo a la geometría de la capa que hicimos con la unión geométrica de las manzanas.
Lo único que hace falta ahora es calcular la densidad de población de cada colonia. Para eso vamos a agregar una columna en donde vamos a calcular el area de cada colonia:
alter table manzanas_union add column densidad_poblacion float
Y calcular la densidad de población:
update manzanas_union set densidad_poblacion = poblacion/st_area(geom)
Finalmente, haz unos mapas con lo que obtuvimos.