Re: how to make lat/long grid for MODIS in IDL [message #86494 is a reply to message #86352] |
Wed, 13 November 2013 16:30   |
dm_gty88
Messages: 8 Registered: October 2013
|
Junior Member |
|
|
Thanks for the reply! I've tried to do it but to no avail. Here's the code I used:
; Opens the file
fileID = HDF_SD_START('MOD09A1.A2002081.h29v07.005.2007136054406.hdf' , /read)
index1 = HDF_SD_NAMETOINDEX(fileID, 'sur_refl_b01')
dataset1 = HDF_SD_SELECT(fileID, index1)
HDF_SD_GETDATA, dataset1, data1
HDF_SD_ENDACCESS, dataset1
; Reverses the file to make it upright
data1r = reverse(data1, 2)
; Setup map projection
smap = MAP_PROJ_INIT('Sinusoidal', LIMIT=[10,111.6969,20,127.7102])
; Forward transform
fmap = MAP_PROJ_FORWARD(data1r, MAP_STRUCTURE=smap)
When I print, fmap all I get is an array of zeroes.
On Wednesday, October 30, 2013 8:53:19 PM UTC+8, David Fanning wrote:
> dm_gty88 writes:
>
>
>
>> Just for an example, using the file in http://www.idlcoyote.com/map_tips/warptomap.php, I have only peruimage. I want to generate peru_lat and peru_lon. How would I do that?
>
>
>
> If you just have an image, you are hosed. But, suppose you have an image
>
> and the lat/lon of the four corner pixels and you know the map
>
> projection the image is in. Then, you are golden!
>
>
>
> Imagine an image printed on a piece of paper. Then, imagine you have a
>
> piece of screen left over from when you repaired the front bedroom
>
> window. By an unbelievable coincidence, the grid of the screen is just
>
> exactly the size of one image pixel.
>
>
>
> Lay the screen down over the image, and overlay the screen grid so that
>
> the edges of the grid are parallel to the sides of the image. Now,
>
> rotate both the paper with the image on it, and the screen that is
>
> aligned to the image so that the sides of the image are vertical from
>
> *your* perspective. Move them both together, don't change the alignment
>
> you had already established between the image paper and the screen.
>
>
>
> What you are looking at now, is a projected meter rectangular grid
>
> overlaying your image. Each grid cell in the screen is overlayed exactly
>
> on an image pixel.
>
>
>
> Now, set up your map projection with Map_Proj_Init (or, I would use
>
> cgMap, because I like to do things the easy way). Take your corner pixel
>
> lat/lon values and forward transform them into projected meter space.
>
> Take these numbers and label the paper with the image on it. Draw some
>
> axes while you are at it along the left and bottom of the image.
>
>
>
> When you are finished, take your pen and connect the four corner pixels
>
> in clockwise order. You are looking, are you not, at a rectangular box
>
> in a XY coordinate system. And, you know the values of all four corners
>
> of the box. If you remember fourth grade math at all, it should be
>
> possible to figure how how to assign [x,y] position values to each of
>
> the screen grid cells inside the rectangle, given that you know
>
> *exactly* how many of them there are. If you don't, ask the image how
>
> big it is.
>
>
>
> Each pixel now has a "location" in the XY grid. But, you want each
>
> pixel's "location" in latitude and longitude. Simply take your handy-
>
> dandy map projection object or function and *inverse* transform those XY
>
> locations back to latitude/longitude locations.
>
>
>
> Whala! Finished!
>
>
>
> Cheers,
>
>
>
> David
>
>
>
> --
>
> David Fanning, Ph.D.
>
> Fanning Software Consulting, Inc.
>
> Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
>
> Sepore ma de ni thue. ("Perhaps thou speakest truth.")
|
|
|