Procesando un shapefile de unidades litológicas.

Esta es la parte 1 de la sección: Integración y exploración. 

En esta serie de post-tutoriales, se integrara y procesará información geológica, geofísica y geoquímica, como normalmente se haría en un proyecto de exploración. Apoyándonos en las ventajas y bondades que nos ofrece el potente lenguaje Python. 

Todos los datos preprocesados se usaran en publicaciones posteriores para tratar temas como geoestadística, modelado, ciencia de datos y aprendizaje automático. 

En esta primera parte decidí procesar datos cartográficos en formato shp usando las librerías GeoPandas y Folium para usar a python como un SIG.

Ubicación. 

Salt Lake, Utha. EUA.

Me decidí por esta zona porque al SW de la ciudad, se encuentra una de las minas más grandes de cobre en Estados Unidos (Mina Bingham Canyon).

Datos.

Todos los datos fueron descargados del USGS-Mineral Resources Program. Los datos se componen de:

• Un set de datos aeromagnéticos (Anomalías magnéticas y radiométricas en formato .xyz)
• Una carta geológica de Utha en formato .shp.

Geología.

Corresponde a una base de datos de unidades geológicas de Utha, con litología, edad y datos de estructuras. 

 

Geofísica.

Titulo: Airbone geophysical survey: Salt Lake City 1° y 2° Quadrangle.

Publicación: 2009.

Forma de presentación de los datos geoespaciales: Tabular digital data.

Espaciado de linea de vuelo: 3 millas.

Altura: 400 ft

Geoquímica del suelo. 

Los datos los descargue del USGS del programa National Geochemical Database-Soil. Cuyo formato es CSV.

Delimitar locación.

Usaremos la geofísica para obtener los limites del área de estudio. 

Lo primero que tenemos que hacer, es exportar todas las librerías que utilizaremos en este trabajo de integración

 

#importar librerias 
import pyproj
import numpy as np
import pandas as pd
import folium
import geopandas as gpd
import geojson
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.colors import to_hex
from shapely import geometry 
from geopandas import geoseries
from pyproj import Transformer 
   
 

De los metadatos de la geofísica (salt_lake_city_meta.txt) copie los encabezados  y los ingrese como un objeto lista, para poder integrarlos al dataframe que haremos con el archivo .xyz, que es donde se encuentran nuestros datos magnéticos y las coordenadas. 

Con el encabezado podemos identificar las columnas que contienen la latitud y la longitud.

#definimos los nombres de las columnas.
nombre_columnas=('line','find','time','day','year','latitude','longitude','radalt','totmag','resmag','diurnal','geology','resmagCM4') 
   
 

Ahora, podemos exportar nuestros datos a un pandas (Dataframe), si queremos utilizar solo algunas columnas, podemos usar el método  usecols . Como nuestro archivo está de limitado por espacios usaremos el argumento  delim_whitespace=True .

  
 
   
#se exportan los datos geofísicos .xyz
magne_data=pd.read_csv('/content/drive/My Drive/paracolab/mapas_Utha/Geofisica/salt_lake_city_mag.xyz', delim_whitespace=True,names=nombre_columnas,usecols=['longitude','latitude','totmag'])
#visualizar los 5 primeros elementos del dataframe
magne_data.head()
   
 

De los metadatos, sabemos que las coordenadas de los datos magnéticos están referenciadas al datum NAD27;  con la librería pyproj transformaremos el datum a WGS84 y UTM Zona 12.

  
#Se transforman las coordenadas de NAD27 a WGS84 Zona 12
transform1=Transformer.from_crs('epsg:4267','epsg:4326')

transform2=Transformer.from_crs('epsg:4267','epsg:32612',always_xy=True)

#creamos cuatro columnas con las conversiones

magne_data['long_WGS84'],magne_data['lat_WGS84']= transform1.transform(magne_data.longitude.values,magne_data.latitude.values)

magne_data['E_utm'],magne_data['N_utm']=transform2.transform(magne_data.longitude.values,magne_data.latitude.values)
  

Gepandas es una librería muy versátil para el manejo de datos geoespaciales. El objeto que nos entrega es un Geodataframe, que se puede manipular como si fuera un dataframe de pandas. 

"GeoPandas es un proyecto de código abierto para facilitar el trabajo con datos geoespaciales en Python. GeoPandas amplía los tipos de datos utilizados por pandas para permitir operaciones espaciales en tipos geométricos. Las operaciones geométricas son realizadas por shapely . Geopandas depende además de fiona para acceder a los archivos; descartes y matplotlib para graficar."

Transformaremos nuestro dataframe a un geodataframe con geopandas, para eso crearemos una columna llamada "geometry" que contendrá las geometrías de nuestro archivo. Lo que haremos es iterar con un ciclo for a través de las filas y transformar los valores de longitud y latitud a puntos por cada entrada. Luego lo agregamos a un geodataframe.

 

#Creamos una columna con la geometria 
magne_data['geometry']=[geometry.Point(x,y) for x,y in zip(magne_data['long_WGS84'],magne_data['lat_WGS84']) ]
  

#creamos el Geodataframe con la geometria que utilizaremos (puntos)
magne_data=gpd.GeoDataFrame(magne_data,geometry='geometry',crs='EPSG:4326')
  

Guardamos nuestro Geodataframe en nuestro disco de almacenamiento, este archivo lo utilizaremos más adelante cuando trabajemos con datos geofísicos. 

#guardamos el geodataframe como un csv
magne_data.to_csv('/content/drive/My Drive/paracolab/mapas_Utha/Geofisica/Datos_magneticos.csv')

#visualisamos los primeros 5 elementos
magne_data.head()

Utilizamos el método envelope para obtener el rectángulo delimitados de todos los puntos de nuestra geometría. 

#construir un rectangulo a apartir de nuestros puntos
puntos=geometry.MultiPoint(magne_data['geometry'])
grada=puntos.envelope

Conseguimos las coordenadas de los límites y calculamos la media de las coordenadas usando numpy, para centrar nuestro rectángulo en nuestro mapa. 

#Obtenemos las coordenadas del límite y centrales para centrar nuesto poligono
coords=np.vstack(grada.boundary.coords.xy)
centro_mapa=list(coords.mean(1))[::-1]
centro_mapa
[40.400639999999996, -111.2]

Ahora visualizaremos nuestra área en un mapa creado gracias a la librería folium que nos ayuda a visualizar mapas rápidamente usando python.




"folium facilita la visualización de datos que se han manipulado en Python en un mapa de folleto interactivo. Permite tanto la vinculación de datos a un mapa para las visualizaciones como vectores / ráster / HTML en marcadores del mapa."

#creamos nuestro mapa 
m=folium.Map(location=centro_mapa, zoom_start=8, control_scale=True)
folium.PolyLine(coords[::-1].T).add_to(m)
folium.LatLngPopup().add_to(m)

Integración de la geología. 

Como ya se mencionó los archivos se descargaron de la base de datos del servicio geológico de los Estados Unidos. Este archivo corresponde a un shapefile que contiene información sobre la litología del estado de Utah.

La librería GeoPandas lee shapefiles directamente con el método read_file. 

#se carga el shp de la carta geológica y con ella se crea un geodataframe georeferenciado
geologia_shp=gpd.read_file('/content/drive/My Drive/paracolab/mapas_Utha/Geologia/utgeol_poly_dd.shp')
geologia=gpd.GeoDataFrame(geologia_shp, geometry='geometry',crs='EPSG:4326')

#podemos ingresar a cualquier dato de nuestra tabla usando el metodo .loc
geologia.loc[132]

Visualizamos nuestro geodataframe con el metodo plot.

#visualizamos nuestro mapa 
ax=geologia.plot(facecolor='w', figsize=(12,8),  edgecolor='black')
fig=ax.figure; fig.tight_layout()

Ahora vamos a comprobar si cada polígono cruza con nuestra área de estudio. Si lo hace, vamos a recortar utilizando el método intersection de acuerdo a los limites y actualizamos la entrada. Si no se cruzan, agregamos su indice a la lista y lo colocamos después del ciclo. 

#como nuestros datos geofisicos no cubren toda el area de la carta
#vamos a recortar solo el area que nos interesa usando el rectangulo de nuestros datos geofisicos
topdrop=[]

#Se icia la iteración
for i, row in geologia.iterrows():
   
   #Se comprueba si se cruzan con los límites
    if row.geometry.intersects(grada):
           
           #Se recorta la intersección
           geologia.loc[i].at['geometry'] = row.geometry.intersection(grada)

    else:
        topdrop.append(i)
        
geologia.drop(geologia.index[topdrop],inplace=True)

Veamos nuestro nuevo geodataframe recortado. 

#visualizamos nuestra area de ineterés
ax=geologia.plot(facecolor='w',figsize=(12,8),edgecolor='black')
fig=ax.figure; fig.tight_layout()

Aplicar colores.

Para aplicar los colores a nuestro mapa, utilizaremos el sistema de códigos RGB para cartografía geológica (tabla de colores litológicos) del USGS, cuyo archivo está en formato txt.

Primero tenemos que cambiar todas las entradas que sean diferentes a la tabla de colores. 

 

#actualizar las entradas del tipo de roca para que coincidan con el txt de las claves de colores del USGS
geologia['ROCKTYPE1'] = geologia['ROCKTYPE1'].replace(['carbonate'], 'carbonate rock')
geologia['ROCKTYPE1'] = geologia['ROCKTYPE1'].replace(['clastic'], 'clastic rock')
geologia['ROCKTYPE1'] = geologia['ROCKTYPE1'].replace(['coarse-grained mixed clastic'], 'coarse-grained mixed clastic rock')
geologia['ROCKTYPE1'] = geologia['ROCKTYPE1'].replace(['eolian'], 'eolian material')

Importamos el sistema de colores en un dataframe y cambiamos los caracteres de la columna text a minúsculas. 

#cargar los codigos de colores para mapas segun el USGS
colores=pd.read_csv('/content/drive/My Drive/paracolab/mapas_Utha/lithrgb.txt',sep=' ',delimiter='\t')
colores.index=colores.text.str.lower()
colores.replace('-','') 
colores.head()

Convertimos las claves RGB a formato hexadecimal y RGBA  mediante una función lambda aplicada a lo largo de cada columna (R-G-B). 

#creamos una funcion lambda para trasformar codigo RGB a RGBA y HEX
colores['rgba']=colores.apply(lambda x:'rgba(%s,%s,%s,%s)'% (x['r'],x['g'],x['b'],255),axis=1)
colores['hex']=colores.apply(lambda x:to_hex(np.asarray([x['r'],x['g'],x['b'],255])/255.0), axis=1)
colores.head()

Como tenemos un mismo polígono para varios tipos de roca. Separamos solo el grupo que nos interesa (ROCKTYPE1) mediante el metodo .groupby. Y verificamos el nombre de las rocas en el grupo. 

#creamos un grupo donde solo se encuentren nuestro tipo de roca con el metodo .groupby
grupoderoca=geologia.groupby('ROCKTYPE1')

#podemos visualizar los tipos de roca presentes en nuestra área
print('--'.join(grupoderoca.groups))
alluvium--arenite--ash-flow tuff--basalt--biogenic sediment--carbonate rock--clay or mud--conglomerate--dacite--dolostone (dolomite)--fine-grained mixed clastic rock--glacial drift--granodiorite--landslide--limestone--medium-grained mixed clastic--metamorphic rock--mudstone--rhyolite--sandstone--shale--siltstone--water

Apliqué los colores al mapa por medio de una función lambda almacenándola en una variable llamada style_function, posteriormente creamos un archivo Geojson para exportar nuestra tabla directamente a un mapa. 

m = folium.Map(location=centro_mapa, zoom_start=8, height='100%', control_scale=True)
fog=folium.FeatureGroup(name='geologia').add_to(m)
 
for rock in grupoderoca.groups.keys():
    rock_df = geologia.loc[geologia.ROCKTYPE1 == rock]
    
    style_function=lambda x: {'color':colores.loc[x['properties']['ROCKTYPE1']]['rgba'],
                             'opacity':0,
                             'fillOpacity':.85,}
    g=folium.GeoJson(rock_df, style_function=style_function, name=rock.title())
    g.add_child(folium.Popup(rock.title()))
    fog.add_child(g)
folium.LayerControl().add_to(m)
m
   
 

¡Tenemos nuestro mapa coloreado!

Calcular el área total por grupo de roca.

La ventaja que tenemos al utilizar python y librerías como pandas y matplotlib, es que podemos manejar estadísticas de nuestros archivos y mapas con gran facilidad.

En esta sección se calculara el área que ocupa cada tipo de roca.

Primero calculé el área de cada polígono para obtener el área en metros cuadrados. Después transformé los polígonos al sistema UTM y dividí toda la columna AREA entre 10^6 para obtener el área en kilómetros cuadrados. 

#calculamos el area de cada tipo de roca en km2
geologia['AREA_KM2']=(geologia.to_crs(32612).area)/10**6
   
 

#visualizamos que todos los datos sean correctos
geologia=geologia.loc[:, ['AREA_KM2','ROCKTYPE1']]
geologia.head()
   
 

Obtenemos la suma de todos los grupos de roca y ordenamos las áreas.

grupoderoca2=geologia.groupby('ROCKTYPE1')
suma_grupoderoca=grupoderoca2.sum()
suma_grupoderoca['ROCA']=suma_grupoderoca.index
suma_grupoderoca=suma_grupoderoca.sort_values('AREA_KM2',ascending=False)
   
 

Y por ultimo utilicé seaborn para gráficar el área ocupada por cada roca en nuestra área de estudio. 

fig, ax = plt.subplots(figsize = (10, 10))
ax=sns.barplot(y='ROCA',x='AREA_KM2',data=suma_grupoderoca)
fig=ax.get_figure() ; fig.tight_layout()
   
 

Podemos ver que en nuestra área de estudio dominan rocas sedimentarias posiblemente de ambientes de transición y plataforma por la abundancia de areniscas, lutitas y conglomerados. Otro grupo importante de rocas pertenece a ambiente glacial. 

Fuentes:

Mapa geológico:

Hintze, L.F., 1980, Geologic map of Utah: Utah Geological and Mineral Survey, scale 1:500,000. 

Datos aeromagnéticos 

By: D.R. Mabey, M.D. Crittenden Jr., H.T. Morris, R.J. Roberts, and E.W. Tooker