I'm trying to get the vegetation index nearest to a particular lat/long. I'm getting the data from here which is 250-meter spatial resolution hdf file. The globally gridded projection is split in 36 (horizontal) by 18 (vertical) tiles, each with a resolution of approximately 10° by 10°.
Now each tile has 4800×4800 values. I can extract the GRINGPOINT
(GRINGPOINTLATITUDE/GRINGPOINTLONGITUDE) from the metadata but I'm not sure how to store values of the VI(vegetation index) corresponding to the lat/long inside the tile.
Let us suppose the GRINGPOINT is like this:
GRINGPOINTLATITUDE = [-0.0068357003079556, 9.99897831672069, 9.9909309627606, 5.67994760615426e-06]
GRINGPOINTLONGITUDE = [-179.999951582871, 179.928473473918, -169.920147289013, -169.99173290556]
Which value should I take in place of latitude and longitude in the below code ?
latitude = -90.0
longitude = -180.0
count = 0
res = []
for i in range(0,4800):
del res[:]
for j in range(0,4800):
data = {"lat":latitude,"lng":longitude,"ndvi":values[i,j]}
res.append(data)
longitude += 0.00225
collection.insert(res)
latitude +=0.00225
longitude = -180
If there is any other way of doing this then please let me know.
Best Answer
Projection of this kind of files is sinusoidal. For this one:
ftp://ladsweb.nascom.nasa.gov/allData/6/MOD13Q1/2016/129/MOD13Q1.A2016129.h07v06.006.2016147112419.hdf
the next code can access to coordinates for 256 values for subDatasets[0][0] (NDVI values).
After running the code at the Python Console of QGIS I got:
When I load the file at the Map Canvas of QGIS it looks as (dataset is part of Mexico; world_borders shapefile acts as context):
You can save the raster as *.tif and with a long/lat projection; after correcting it (see metadata to find out correction factor) to get true NDVI values.
Editing note: If your area has more than one tile is preferable that you before merge all NDVI subsets by using, for example, QGIS (Raster -> Miscellaneous -> Merge). Afterward, it can be done raster reprojection to EPSG:4326 (long/lat WGS 84). In my example, I used (with QGIS) this gdalwarp command (by default, no data are -3000):
I got this raster:
You can observe, at the cursor position, that the status bar is already showing coordinates in degrees. Finally, to get NDVI values, by using Map Algebra, correction factor is 0.0001 (and, I think, that no data it would be -0.0003).