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

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

For others looking to download this data without an account, it is also available via the NASA LAADS archive at https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/ — Heather QC, Jan 17 '18 at 19:19
Can you be more specific about your use case? Are you needing to interpolate 1km MODIS geolocation to 500m for use with this data? Or are you creating a 500m global grid? Also, Is there a specific reason why you are using Collection 5 data (.005. in your filename) over Collection 6? — Heather QC, Jan 17 '18 at 19:21
Thanks for the site. The article I'm following uses version 5 of MODIS. I would like to make 1km lat/lon to get the lat and lon from the station to check if the data is right and then do a validation. — Lucas Fagundes, Jan 17 '18 at 19:47
Thanks. So, am I understanding correctly that you would like to get the lat/lon bounding box from the HDF grid attributes and then use those to create a latitude and longitude array of the same shape as the data? — Heather QC, Jan 17 '18 at 22:13
You could consider trying [tag:satpy] for easier reading and processing. — gerrit, Jan 14 '22 at 08:53

score 1 · Answer 1 · answered Jan 18 '18 at 16:24

1

The latitude and longitude bounds are stored as part of a massive string in the file metadata. You can access the full string similarly to this, specifically calling the ArchiveMetadata.0 file attribute:

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

You'll need to parse the resulting string to get the values of NORTHBOUNDINGCOORDINATE, SOUTHBOUNDINGCOORDINATE, EASTBOUNDINGCOORDINATE, and WESTBOUNDINGCOORDINATE. The re module will be able to do it. Once you have those, you can use them to create 1D latitude and longitude arrays using numpy arange and if you need them to be 2D to match your data, you can create a numpy meshgrid.

answered Jan 18 '18 at 16:24

Heather QC

680
8
11

I created an array with latitude and another with longitude, but I could not join with meshgrid. I tried using .loc[] for latitude and longitude, but it does not work in array only works in df. – Lucas Fagundes Jan 19 '18 at 12:25
numpy.where is the function for finding indices based on certain conditions in numpy arrays – Heather QC Jan 19 '18 at 12:42
Yes, but I need to pick up the number for the latitude and longitude that i need. In np.where I catch the latitude and longitude separated. – Lucas Fagundes Jan 19 '18 at 13:23
Yes, and those indices give you the x and y pixel address of the data that corresponds to those lat/lon. I'm glad you got it to work in a dataframe. Good luck with your analysis! – Heather QC Jan 19 '18 at 13:28
I am not able to use np.where. What I need to do is access the latitude which is -30.275 and longitude equal -53.147. But I have an array with lat and another with lon the two with (2400,2400). And then access my array with the data and figure out what the pixel value is. But i can not. Thanks! – Lucas Fagundes Jan 19 '18 at 15:00
In [130]: lon Out[130]: 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 ] ..., [-65.26239707, -65.25695627, -65.25151547, ..., -52.22079926, -52.21535845, -52.20991765], [-65.27036446, -65.26492299, -65.25948153, ..., -52.22717449, -52.22173303, -52.21629157]]) – Lucas Fagundes Jan 19 '18 at 15:31

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

1 Answers1