Geometric mismatch using cgImage and cgDrawShapes [message #92042] |
Fri, 02 October 2015 14:10  |
Jonas Ardo
Messages: 12 Registered: April 2014
|
Junior Member |
|
|
Hi
I have shape file (covering Senegal) in Sinusoidal projection I like to overlay on an image with the same projection.
Plotting the shape file only using cgDrawShapes works fine (http://postimg.org/image/8ax6pgm3x/) but when first using
cgImage to draw an image with the shape file on top, the data sets don't line up correctly and there is a geometric
mismatch ( http://postimg.org/image/s796ijbjj/).
The image and the shape file fit perfectly when displayed using ENVI.
What have I missed in order to get a good geometric match between the shape file and the image?
Code used:
pro mapmap
compile_opt idl2
e = ENVI(/HEADLESS)
file1 = 'D:\forskning\2015\MOD09A1\XPP\Senegal\MEAN_NPP_SENEGAL.img'
raster1 = e.OpenRaster(file1)
meannpp = raster1.GetData()
cgWindow, WAspect=1
shapefile = 'D:\forskning\2015\MBOW_etal\Brandt\data\Senegal_long_rot.sh p'
mapCoord = Obj_New('cgMap', 'SINUSOIDAL', LIMIT=[12, -18, 17, -11], $
Position=[0.05, 0.05, 0.95, 0.95], CENTER_LATITUDE=0, CENTER_LONGITUDE=0, $
FALSE_EASTING=0, FALSE_NORTHING = 0, SPHERE_RADIUS=6371007.181 )
mapCoord -> AddCmd
cgloadct,53, /Window, /Addcmd
cgImage, meannpp, /KEEP_ASPECT_RATIO, /ORDER, /WINDOW,MINVALUE=1,maxvalue=2000, /addcmd
cgDrawShapes, shapefile, Colors='black', MapCoord=mapCoord, /projected_XY, /AddCmd
cgMap_Grid, LatDel = 1.0, LonDel = 1.0, /Box_Axes, Color=cgColor('blue'), $
Map_Structure=mapCoord, /AddCmd
END
--
Jonas Ardö
Lund University, Sweden
|
|
|
Re: Geometric mismatch using cgImage and cgDrawShapes [message #92043 is a reply to message #92042] |
Fri, 02 October 2015 14:47   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Jonas Ardo writes:
> I have shape file (covering Senegal) in Sinusoidal projection I like to overlay on an image with the same projection.
> Plotting the shape file only using cgDrawShapes works fine (http://postimg.org/image/8ax6pgm3x/) but when first using
> cgImage to draw an image with the shape file on top, the data sets don't line up correctly and there is a geometric
> mismatch ( http://postimg.org/image/s796ijbjj/).
> The image and the shape file fit perfectly when displayed using ENVI.
>
> What have I missed in order to get a good geometric match between the shape file and the image?
A couple of things, I think. First, it doesn't look to me like your
image and your map have the same aspect ratio. So, because you keep the
aspect ratio of your image, it doesn't fit properly into the map
projection. Second (and possibly related to the first), you use the
POSITION keyword on your map, but neglect to use it for your image. I
would remove the KEEP_ASPECT keyword and use the POSITION keyword
itself. Or, more likely, I would create a window with the aspect ratio
of the image, since you say it has been created with a map projection.
Also, don't do this:
Color=cgColor('blue')
on Coyote Graphics commands. Just do this:
Color='blue'
Only use cgColor if you are trying to get "normal" IDL commands to
behave properly. This will avoid problems for people who are still using
the indexed color model.
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: Geometric mismatch using cgImage and cgDrawShapes [message #92048 is a reply to message #92043] |
Sat, 03 October 2015 12:10   |
Jonas Ardo
Messages: 12 Registered: April 2014
|
Junior Member |
|
|
Thanks for suggestions (and for providing cg*).
Removing /KEEP_ASPECT_RATIO and Position don't make a big difference (http://postimg.org/image/7gthbxqfh/).
Adding
cgMap_Continents, /ADDCMD, MAP_STRUCTURE=mapCoord, Color='red', /hires, LINESTYLE=3, THICK=2, /COUNTRIES, /COASTS
show that my own shape file line up with the vectors drawn by cgMap_Continents but the image is still geometrically off.
The datum of the image use is from MODIS (a sphere with radius = 6371007.181) and can perhaps be the source of the
error? Don't know if this is a standard supported datum or if there are any IDL-ENVI differences in handling geometry?
Regards
/Jonas
Code used:
pro mapmap
compile_opt idl2
e = ENVI(/HEADLESS)
file1 = 'D:\forskning\2015\MOD09A1\XPP\Senegal\MEAN_NPP_SENEGAL.img'
raster1 = e.OpenRaster(file1)
meannpp = raster1.GetData()
cgWindow
shapefile = 'D:\forskning\2015\MBOW_etal\Brandt\data\Senegal_long_rot.sh p'
mapCoord = Obj_New('cgMap', 'SINUSOIDAL', LIMIT=[11, -18, 18, -11], $
CENTER_LATITUDE=0, CENTER_LONGITUDE=0, $
FALSE_EASTING=0, FALSE_NORTHING = 0, SPHERE_RADIUS=6371007.181 )
mapCoord -> AddCmd
cgloadct,53, /Window, /Addcmd
cgImage, meannpp, /ORDER, /WINDOW,MINVALUE=1,maxvalue=2000, MapCoord=mapCoord, /addcmd
cgDrawShapes, shapefile, Colors='black', MapCoord=mapCoord, /projected_XY, /AddCmd
cgMap_Continents, /ADDCMD, MAP_STRUCTURE=mapCoord, Color='red', /hires, LINESTYLE=3, THICK=2, /COUNTRIES, /COASTS
cgMap_Grid, LatDel = 1.0, LonDel = 1.0, /Box_Axes, Color='blue', $
Map_Structure=mapCoord, /AddCmd
END
On 02/10/2015 23:47, David Fanning wrote:
> Jonas Ardo writes:
>
>> I have shape file (covering Senegal) in Sinusoidal projection I like to overlay on an image with the same projection.
>> Plotting the shape file only using cgDrawShapes works fine (http://postimg.org/image/8ax6pgm3x/) but when first using
>> cgImage to draw an image with the shape file on top, the data sets don't line up correctly and there is a geometric
>> mismatch ( http://postimg.org/image/s796ijbjj/).
>> The image and the shape file fit perfectly when displayed using ENVI.
>>
>> What have I missed in order to get a good geometric match between the shape file and the image?
>
> A couple of things, I think. First, it doesn't look to me like your
> image and your map have the same aspect ratio. So, because you keep the
> aspect ratio of your image, it doesn't fit properly into the map
> projection. Second (and possibly related to the first), you use the
> POSITION keyword on your map, but neglect to use it for your image. I
> would remove the KEEP_ASPECT keyword and use the POSITION keyword
> itself. Or, more likely, I would create a window with the aspect ratio
> of the image, since you say it has been created with a map projection.
>
> Also, don't do this:
>
> Color=cgColor('blue')
>
> on Coyote Graphics commands. Just do this:
>
> Color='blue'
>
> Only use cgColor if you are trying to get "normal" IDL commands to
> behave properly. This will avoid problems for people who are still using
> the indexed color model.
>
> Cheers,
>
> David
>
>
--
Jonas Ardö
Lund University, Sweden
|
|
|
Re: Geometric mismatch using cgImage and cgDrawShapes [message #92050 is a reply to message #92048] |
Sat, 03 October 2015 18:24   |
David Fanning
Messages: 11724 Registered: August 2001
|
Senior Member |
|
|
Jonas Ardo writes:
> Removing /KEEP_ASPECT_RATIO and Position don't make a big difference (http://postimg.org/image/7gthbxqfh/).
Well, again, you have to use the SAME position keyword on BOTH the
Map_Coord command AND the cgImage command. The two graphics are going
into different locations in the window, clearly.
mapCoord = Obj_New('cgmap', ..., POSITION=thisPosition)
cgImage, image, POSITION=thisPosition
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: Geometric mismatch using cgImage and cgDrawShapes [message #92083 is a reply to message #92050] |
Wed, 07 October 2015 10:51   |
Jonas Ardo
Messages: 12 Registered: April 2014
|
Junior Member |
|
|
OK, thanks, got it working OK (http://postimg.org/image/pwpsaj6n3/) using
http://www.idlcoyote.com/map_tips/modis_overlay.html
But, when looking closely there is an extra horizontal line in the colorbar (at about1100-1200). How can I avoid that?
Cheers
/Jonas
On 04/10/2015 03:24, David Fanning wrote:
> Jonas Ardo writes:
>
>> Removing /KEEP_ASPECT_RATIO and Position don't make a big difference (http://postimg.org/image/7gthbxqfh/).
>
> Well, again, you have to use the SAME position keyword on BOTH the
> Map_Coord command AND the cgImage command. The two graphics are going
> into different locations in the window, clearly.
>
> mapCoord = Obj_New('cgmap', ..., POSITION=thisPosition)
> cgImage, image, POSITION=thisPosition
>
> Cheers,
>
> David
>
--
Jonas Ardö
Lund University, Sweden
|
|
|
|