How i create lat and lon in my MODIS file with python3

Question

I have some accounts to do for a certain latitude and longitude, but in MODIS I have neither latitude nor longitude. I roamed the internet to create latitude and longitude and then after that calculate what is needed in my work. I was told that the best program to do this kind of thing is python, and I'm using python3.

I downloaded the data on this site (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mod09a1), more specifically this one (https://e4ftl01.cr.usgs.gov/MOLT/MOD09A1.005/2017.03.22/MOD09A1.A2017081.h13v12.005.2017090201940.hdf), account is required.

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}

UPDATE!

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[+-]?\d+\.\d+)
                                  ,
                                  (?P[+-]?\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[+-]?\d+\.\d+)
                                  ,
                                  (?P[+-]?\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)

Yes, I was able to create an array of latitude and another array of longitude, but I can not join with the data to get a latitude and longitude of the pixel that interests me.

UPDATE!

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]])

How i create lat and lon in my MODIS file with python3

Answers (1)

Related Questions