2

对于某个纬度和经度，我有一些帐户要做，但在 MODIS 中，我既没有纬度也没有经度。我在互联网上漫游以创建纬度和经度，然后计算我的工作需要什么。有人告诉我，做这种事情的最佳程序是 python，我正在使用 python3。

我下载了这个网站上的数据（https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod09a1），更具体地说是这个（https://e4ftl01.cr.usgs.gov/MOLT/MOD09A1.005/ 2017.03.22/MOD09A1.A2017081.h13v12.005.2017090201940.hdf )，需要账号。

import numpy as np
import pyhdf
import matplotlib.pyplot as plt
from pyhdf.SD import SD, SDC
import satpy
import pyhdf.SD

# Read
#MOD09Q1 = './MOD09Q1.A2017081.h13v12.005.2017090201940.hdf'
MOD09A1 = './MOD09A1.A2017081.h13v12.005.2017090201940.hdf'
MOD09A1 = SD(MOD09A1, SDC.READ)
#MOD09Q1 = SD(MOD09Q1, SDC.READ)

datasets_dic = MOD09A1.datasets()

for idx,sds in enumerate(datasets_dic.keys()):
    print (idx,sds)

# function
def calibrate(h5):
    attr = h5.attributes()
    FillValue = attr['_FillValue']
    scale_factor = attr['scale_factor']
    valid_range = attr['valid_range']
    data = h5.get()
    data_c = data * scale_factor
    data_c[data==FillValue] = np.nan
    data_c = np.maximum(data_c, 0.0)
    data_c = np.minimum(data_c, 1.0)
    return data_c

#select
#b01 = MOD09Q1.select('sur_refl_b01')
b02 = MOD09Q1.select('sur_refl_b02')
b03 = MOD09A1.select('sur_refl_b03')
b04 = MOD09A1.select('sur_refl_b04')
b05 = MOD09A1.select('sur_refl_b05')
b06 = MOD09A1.select('sur_refl_b06')
b07 = MOD09A1.select('sur_refl_b07')
T_sup = MOD11A2.select('sur_refl_b07')
#b07 = MOD09A1.select('sur_refl_b07')
#b07 = MOD09A1.select('sur_refl_b07')

band02 = calibra(b02)

In [10]: b01.attributes()
Out[10]: 
{'calibrated_nt': 5,
 'valid_range': [-100, 16000],
 'units': 'reflectance',
 'add_offset': 0.0,
 'scale_factor_err': 0.0,
 'long_name': 'Surface_reflectance_for_band_1',
 '_FillValue': -28672,
 'add_offset_err': 0.0,
 'scale_factor': 0.0001}

更新！

import os
import re

import matplotlib as mpl
import matplotlib.pyplot as plt
#from mpl_toolkits.basemap import Basemap
#import mpl_toolkits.basemap.pyproj as pyproj
import numpy as np
from pyhdf.SD import SD, SDC
import gdal

USE_GDAL = False

    # Identify the data field.
    DATAFIELD_NAME = 'sur_refl_b01'

        GRID_NAME = 'MOD_Grid_500m_Surface_Reflectance'

        from pyhdf.SD import SD, SDC
        hdf = SD(FILE_NAME, SDC.READ)

        # Read dataset.
        data2D = hdf.select(DATAFIELD_NAME)
        data = data2D[:,:].astype(np.float64)

        # Read global attribute.
        fattrs = hdf.attributes(full=1)
        ga = fattrs["StructMetadata.0"]
        gridmeta = ga[0]

        # Construct the grid.  The needed information is in a global attribute
        # called 'StructMetadata.0'.  Use regular expressions to tease out the
        # extents of the grid. 
        ul_regex = re.compile(r'''UpperLeftPointMtrs=\(
                                  (?P<upper_left_x>[+-]?\d+\.\d+)
                                  ,
                                  (?P<upper_left_y>[+-]?\d+\.\d+)
                                  \)''', re.VERBOSE)
        match = ul_regex.search(gridmeta)
        x0 = np.float(match.group('upper_left_x')) 
        y0 = np.float(match.group('upper_left_y')) 

        lr_regex = re.compile(r'''LowerRightMtrs=\(
                                  (?P<lower_right_x>[+-]?\d+\.\d+)
                                  ,
                                  (?P<lower_right_y>[+-]?\d+\.\d+)
                                  \)''', re.VERBOSE)
        match = lr_regex.search(gridmeta)
        x1 = np.float(match.group('lower_right_x')) 
        y1 = np.float(match.group('lower_right_y')) 
        ny, nx = data.shape
        xinc = (x1 - x0) / nx
        yinc = (y1 - y0) / ny

x = np.linspace(x0, x0 + xinc*nx, nx)
y = np.linspace(y0, y0 + yinc*ny, ny)
xv, yv = np.meshgrid(x, y)

# In basemap, the sinusoidal projection is global, so we won't use it.
# Instead we'll convert the grid back to lat/lons.
sinu = pyproj.Proj("+proj=sinu +R=6371007.181 +nadgrids=@null +wktext")
wgs84 = pyproj.Proj("+init=EPSG:4326") 
lon, lat= pyproj.transform(sinu, wgs84, xv, yv)

是的，我能够创建一个纬度数组和另一个经度数组，但我无法加入数据来获得我感兴趣的像素的纬度和经度。

更新！

Out[14]: lat
array([[-30.        , -30.        , -30.        , ..., -30.        ,
        -30.        , -30.        ],
       [-30.0041684 , -30.0041684 , -30.0041684 , ..., -30.0041684 ,
        -30.0041684 , -30.0041684 ],
       [-30.0083368 , -30.0083368 , -30.0083368 , ..., -30.0083368 ,
        -30.0083368 , -30.0083368 ],
       ..., 
       [-39.99166319, -39.99166319, -39.99166319, ..., -39.99166319,
        -39.99166319, -39.99166319],
       [-39.99583159, -39.99583159, -39.99583159, ..., -39.99583159,
        -39.99583159, -39.99583159],
       [-40.        , -40.        , -40.        , ..., -40.        ,
        -40.        , -40.        ]])


Out[15]: lon
array([[-57.73502691, -57.73021365, -57.7254004 , ..., -46.19764805,
        -46.19283479, -46.18802153],
       [-57.73745225, -57.73263879, -57.72782533, ..., -46.19958872,
        -46.19477526, -46.1899618 ],
       [-57.73987809, -57.73506443, -57.73025076, ..., -46.2015298 ,
        -46.19671613, -46.19190247],
       ..., 
       [-65.26239707, -65.25695627, -65.25151547, ..., -52.22079926,
        -52.21535845, -52.20991765],
       [-65.26638035, -65.26093921, -65.25549808, ..., -52.22398654,
        -52.21854541, -52.21310428],
       [-65.27036446, -65.26492299, -65.25948153, ..., -52.22717449,
        -52.22173303, -52.21629157]])
4

1 回答 1

1

纬度和经度界限存储为文件元数据中大量字符串的一部分。您可以访问与此类似的完整字符串，特别是调用 ArchiveMetadata.0 文件属性：

all_metadata=indata.attributes()["ArchiveMetadata.0"]

您需要解析生成的字符串以获取 NORTHBOUNDINGCOORDINATE、SOUTHBOUNDINGCOORDINATE、EASTBOUNDINGCOORDINATE 和 WESTBOUNDINGCOORDINATE 的值。re模块将能够做到这一点。一旦你有了这些，你可以使用它们使用numpy arange创建一维纬度和经度数组，如果你需要它们是二维的以匹配你的数据，你可以创建一个numpy meshgrid。

于 2018-01-18T16:24:33.967 回答