| how to make lat/long grid for MODIS in IDL [message #86122] |
Tue, 08 October 2013 20:45  |
dm_gty88
Messages: 8
Registered: October 2013
|
Junior Member |
|
|
Hi,
I've been trying to make a 2-dimensional array containing lat and long values for the MOD09A1 (basically georeferencing). I've tried using the map_proj_init and map_proj_inverse functions but the output is way off the bounds of my tile. Does anyone have a way to do this in IDL? Other software/methods would be okay but it'll be great if I could use IDL for this. The only data I could supply is are the corners of the tile and other metadata in the MODIS file, like sphere radius. Thanks!
|
|
|
|
|
|
|
|
|
|
| Re: how to make lat/long grid for MODIS in IDL [message #86351 is a reply to message #86350] |
Wed, 30 October 2013 00:48  |
dm_gty88
Messages: 8
Registered: October 2013
|
Junior Member |
|
|
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?
On Wednesday, October 30, 2013 1:52:58 PM UTC+8, dm_gty88 wrote:
> Thanks for the link but I'm trying to make an array that would contain the latitude values and another to hold the longitude values. Somewhat like given a MODIS 1200x1200 image and its the coordinates for its corners, I will produce two variables containing 1200x1200 arrays with the lat and long values. Can't seem to get a result like that from the procedure you posted.
>
>
>
> On Thursday, October 10, 2013 9:35:38 PM UTC+8, David Fanning wrote:
>
>> dm_gty88 writes:
>
>>
>
>>
>
>>
>
>>> I've been trying to make a 2-dimensional array containing lat and long values for the MOD09A1 (basically georeferencing). I've tried using the map_proj_init and map_proj_inverse functions but the output is way off the bounds of my tile. Does anyone have a way to do this in IDL? Other software/methods would be okay but it'll be great if I could use IDL for this. The only data I could supply is are the corners of the tile and other metadata in the MODIS file, like sphere
>
>>
>
>> radius. Thanks!
>
>>
>
>>
>
>>
>
>> You might be interested in this article:
>
>>
>
>>
>
>>
>
>> http://www.idlcoyote.com/map_tips/modis_overlay.html
>
>>
>
>>
>
>>
>
>> 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.")
|
|
|
|
| Re: how to make lat/long grid for MODIS in IDL [message #86352 is a reply to message #86351] |
Wed, 30 October 2013 05:53 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
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.")
|
|
|
|
| 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.")
|
|
|
|
| Re: how to make lat/long grid for MODIS in IDL [message #86495 is a reply to message #86494] |
Wed, 13 November 2013 18:03 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
dm_gty88 writes:
> 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.
I can imagine, since the data values are not in EITHER latitude OR
longitude! Where are the latitude and longitude values of your corner
pixels? You told us you had those.
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.")
|
|
|
|
| Re: how to make lat/long grid for MODIS in IDL [message #86496 is a reply to message #86495] |
Thu, 14 November 2013 03:35 |
dm_gty88
Messages: 8
Registered: October 2013
|
Junior Member |
|
|
It's in the limits keyword. Did I use it correctly?
MinLon = 111.6969
MaxLon = 127.7102
MinLat = 10
MaxLat = 20
Or alternatively,
Upper Left = 20, 117.0596
Upper Right = 20, 127.7012
Lower Left = 10, 111.6969
Lower Right = 10, 121.8512
On Thursday, November 14, 2013 10:03:01 AM UTC+8, David Fanning wrote:
> dm_gty88 writes:
>
>
>
>> 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.
>
>
>
> I can imagine, since the data values are not in EITHER latitude OR
>
> longitude! Where are the latitude and longitude values of your corner
>
> pixels? You told us you had those.
>
>
>
> 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.")
|
|
|
|
| Re: how to make lat/long grid for MODIS in IDL [message #86497 is a reply to message #86496] |
Thu, 14 November 2013 05:55 |
David Fanning
Messages: 11724
Registered: August 2001
|
Senior Member |
|
|
dm_gty88 writes:
> It's in the limits keyword. Did I use it correctly?
No, you don't need the LIMITS keyword. Get rid of it. :-)
> MinLon = 111.6969
> MaxLon = 127.7102
> MinLat = 10
> MaxLat = 20
Writing values like this suggests these are NOT the corner pixels. Are
they? How do you know this? You need to get the image into a rectangular
box. You can do this if you know the CORNER PIXEL values of your image
in latitute and longitude. These are EXTREMELY unlikely to be the
minimum and maximum values of latitude and longitude for a sinusoidal
projection. Once you know these, here are the rest of the directions
again.
;*********************************************************** **
Now, set up your map projection with Map_Proj_Init. Take your CORNER
PIXEL lat/lon values and forward transform them into projected meter
space. Take these numbers (four lats and four lons, converted into
projected meters) 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.
;*********************************************************** **
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thue. ("Perhaps thos speakest truth.")
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
