Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
Mercator Map Projection Image and Contour "Bleeding" Problem
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