Leer imagenes .BIL del INEGI con Python

Los archivos de Banda Intercalada por Linea (BIL) es un formato de archivo binario en el que las distintas bandas de una imagen se concatenan formando una única matriz de datos. Los disponibles en la página del INEGI suelen tener una sola banda por lo que su extracción resulta aún más sencilla.

Para este ejemplo vamos utilizar fotografías aéreas. Vamos a ir a la pestaña de Datos, luego Mapas y finalmente a Topografía que nos mostrará este mapa. Seleccionamos la escala de 1:20,000 y nos aparecerán cuadriculas sobre el territorio y seleccionamos la que sea de nuestro interés. Al momento en que escribo esta entrada la página tiene un bug y muestra el mensaje de "No existen mapas relacionados con esta carta" cuando se da clic en cualquier zona y se podría pensar que no hay datos disponibles para ella, pero los datos aparecerán más abajo. Voy a utilizar una fotografía que cubre parte del municipio de Jojutla de Juarez, Morelos [archivo E14A69e]. Los archivos vienen en una carpeta comprimida ZIP. La información requerida para extraer los datos (dimensiones y formato) viene en un archivo .txt. En este caso el archivo es e14a69e.txt:

	ORTOFOTO DIGITAL: E14A69E
	FUENTE: Fotografías aéreas escala 1:75,000 de Noviembre de 1995
	PROCESAMIENTO: Rectificación de fotografías aéreas, con auxilio
	de puntos de control geodésico y Modelo Digital de Elevación.
	PROYECCION: Universal Transversa de Mercator (UTM)
	DATUM: ITRF92
	ELIPSOIDE: GRS 80
	DIMENSIONES DE LA IMAGEN:
	Columnas : 5950
	Renglones: 6999
	ZONA UTM: 14
	COORDENADAS DE LA ESQUINA NOROESTE:
	Este: 476480.0
	Norte: 2059450.0
	COORDENADAS DE LA ESQUINA SURESTE:
	Este: 488380.0
	Norte: 2045452.0
	DIMENSIONES DEL PIXEL X,Y: 2.0 metros
	FORMATO: Datos binarios crudos: 1 byte por pixel

view raw e14a69e.txt hosted with ❤ by GitHub

El programa de lectura y visualización es el siguiente:

	#!/usr/bin/env python3
	# -*- coding: utf-8 -*-
	"""
	Created on Mon Jan 14 20:23:30 2019
	@author: Rodolfo Escobar
	"""
	import matplotlib.pyplot as plt
	from matplotlib.ticker import StrMethodFormatter
	import numpy as np
	#Dimensiones y tipo de datos
	X = 5950 # Columnas
	Y = 6999 # Renglones
	dim = (Y,X)
	dtype = np.dtype('uint8') # 1 byte por pixel
	# Leer .bil
	fid = open("e14a69e.bil", 'rb')
	data = np.fromfile(fid, dtype)
	image = data.reshape(dim)
	#Recorte de imagen
	xMin = 1650
	xMax = 2900
	yMin = 0
	yMax = 1250
	sub_image = image[yMin:yMax,xMin:xMax]
	#
	#Conversión Pixel-Coordenada
	def deg2dec(g,m,s):
	"""
	función de conversión de formato gms
	a decimal
	"""
	return g+m/60.+s/(3600.)
	minLat = deg2dec(20,37,30)
	maxLat = deg2dec(20,45,00)
	minLon = deg2dec(89,53,20)
	maxLon = deg2dec(89,46,40)
	bx = minLon
	ax = (maxLon - bx)/X
	by = minLat
	ay = (maxLat - by)/Y
	def Lon(pixel):
	return ax*pixel+bx
	def Lat(pixel):
	return ay*pixel+by
	#
	#Grafica
	plt.imshow(sub_image, cmap = 'gray',extent=([Lon(xMin),Lon(xMax),Lat(yMin),Lat(yMax)]))
	plt.xlabel("Longitud")
	plt.ylabel("Latitud")
	plt.title("Jojutla (1995)")
	plt.gca().xaxis.set_major_formatter(StrMethodFormatter('{x:,.2f}')) # 2 decimales en x
	plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.2f}')) # 2 decimales en y
	plt.grid(True)
	plt.show()

view raw Lectura_BIL.py hosted with ❤ by GitHub

La visualización de la imagen resultante es esta: