Extraindo valores em latitude e longitude específica dos dados da faixa MODIS

9

Estou tentando determinar a quantidade de vapor de água precipitável (PWV), ozônio e aerossóis em função do tempo em pontos específicos da Terra, a saber, nossos observatórios astronômicos. Para fazer isso, eu já tenho um código Python usando o modapsclientqual fará o download dos produtos MODIS Aqua e Terra MYDATML2 e MODATML2 duas vezes por dia que cobrem a latitude e longitude específicas em que estou interessado.

O que não tenho certeza é como extrair as quantidades específicas que quero, como o tempo em que os dados do MODIS foram obtidos e o PWV para a posição específica de latitude e longitude do meu observatório para transformá-los em uma série temporal de valores. Os produtos MYDATML2 parecem conter 2D de latitude e longitude grades de Cell_Along_Swath_5kme Cell_Across_Swath_5kmentão eu acho que isso faz com que faixa de dados em oposição a azulejos ou grade de dados? As quantidades que eu quero, como Precipitable_Water_Infrared_ClearSkytambém parecem ser contra a Cell_Along_Swath_5kme Cell_Across_Swath_5kmbutI'm não sei como obter esse valor VOP na específica lat, long estou interessado. Ajuda por favor?

astrosnapper
fonte
Você pode fornecer um link para as imagens ou uma amostra delas?
Andrea Massetti
Claro, aqui está um exemplo de arquivo no arquivo MODIS: ladsweb.modaps.eosdis.nasa.gov/archive/allData/61/MODATML2/2018/…
astrosnapper
Olá, você teve a chance de tentar minha resposta?
Andrea Massetti 23/10
11
Desculpe, participei de uma conferência apresentando trabalhos baseados em determinações semelhantes de PWV a partir de dados sat ... 'off by 1' diferença no início do índice de matriz)
astrosnapper

Respostas:

1

[EDIT 1 - Alterei a pesquisa por coordenadas de pixel]

Usando esta amostra de MODATML que você forneceu e usando a biblioteca gdal. Vamos abrir o hdf com gdal:

import gdal
dataset = gdal.Open(r"E:\modis\MODATML2.A2018182.0800.061.2018182195418.hdf")

Em seguida, queremos ver como os subdadosets são nomeados para importar corretamente os que precisamos:

datasets_meta = dataset.GetMetadata("SUBDATASETS")

Isso retorna um dicionário:

datasets_meta
>>>{'SUBDATASET_1_NAME': 'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Optical_Thickness', 
'SUBDATASET_1_DESC': '[406x271] Cloud_Optical_Thickness atml2 (16-bit integer)',
'SUBDATASET_2_NAME':'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Effective_Radius',
'SUBDATASET_2_DESC': '[406x271] Cloud_Effective_Radius atml2 (16-bit integer)',
[....]}

Digamos que queremos obter a primeira variável, a espessura óptica da nuvem, podemos acessar seu nome por:

datasets_meta['SUBDATASET_1_NAME']
>>>'HDF4_EOS:EOS_SWATH:"E:\\modis\\MODATML2.A2018182.0800.061.2018182195418.hdf":atml2:Cloud_Optical_Thickness'

Agora podemos carregar a variável na memória chamando novamente o método .Open ():

Cloud_opt_th = gdal.Open(datasets_meta['SUBDATASET_1_NAME'])

Por exemplo, você pode acessar Precipitable_Water_Infrared_ClearSky de seu interesse fornecendo 'SUBDATASET_20_NAME'. Basta dar uma olhada no dicionário datasets_meta.

No entanto, a variável extraída não possui uma geoprojeção (var.GetGeoprojection ()), como seria de esperar de outros tipos de arquivo, como GeoTiff. Você pode carregar a variável como uma matriz numpy e plotar a variável 2d sem projeção:

Cloud_opt_th_array = Cloud_opt_th.ReadAsArray()
import matplotlib.pyplot as plt
plt.imshow(Cloud_opt_th_array)

Agora, como não há geoprojeção, examinaremos os metadados da variável:

Cloud_opt_th_meta = Cloud_opt_th.GetMetadata()

Este é outro dicionário que inclui todas as informações necessárias, incluindo uma descrição longa da subamostragem (notei que isso é fornecido apenas com o primeiro subdataset), que inclui a explicação desses Cell_Along_Swath:

Cloud_opt_th_meta['1_km_to_5_km_subsampling_description']
>>>'Each value in this dataset does not represent an average of properties over a 5 x 5 km grid box, but rather a single sample from within each 5 km box. Normally, pixels in across-track rows 4 and 9 (counting in the direction of increasing scan number) out of every set of 10 rows are used for subsampling the 1 km retrievals to a 5 km resolution. If the array contents are determined to be all fill values after selecting the default pixel subset (e.g., from failed detectors), a different pair of pixel rows is used to perform the subsampling. Note that 5 km data sets are centered on rows 3 and 8; the default sampling choice of 4 and 9 is for better data quality and avoidance of dead detectors on Aqua. The row pair used for the 1 km sample is always given by the first number and last digit of the second number of the attribute Cell_Along_Swath_Sampling. The attribute Cell_Across_Swath_Sampling indicates that columns 3 and 8 are used, as they always are, for across-track sampling. Again these values are to be interpreted counting in the direction of the scan, from 1 through 10 inclusively. For example, if the value of attribute Cell_Along_Swath_Sampling is 3, 2028, 5, then the third and eighth pixel rows were used for subsampling. A value of 4, 2029, 5 indicates that the default fourth and ninth rows pair was used.'

Eu acho que isso significa que, com base nesses pixels de 1 km, os 5 km foram construídos levando exatamente os valores de pixels em uma determinada posição na matriz de detecção 5x5 (a posição é indicada nos metadados, acho que isso é um instrumento para reduzir falhas).

De qualquer forma, neste ponto, temos uma matriz de células 1x1 km (veja a descrição da subamostragem acima, não tenho certeza sobre a ciência por trás dela). Para obter as coordenadas de cada centróide de pixel, precisamos carregar os subdatasets de latitude e longitude.

Latitude = gdal.Open(datasets_meta['SUBDATASET_66_NAME']).ReadAsArray()
Longitude = gdal.Open(datasets_meta['SUBDATASET_67_NAME']).ReadAsArray()

Por exemplo,

Longitude
>>> array([[-133.92064, -134.1386 , -134.3485 , ..., -154.79303, -154.9963 ,
    -155.20723],
   [-133.9295 , -134.14743, -134.3573 , ..., -154.8107 , -155.01431,
    -155.2256 ],
   [-133.93665, -134.1547 , -134.36465, ..., -154.81773, -155.02109,
    -155.23212],
   ...,
   [-136.54477, -136.80055, -137.04684, ..., -160.59378, -160.82101,
    -161.05663],
   [-136.54944, -136.80536, -137.05179, ..., -160.59897, -160.8257 ,
    -161.06076],
   [-136.55438, -136.81052, -137.05714, ..., -160.6279 , -160.85527,
    -161.09099]], dtype=float32)        

Você pode perceber que as coordenadas de latitude e longitude são diferentes para cada pixel.

Digamos que seu observatório esteja localizado nas coordenadas lat_obs, long_obs, do que você minimiza a diferença de coordenadas x, y:

coordinates = np.unravel_index((np.abs(Latitude - lat_obs) + np.abs(Longitude - long_obs)).argmin(), Latitude.shape)

e extraia seu valor

Cloud_opt_th_array[coordinates]
Andrea Massetti
fonte
Obrigado pela informação, mas estou tendo problemas com a parte de conversão coordenada; as Longitude_pxe Latitude_pxsão matrizes de comprimento zero. Também existe uma maneira de lidar com a conversão usando- gdalse? (. invés de depender de uma aproximação de 1 grau é X nº de milhas e, em seguida, re-aproximando-se ao km)
astrosnapper
Latitude e Longitude são fornecidas como subdatasets, ou seja, 66 e 67. Atualizarei a segunda parte.
Andrea Massetti
@astrosnapper agora a resposta deve responder completamente à sua pergunta.
Andrea Massetti