Hello everyone,
What I'm trying to do is to create a contour/image and have it appear
on a map projection. The latitudes in my dataset are seperated by 5
degrees, beginning at -87.5, and ending at +87.5. Similarly, the
longitudes that I'm using are seperated 5 degrees, beginning at -180,
and ending at +180.
I first noticed the bleeding when I attempted to contour my data to a
map projection, and while the contours would appear in the northern
hemisphere, they would be overwritten in the southern hemisphere from
what appeared to be northern hemisphere equitorial data. Frustrated,
I decided to try the map_patch() routine instead, and use it to create
an image of my data and overlay it onto the projection. Strangely, I
got the same results. However, the image will display correctly for
any other map projection that I choose.
In case I haven't been clear enough describing my problem, I have
attached below a program that demonstrates it. Simply run the
program, and an image will be displayed showing my issue.
Unfortunately, I've only been able to test it on IDL Versions 5.5 and
5.6, with the same results.
ANY help, insight, advice, or experience with a similar problem would
be greatly appreciated, and I would be eternally grateful.
Thanks,
Justin Bronn
;;BEGIN TEST_MERCATOR PROCEDURE
pro test_mercator $
, MERCATOR=mercator, CYLINDRICAL=cylindrical $
, NORTH_POLAR=north_polar, SOUTH_POLAR=south_polar
;;Setting the lats array, and associated variables
lats = [-87.5, -82.5, -77.5, -72.5, -67.5, -62.5, -57.5, -52.5, $
-47.5, -42.5, -37.5, -32.0, -27.5, -22.5, -17.5, -12.5, $
-7.5, -2.5, 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, $
37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, $
82.5, 87.5]
num_lats = n_elements(lats)
min_lat = min(lats)
max_lat = max(lats)
;;Setting the lons array and associated variables
lons = [-180.0, -175.0, -170.0, -165.0, -160.0, -155.0, -150.0, $
-145.0, -140.0, -135.0, -130.0, -125.0, -120.0, -115.0, $
-110.0, -105.0, -100.0, -95.0, -90.0, -85.0, -80.0, $
-75.0, -70.0, -65.0, -60.0, -55.0, -50.0, -45.0, -40.0, $
-35.0, -30.0, -25.0, -20.0, -15.0, -10.0, -5.0, 0.0, $
5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, $
50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0, $
95.0, 100.0, 105.0, 110.0, 115.0, 120.0, 125.0, 130.0, $
135.0, 140.0, 145.0, 150.0, 155.0, 160.0, 165.0, 170.0, $
175.0, 180.0]
num_lons = n_elements(lons)
min_lon = min(lons)
max_lon = max(lons)
;;The projection type, change to CYLINDRICAL to see what I'm sort of
;;expecting
if keyword_set(mercator) then proj_type = 'MERCATOR'
if keyword_set(cylindrical) then proj_type = 'CYLINDRICAL'
if keyword_set(north_polar) then proj_type = 'NORTH_POLAR'
if keyword_set(south_polar) then proj_type = 'SOUTH_POLAR'
;;Setting the projection type if it hasn't been set
if n_elements(proj_type) eq 0 then proj_type = 'MERCATOR'
;;Using the dist() function to create a two-dimensional array
;;corresponding to the number of lats and lons
plot_data = dist(num_lons, num_lats)
;;Setting the scale minimums and maximums
scale_min = min(plot_data)
scale_max = max(plot_data)
;;Loading a color table and setting up the device
loadct, 39
device, DECOMPOSED=0, RETAIN=2
;;Creating and scaling the image
image = bytscl(plot_data, TOP=!D.TABLE_SIZE-3 $
, MIN=scale_min, MAX=scale_max) + 1
;;Setting the map projection with map_set to a mercator projection
case proj_type of
'CYLINDRICAL' : begin
map_set, 0.0, 0.0, /NOBORDER, /CONTINENTS
end
'NORTH_POLAR' : begin
map_set, 90.0, 0.0 $
, /STEREO, /NOBORDER $
, LIMIT=[37.5, -180.0, 90.0, 180]
end
'SOUTH_POLAR' : begin
map_set, -90.0, 180.0 $
, /STEREO, /NOBORDER $
, LIMIT=[-90.0, -180.0, -37.5, 180]
end
'MERCATOR' : begin
map_set, 0.0, 0.0, 90.0, /MERCATOR, CENTRAL_AZIMUTH=90.0 $
, /HORIZON, /NOBORDER, /CONTINENTS
end
else : stop, 'I do not know that projection type!'
endcase
;;Creating the image that will be warped to the map projection using
;;the map_patch() function
map_image = map_patch(image, lons, lats1 $
, XSTART=xstart, YSTART=ystart $
, XSIZE=xsize, YSIZE=ysize $
, LON0=min_lon, LON1=max_lon $
, LAT0=min_lat, LAT1=max_lat)
;;Displaying the image with IDL's TV procedure, as well as adding
;;the map outlines and grid
tv, map_image, xstart, ystart, XSIZE=xsize, YSIZE=ysize
map_continents
map_grid
end
;;END TEST_MERCATOR PROCEDURE
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
