Tengo un objetivo "simple" en mente que resultó ser bastante difícil de lograr. Tengo tres matrices bidimensionales lat, lon y data todas con dimensiones (24576, 24576). Los dos primeros son coordenadas de latitudes y longitudes de la variable data que estoy tratando de trazar sobre un mapa. Estoy leyendo todos estos datos de una combinación de archivos binarios y archivos de texto, por lo que cualquier operación de preprocesamiento es prácticamente imposible y debe realizarse en el script de Python.

Dada la dimensionalidad de la matriz, es prácticamente imposible, debido a las limitaciones de memoria, trazar directamente los datos incluso al elegir una proyección de mapa base en un área pequeña del globo. Ya lo intenté y aparece un error de memoria al intentar trazar con basemap.contourf.

Por lo tanto, necesito un subconjunto de la matriz ANTES de pasarla a la función de contorno. He intentado muchas cosas, pero nada parece funcionar. Mi idea era hacer algo como esto

lat_bnds, lon_bnds = [35, 50], [5, 20]
condition=((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & (lons > lon_bnds[0]) & (lons < lon_bnds[1])

Donde lats y lons son las matrices de coordenadas 2-D. Esto produce una matriz booleana 2-D que tiene la misma dimensión que lats (o equivalentemente lons) que luego puedo usar para enmascarar los datos originales, pero no para crear subconjuntos.

Usar la misma condición en numpy.where produce una matriz con forma de (2, 2564856) que no puedo usar. Creo que el problema aquí es que hay múltiples condiciones que deben cumplirse en cada punto de la matriz 2-D y no hay garantía de que esto conduzca a una submatriz rectangular contigua. Sin embargo, por el bien de trazar esos datos, realmente necesito subconjuntos, o encontrar otro método para elaborarlos.

¿Me estoy perdiendo algo obvio? ¿Existe alguna otra forma inteligente de trazar los datos originales sin encontrar ningún error?

Fuentes de los datos: http://nsidc.org/data/G02156

Breve script para leer los datos en el supuesto de que descargó los archivos necesarios:

import numpy as np
import pandas as pd
from mpl_toolkits.basemap import Basemap # Import the Basemap toolkit

with open('IMS1kmLats.24576x24576x1.double', 'rb') as f:
 data = np.fromfile(f, dtype='<d', count=24576*24576)
 lats = np.reshape(data, [24576, 24576], order='F')

with open('IMS1kmLons.24576x24576x1.double', 'rb') as f:
 data = np.fromfile(f, dtype='<d', count=24576*24576)
 lons = np.reshape(data, [24576, 24576], order='F')

widths=np.full((24576), 1, dtype=int).tolist()
data=np.array(pd.read_fwf('ims2017312_1km_v1.3.asc', skiprows=30, 
widths=widths, lineterminator='\n', header=None))
1
Guido Cioni 9 nov. 2017 a las 03:18

2 respuestas

La mejor respuesta

En caso de que alguien se encuentre con esta publicación mientras intenta trazar datos de resolución de 1 km de IMS con Python, tengo una solución para pobres que funcionará bien.

El error de memoria es generado por las rutinas de trazado, pero la matriz aún se puede almacenar en Python. Entonces, en lugar de crear subconjuntos usando una función where, solo intenté usar una indexación explícita como

lat_subset=lat[imin:imax, jmin:jmax]

Y luego trazar el resultado con plot.imshow () sin hacer un trazado de contorno o usar una proyección de mapa para tener una idea de cómo se veían los datos. Esto me permitió elegir el intervalo de índices de manera que mi región de interés estuviera dentro. Ahora pude hacer el trazado de contorno sin tener un error de memoria.

Tengo un pequeño repositorio con un cuaderno que muestra cómo trazar esos datos: https: //github.com/guidocioni/snow_ims/blob/mistral/plot_ims.ipynb. Tiene la ventaja de leer el archivo directamente en línea, aunque los archivos de coordenadas aún deben descargarse dado su tamaño.

0
Guido Cioni 10 nov. 2017 a las 16:58

No estoy seguro de cómo obtuvo una matriz de esa forma, pero considere que si desea un subconjunto de la matriz original y tiene sus condiciones, puede cortar de la lista original.

a = np.random.randint(10,50, size=(50,2))  # ---- generate coordinates

array([[44, 11],
       [40, 36],
       [19, 26],
       ..., 
       [33, 26],
       [42, 12],
       [15, 25]])

lats = a[:, 1]  # ---- latitude is Y-axis
lons = a[:, 0]  # ---- longitude is X-axis
lat_bnds, lon_bnds = [35, 50], [5, 20]  # ---- your desired bounds

condition =((lats > lat_bnds[0]) & (lats < lat_bnds[1])) & 
            (lons > lon_bnds[0]) & (lons < lon_bnds[1])

condition
array([False, False, False, ..., False, False, False], dtype=bool)

condition.shape  # => (50,)

a[condition]     # ---- slice the original array

array([[11, 45],
       [15, 40],
       [15, 43],
       [15, 49]])

Espero que te lleve en la dirección que quieres ir.

0
NaN 9 nov. 2017 a las 04:09