2

I have an HDF4 file whose StructMetadata.0 contains the following attributes:

UpperLeftPointMtrs = (-20015109.354000,1111950.519667)
LowerRightMtrs     = (-18903158.834333,0.000000)

These are X and Y distances in meters of the MODIS Tile for L3 Gridded product (Sinusoidal Projection). I want to extract/create the coordinates of all the pixels (240 x 240) in this tile given the pixel resolution is 5km. How can I acheive this in Python?

  • Would the coordinates of each pixel be located in the pixel center or one of its corners? – Marcelo Villa-Piñeros Jun 20 '19 at 21:19
  • In the pixel centre. @MarceloVilla – MNK Jun 20 '19 at 21:31
  • Can you post a sample? – alphabetasoup Jun 21 '19 at 0:53
2

You can use numpy to create one 2D array with the x's coordinates and another 2D array with the y's coordinates.

import numpy as np

# define origin and end point of the raster
ox, oy = UpperLeftPointMtrs
ex, ey = LowerRightMtrs

# define pixel width and height
pw = 5000 # 5km
ph = 5000 # 5km

# define columns and rows
cols = 240
rows = 240

# you can also compute the right lower corner if you didn't have that information
ex = ox + (cols * pw)
ey = oy - (rows * ph)

# create 1D arrays with the coordinates of each axis (shifted to the center)
x = np.arange(ox, ex, pw) + pw/2
y = np.arange(ey, oy, ph) - ph/2

# create the 2D coordinates arrays
xx, yy = np.meshgrid(x, y)

# flip yy so coordinates are descending and not ascending
yy = np.flip(yy)

Note: If you raster is skewed (i.e. pixels are not rectangular) you would have to adjust the coordinates to account for this.

|improve this answer|||||
  • I tried this code but it is not giving me the right LowerRight coordinates for 5KM. But it gives approximate results for 4.632 KM. Is it possible that the resolution is given as 5KM in MODIS documentation but it may actually be less than 5km? Or am I doing something wrong? (How can I adjust the coordinates for skewed raster?) – MNK Jun 20 '19 at 21:55
  • Also I need the coordinates in degrees rather than distances. How can that be done? – MNK Jun 20 '19 at 21:57
  • Do you have a link to download one sample image so I can test it myself? Furthermore, what do you mean with degrees rather than distances? Do you want geographic coordinates? These are planar coordinates. – Marcelo Villa-Piñeros Jun 20 '19 at 22:02
  • I have the HDF4 file for 1 modis tile. I do not need an image what I actually need is the lat/lon of all the 240x240 pixels in one modis tile. The upper left and lower right corners of this tile have coordinates in x and y distances (meters). I know the projection used is sinusoidal. – MNK Jun 20 '19 at 22:06
  • Can you use the netcdf4 library to get the coordinates directly? There's actually not enough information to even know if your grid is regular. – alphabetasoup Jun 21 '19 at 0:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.