Extraer población de una cobertura de radio FM o de televisión (INEGI 2020)

Extraer población de una cobertura de radio FM o de televisión (INEGI 2020)

mario.hernandez 14 March 2021

A continuación muestro el procedimiento para saber que poblaciones del INEGI 2020 cubre una cobertura de televisión o de radio FM generada con el Sistema de Preanálisis del Instituto Federal de Telecomunicaciones - IFT.

Primero necesitamos el siguiente software libre (Puede ser descargado o utilizado gratuitamente)

  1. Sistema de Preanálisis del IFT. Genera tu cobertura de radio FM o TV en la plataforma http://mapasradiodifusion.ift.org.mx/CPCREL-web;
  2. QGIS3. Descarga en https://www.qgis.org/ ;
  3. Python3. Descarga en https://www.python.org/

Dividiremos el proceso en los siguientes pasos:

  1. Descarga de las poblaciones del censo INEGI 2020;
  2. Convertir poblaciones CSV del INEGI 2020 a SHP (formato ShapeFile);
  3. Generar cobertura Radio o TV con la Herramienta de Preanálisis del IFT;
  4. Conversión de cobertura KMZ a SHP;
  5. Extracción de las poblaciones del INEGI 2020 dónde tendría servicio la cobertura generada;
  6. Exportar poblaciones a Excel y cálculo de población total.

Antes de ejecutar los scripts, deben tener instaladas las siguientes librerías de python3:

  • requests (Descargar CSV de INEGI);
  • Pillow (Procesamiento de Imágenes);

Esto lo realizan desde la terminal con pip3:

pip3 install requests
pip3 install Pillow

Y definir un directorio de trabajo. Yo lo definí en la siguiente dirección:

work_dir = '/Users/m-ario/Documents/get_pop'

Los scripts serán ejecutados desde QGIS3, en la siguiente sección (View/Panels/Processing ToolBox/Create New Script...). Además, activar la consola de python:

1. Descarga de las poblaciones del censo INEGI 2020

Para descargar todas las poblaciones de las entidades federativas de México y juntarlas en un solo archivo CSV, utilizaremos el siguiente script en python:

import requests
import zipfile
import os
import csv

def del_tmp(work_dir):
    '''
    CREA CARPETA TEMPORAL DE TRABAJO tmp DENTRO DEL DIRECTORIO DE TRABAJO
    '''
    if not os.path.exists(f'{work_dir}/tmp'):
        os.system(f'mkdir {work_dir}/tmp')
    else:
        os.system(f'rm -rf {work_dir}/tmp && mkdir {work_dir}/tmp')


def download_censo_2020(work_dir):
    '''
    DESCARGA LAS POBLACIONES DEL INGEGI 2020 Y LAS JUNTA EN UN SOLO ARCHIVO
    ingegi_2020.csv
    '''
    pop_nacional = [['ENTIDAD', 'NOM_ENT', 'MUN', 'NOM_NUM',
                     'LOC', 'NOM_LOC', 'LONGITUD', 'LATITUD',
                     'ALTITUD', 'POBTOT']]

    for i in range(33):
        del_tmp(work_dir)
        print(f'Descargando iter_{ "0" + str(i) if i < 10 else i }_cpv2020')
        r = requests.get(f'https://www.inegi.org.mx/contenidos/programas/ccpv/2020/datosabiertos/iter/iter_{ "0" + str(i) if i < 10 else i }_cpv2020_csv.zip')
        open(f'{work_dir}/tmp/{i}.zip', 'wb').write(r.content)

        with zipfile.ZipFile(f'{work_dir}/tmp/{i}.zip', 'r') as zip_ref:
            zip_ref.extractall(f'{work_dir}/tmp')

        if os.path.exists(f'{work_dir}/tmp/conjunto_de_datos'):
            conjunto_datos = f'{work_dir}/tmp/conjunto_de_datos/conjunto_de_datos_iter_{ "0" + str(i) if i < 10 else i }_cpv2020.csv'
        else:
            conjunto_datos = f'{work_dir}/tmp/iter_{ "0" + str(i) if i < 10 else i }_cpv2020'
            conjunto_datos = conjunto_datos + f'/conjunto_de_datos/conjunto_de_datos_iter_{ "0" + str(i) if i < 10 else i }_cpv2020.csv'

        with open(conjunto_datos, encoding='UTF-8') as csvfile:
            file = csv.reader(csvfile, delimiter=',')

            csv_filter = filter(lambda x: x[6] != '', file)  # Retira registros que no tienen coordenadas
            csv_filter = list(map(lambda x: x[:10], csv_filter))  # Solo las primeras 9 columnas

        pop_nacional = pop_nacional + csv_filter[1:]

    with open(f'{work_dir}/inegi_2020.csv', 'w', encoding='UTF-8') as csvfile:
        csv_writer = csv.writer(csvfile, delimiter=',')
        csv_writer.writerows(pop_nacional)


work_dir = '/Users/m-ario/Documents/get_pop'

download_censo_2020(work_dir)

Lo anterior ejecutarlo desde QGIS3:

Al finalizar la ejecución, en el directorio de trabajo tendran un archivo CSV que contiene todas las poblaciones del censo INEGI 2020.

inegi censo 2020

2. Convertir poblaciones CSV del INEGI 2020 a SHP

Una vez que se termine de ejecutar el script, hacer el siguiente procedimiento "manual", a fin de convertir el archivo CSV a SHP. Es importante que los archivos del SHP estén en el directorio de trabajo:

Cargar archivo CSV:

Guardar capa como SHP:

Cuando finalice el proceso de conversión, aseguren que los archivos asociados al formato SHP se encuentren en el directorio de trabajo:

Archivo SHP Censo INEGI 2020

3. Generar cobertura Radio o TV con la Herramienta de Preanálisis del IFT;

Para obtener la cobertura de radio, generarla con la herramienta http://mapasradiodifusion.ift.org.mx/CPCREL-web y descargar el archivo KMZ (Y colocarlo en el directorio de trabajo):

Además deben colocar el nombre del archivo KMZ en la variable kmz_file, en mi caso (Ver al final del Post en el código Final):

kmz_file = 'TDT_Canal_22.kmz'

4. Conversión de cobertura KMZ a SHP;

Para realizar esta conversión utilizo la misma metodoloía de este Post (Leerlo para más claridad de lo que realiza el código de abajo), sin embargo, se realiza de forma automática (Toda la ejecución es en serie, no obstante dejo la explicación de cada función. Al final ustedes solo veran las poblaciones que cubre la cobertura que generaron con la herramienta del IFT:

4.1 Función get_georeference(work_dir, kmz_file) . Extrae los archivos que contiene el KMZ en la carpeta tmp y obtiene la georreferencia contenida en el archivo KML:

def get_georeference(work_dir, kmz_file):
    '''
    Obtiene georreferenciacion de la cobertura generada con la herramienta del IFT
    http://mapasradiodifusion.ift.org.mx/CPCREL-web
    '''
    import re

    del_tmp(work_dir)  # Borra la carpeta temporal tmp

    with zipfile.ZipFile(f'{work_dir}/{kmz_file}', 'r') as zip_ref:
        zip_ref.extractall(f'{work_dir}/tmp')

    with open(f'{work_dir}/tmp/{kmz_file[:-4]}_i.kml', encoding='UTF-8') as kml:
        kml_info = kml.read()
        north = re.findall('north>(.*)<', kml_info)[0]
        south = re.findall('south>(.*)<', kml_info)[0]
        east = re.findall('east>(.*)<', kml_info)[0]
        west = re.findall('west>(.*)<', kml_info)[0]

    return west, north, east, south

Después de ejecutar esta función se visualizarán los siguiente archivos:

KML Georreferencia

4.2 Función pixels_to_red(work_dir, png_dir) . Convierte todos los pixeles de la cobertura a Rojo (Excepto los pixeles transparentes del fondo):

def pixels_to_red(work_dir, png_dir):
    '''
    Convierte todos los pixeles de la cobertura PNG contenida en el archivo KMZ
    a color Rojo
    '''
    img = Image.open(png_dir)
    img = img.convert("RGBA")

    pixdata = img.load()

    for y in range(img.size[1]):
        for x in range(img.size[0]):
            if pixdata[x, y] != (255, 255, 255, 0):  # No cambia Pixeles Transparentes
                pixdata[x, y] = (255, 0, 0, 255)  # Cualquier otro color lo cambia a rojo

    img.save(f'{work_dir}/tmp/out_red.png')

Después de ejecutar esta función se visualizarán los siguiente archivos:

4.3 Función kml2shp(work_dir, kmz_file, georeference) . Realiza la conversión del archivo KMZ a SHP. El paso de arreglo de errores de Geometría es necesario, ya que cuando utilizo la librería Pillow, la imagen queda en formato RGBA, es decir, con 4 bandas de colores. Al realizar la conversión de TIF a SHP lo realizo sobre la banda 1, no obstante, la imagen TIF de salida aún tiene 3 bandas (RGB), esto por una alguna razón que aún no logro comprender, provoca un error en la geometría del archivo SHP (Quizás es mejor convertirla a blanco y negro?):

def kml2shp(work_dir, kmz_file, georeference):
    '''
    Convierte archivo KMZ a SHP
    '''

    # CONVIERTE IMAGEN PNG A TIF GEOREFERENCIADA
    processing.run('gdal:translate', 
        {'INPUT': f'{work_dir}/tmp/out_red.png', 
        'TARGET_CRS': 'EPSG:4326', 
        'OUTPUT': f'{work_dir}/tmp/out.tif', 
        'EXTRA': f'-a_ullr {georeference[0]} {georeference[1]} {georeference[2]} {georeference[3]}'})
    
    # CONVIERTE IMAGEN TIF A SHP
    processing.run('gdal:polygonize', 
        {'INPUT': f'{work_dir}/tmp/out.tif', 
        'OUTPUT': f'{work_dir}/tmp/out.shp', 
        'BAND': 1})
    
    # SI EXISTEN ERRORES DE GEOMETRIA LOS ARREGLA
    processing.run('qgis:fixgeometries', 
        {'INPUT': f'{work_dir}/tmp/out.shp', 
        'OUTPUT': f'{work_dir}/tmp/out_fixed.shp'})

Después de ejecutar esta función se visualizarán los siguiente archivos:

SHP final sin errores de geometría

5. Poblaciones INEGI 2020 donde brinda servicio la cobertura;

Con la siguiente función obetnemos en un archivo CSV "cobetura_pop2020" (Puede ser abierto con excel) las poblaciones del INEGI 2020 que se encuentran dentro de la cobertura generada. Los archivos serán depositados en la carpeta out, y el nombre de la carpeta del estudio realizado será la hora de ejecución del script:

def get_pop(work_dir):
    '''
    Extrae las poblaciones contenidas en la cobertura SHP out_fixed.shp
    '''
    
    # Crea carpeta para depositar los archivos de salida
    if not os.path.exists(f'{work_dir}/out'):
        os.system(f'mkdir {work_dir}/out')

    study_name = datetime.now().strftime("%I%p-%M-%S-%f__%d-%b-%Y")
    os.system(f'mkdir {work_dir}/out/{study_name}')
    
    processing.run('qgis:clip', 
        {'INPUT': f'{work_dir}/inegi_2020.shp', 
        'OVERLAY': f'{work_dir}/tmp/out_fixed.shp',
        'OUTPUT': f'{work_dir}/out/{study_name}/cobetura_pop2020.csv'})
    
    # Copia archivos SHP al directorio de salida
    shp_files = ['.shp', '.dbf', '.shx', '.prj']
    for file in shp_files:
        os.system(f'cp {work_dir}/tmp/out_fixed{file} {work_dir}/out/{study_name}/out_fixed{file}')

Después de ejecutar esta función se visualizarán los siguiente archivos:

Poblaciones dentro de la cobertura

Para cargar los archivos finales a QGIS3:

6. Población total de la cobertura;

Existen diversos métodos, el más sencillo es abrir el archivo CSV con EXCEL y sumar las poblaciones de la columna POBTOT (Como se observa al final del video esta cobertura cubre en total 4,674,042 habitantes):

Finalmente dejo todo el código para tener el script completo. Es necesario cambiar dos variables que dependen de su directorio de trabajo y el nombre de la cobertura KMZ descargada de la pág del IFT. En mi caso estan configuradas como:

work_dir = '/Users/m-ario/Documents/get_pop'
kmz_file = 'TDT_Canal_22.kmz'

Código python3 Final (Debe estar dentro de QGIS3. Ver primer video del post)

import requests
import zipfile
import os
import csv
import time

from qgis import processing
from PIL import Image

from datetime import datetime


def del_tmp(work_dir):
    '''
    CREA CARPETA TEMPORAL DE TRABAJO tmp DENTRO DEL DIRECTORIO DE TRABAJO
    '''
    if not os.path.exists(f'{work_dir}/tmp'):
        os.system(f'mkdir {work_dir}/tmp')
    else:
        os.system(f'rm -rf {work_dir}/tmp && mkdir {work_dir}/tmp')


def download_censo_2020(work_dir):
    '''
    DESCARGA LAS POBLACIONES DEL INGEGI 2020 Y LAS JUNTA EN UN SOLO ARCHIVO
    ingegi_2020.csv
    '''
    pop_nacional = [['ENTIDAD', 'NOM_ENT', 'MUN', 'NOM_NUM',
                     'LOC', 'NOM_LOC', 'LONGITUD', 'LATITUD',
                     'ALTITUD', 'POBTOT']]

    for i in range(33):
        del_tmp(work_dir)
        print(f'Descargando iter_{ "0" + str(i) if i < 10 else i }_cpv2020')
        r = requests.get(f'https://www.inegi.org.mx/contenidos/programas/ccpv/2020/datosabiertos/iter/iter_{ "0" + str(i) if i < 10 else i }_cpv2020_csv.zip')
        open(f'{work_dir}/tmp/{i}.zip', 'wb').write(r.content)

        with zipfile.ZipFile(f'{work_dir}/tmp/{i}.zip', 'r') as zip_ref:
            zip_ref.extractall(f'{work_dir}/tmp')

        if os.path.exists(f'{work_dir}/tmp/conjunto_de_datos'):
            conjunto_datos = f'{work_dir}/tmp/conjunto_de_datos/conjunto_de_datos_iter_{ "0" + str(i) if i < 10 else i }_cpv2020.csv'
        else:
            conjunto_datos = f'{work_dir}/tmp/iter_{ "0" + str(i) if i < 10 else i }_cpv2020'
            conjunto_datos = conjunto_datos + f'/conjunto_de_datos/conjunto_de_datos_iter_{ "0" + str(i) if i < 10 else i }_cpv2020.csv'

        with open(conjunto_datos, encoding='UTF-8') as csvfile:
            file = csv.reader(csvfile, delimiter=',')

            csv_filter = filter(lambda x: x[6] != '', file)  # Registros que no tienen coordenadas
            csv_filter = list(map(lambda x: x[:10], csv_filter))  # Solo las primeras 9 columnas

        pop_nacional = pop_nacional + csv_filter[1:]

    with open(f'{work_dir}/inegi_2020.csv', 'w', encoding='UTF-8') as csvfile:
        csv_writer = csv.writer(csvfile, delimiter=',')
        csv_writer.writerows(pop_nacional)


def get_georeference(work_dir, kmz_file):
    '''
    Obtiene georreferenciacion de la cobertura generada con la herramienta del IFT
    http://mapasradiodifusion.ift.org.mx/CPCREL-web
    '''
    import re
    
    del_tmp(work_dir)  # Borra la carpeta temporal tmp

    with zipfile.ZipFile(f'{work_dir}/{kmz_file}', 'r') as zip_ref:
        zip_ref.extractall(f'{work_dir}/tmp')

    with open(f'{work_dir}/tmp/{kmz_file[:-4]}_i.kml', encoding='UTF-8') as kml:
        kml_info = kml.read()
        north = re.findall('north>(.*)<', kml_info)[0]
        south = re.findall('south>(.*)<', kml_info)[0]
        east = re.findall('east>(.*)<', kml_info)[0]
        west = re.findall('west>(.*)<', kml_info)[0]

    return west, north, east, south


def pixels_to_red(work_dir, png_dir):
    '''
    Convierte todos los pixeles de la cobertura PNG contenida en el archivo KMZ
    a color Rojo
    '''
    img = Image.open(png_dir)
    img = img.convert("RGBA")

    pixdata = img.load()

    for y in range(img.size[1]):
        for x in range(img.size[0]):
            if pixdata[x, y] != (255, 255, 255, 0):  # No cambia Pixeles Transparentes
                pixdata[x, y] = (255, 0, 0, 255)  # Cualquier otro color lo cambia a rojo

    img.save(f'{work_dir}/tmp/out_red.png')


def kml2shp(work_dir, kmz_file, georeference):
    '''
    Convierte archivo KMZ a SHP
    '''

    # CONVIERTE IMAGEN PNG A TIF GEOREFERENCIADA
    processing.run('gdal:translate', 
        {'INPUT': f'{work_dir}/tmp/out_red.png', 
        'TARGET_CRS': 'EPSG:4326', 
        'OUTPUT': f'{work_dir}/tmp/out.tif', 
        'EXTRA': f'-a_ullr {georeference[0]} {georeference[1]} {georeference[2]} {georeference[3]}'})
    
    # CONVIERTE IMAGEN TIF A SHP
    processing.run('gdal:polygonize', 
        {'INPUT': f'{work_dir}/tmp/out.tif', 
        'OUTPUT': f'{work_dir}/tmp/out.shp', 
        'BAND': 1})
    
    # SI EXISTEN ERRORES DE GEOMETRIA LOS ARREGLA
    processing.run('qgis:fixgeometries', 
        {'INPUT': f'{work_dir}/tmp/out.shp', 
        'OUTPUT': f'{work_dir}/tmp/out_fixed.shp'})
        

def get_pop(work_dir):
    '''
    Extrae las poblaciones contenidas en la cobertura SHP out_fixed.shp
    '''
    
    # Crea carpeta para depositar los archivos de salida
    if not os.path.exists(f'{work_dir}/out'):
        os.system(f'mkdir {work_dir}/out')

    study_name = datetime.now().strftime("%I%p-%M-%S-%f__%d-%b-%Y")
    os.system(f'mkdir {work_dir}/out/{study_name}')
    
    processing.run('qgis:clip', 
        {'INPUT': f'{work_dir}/inegi_2020.shp', 
        'OVERLAY': f'{work_dir}/tmp/out_fixed.shp',
        'OUTPUT': f'{work_dir}/out/{study_name}/cobetura_pop2020.csv'})
    
    # Copia archivos SHP al directorio de salida
    shp_files = ['.shp', '.dbf', '.shx', '.prj']
    for file in shp_files:
        os.system(f'cp {work_dir}/tmp/out_fixed{file} {work_dir}/out/{study_name}/out_fixed{file}')


work_dir = '/Users/m-ario/Documents/get_pop'
kmz_file = 'TDT_Canal_22.kmz'

# Si ya realizaron el proceso 1, no se vuelve a ejecutar
if not os.path.exists(f'{work_dir}/inegi_2020.csv'): download_censo_2020(work_dir)

print('OBTENIENDO GEORREFERECIA DE LA COBERTURA ..')
georeference = get_georeference(work_dir, kmz_file)

print('CONVIRTIENDO PIXELES A ROJO ..')
pixels_to_red(work_dir, f'{work_dir}/tmp/{kmz_file[:-4]}.png')

print('CONVIRTIENDO A SHP ..')
kml2shp(work_dir, kmz_file, georeference)

print('EXTRAYENDO POBLACIÓN ..')
get_pop(work_dir)