Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
Re: Georeferencing satellite data in batch-mode ENVI?
Hi, this is not a trivial task. It begins with reading the data from the HDF file. In case it's HDF5, I have attached a smal routine, READ_HDF, which can help you with that task. Once you have extracted your data and longitude and latitude arrays, you want to display them in a chosen projection. The classical routine for this task is IMAGEMAP by Liam Gumley: http://cimss.ssec.wisc.edu/~gumley/imagemap.html However, I once wrote my own routine based on IMAGEMAP (and a lot of ideas from Rene Preusker and Chris O'Dell) which gives you an object, named, well, MAP_IMAGE, which has the same features as IMAGEMAP and some more. It adds e.g. a colour bar (using David Fannings COLORBAR routine) and lets you choose any of the IDL built-in colour tables plus some more, including one for elevation data and one with blue colours for negative values and red colours for positive values, which is nice for difference plots. You can also choose to use a discrete colour bar. It also allows you to add multiple satellite overpasses or granules subsequently into the same projection. And, above all, it can create plots with black colour onto a white background! :-) The complete file is also attached below. All keywords accepted by MAP_SET are passed through via _EXTRA, e.g. LONDEL, LATDEL, or the projection types STEREO, SINUSOIDAL, etc. The file comes with some descriptions on top and some examples at the end, so I would be pleased if you gave it a try. Besides COLORBAR, it also uses a little routine WIN, which is an extension to the built-in WINDOW command. WIN just checks whether a chosen window already exists, thus avoiding destroying and creating the same window over and over again. Once you arranged the windows on the screen, they stay where they are. The file is also attached below. Now, assuming that the data is in HDF5 and that the appropriate datasets are named "data", "longitude" and "latitude", an IDL code fragment would look like this: <IDL> data = read_hdf(hdf_file, "data") lon = read_hdf(hdf_file, "longitude") lat = read_hdf(hdf_file, "latitude") m=obj_new("map_image", data, lat, lon) obj_destroy, m </IDL> Now, this is a boring-looking black and white colour bar, so let's try <IDL> m=obj_new("map_image", data, lat, lon, /rainbow) obj_destroy, m </IDL> Which at least adds some colours. If you now set <IDL> m=obj_new("map_image", data, lat, lon, /rainbow, /sinusoidal) obj_destroy, m </IDL> you get a different projection. <IDL> m=obj_new("map_image", data, lat, lon, /rainbow, /sinusoidal, $ /horizontal) obj_destroy, m </IDL> gives you a horizontal colour bar. Now you might want to use the TITLE keyword, and REVERT_BACKGROUND for black on white, or CTABLE for a different built-in colour bar, etc. ... Best regards, Peter P.S.: After plotting the images to the screen, you might want to produce jpg or png files. For that task I'd again suggest Liam Gumley and his SAVEIMAGE routine: http://cimss.ssec.wisc.edu/~gumley/imagetools.html P.P.S.: I guess that many of the lines in the files are too long and will be cropped in this thread, so I'll try to find some webspace where I can put them. I'll let you know when that happened. ############################################################ ####### # read_hdf.pro ############################################################ ####### ;$Id: read_hdf.pro,v 1.00 2004/06/16 ; ;+ ; NAME: ; READ_HDF ; ; PURPOSE: ; Read one dataset and attributes from a hdf-file ; ; INPUT: ; filename: HDF-filename ; dataset_name: Name of the HDF dataset ; ; KEYWORDS ; attributes: If set to a named variables, it will return the attributes ; associated with the dataset ; scale: If set, the data will be scaled with the appropriate gain ; and offsets, if avaialbe in the data. Existing nodata values ; will also be scaled ; list: If set, and if dataset_name is actually the name of a valid ; dataset group, a list of all datasets and subgroups of this group ; is returned. If no datset_name is set, all datasets in the parent ; group "/" are returned. ; If LIST is set to 1 (equivalent to /LIST) the complete names ; are returned, e.g. ; "/Group1/Subgroup2/Dataset1" ; "/Group1/Subgroup2/Dataset2" ; "/Group1/Subgroup2/Dataset3". ; If LIST is set to a value greater then 1, only the lowermost level ; of the all names is returned, e.g. ; "Dataset1" ; "Dataset2" ; "Dataset3". ; verbose: Let read_hdf be a bit more talkative about what it does and ; tried to do. ; Setting verbose to 1 will also include the latest ; error message (if any), setting verbose to 2 will ; include all error messages (if any) ; ; ; PROCEDURES: ; needs h5o_is_group.pro ; ; ; MODIFICATION HISTORY: ; written by: Peter Albert, 16.06.2004 ; Added option for choosing dataset from "/" group members (PA 20050804) ; Added keyword LIST (PA 20050816) ; Forgot to close file when LIST was set. Fixed. (PA 20050822) ; Choosing a group instead of a dataset now leads to recursive calls ; to read_hdf until a valid dataset is chosen. If the keyword ; LIST is set, the list of all subgroups and datasets is returned. ; This worked previously just with "/" as base group. (PA 20051004) ; Included printing of error messages when VERBOSE is set (PA 20051005) ; ; - ; ==================================================== function read_hdf, filename, dataset_name, $ list = list, $ attributes = attributes, $ scale = scale, $ verbose = verbose ; Catch errors from e.g. non-existing datasets catch, error if error ne 0 then begin if size(file_id, /type) ne 0 then begin h5f_close, file_id if keyword_set(verbose) then begin catch, /cancel help, /Last_Message, Output = qwertz_theErrorMessage if verbose eq 1 then $ for j = 0, 1 do $ print, "[read_hdf: Error]" + qwertz_theErrorMessage[j] $ else $ for j=0,n_elements(qwertz_theErrorMessage) - 1 do $ print, "[read_hdf: Error]" + qwertz_theErrorMessage[j] endif endif return, -1 endif ; ------------------------------------------------------------ -- ; If no parameter is given and no keyword is set, ; display the list of datasets in the "root directory" ; of the HDF file. The dataset is thus "/". if n_params() eq 1 then dataset_name = "/" file_id = h5f_open(filename) ; Check whether the dataset_name given is already a dataset ; or rather a group of datasets. n_data = h5o_is_group(file_id, dataset_name) $ ? h5g_get_nmembers(file_id, dataset_name) $ : 0 ; ------------------------------------------------------------ -- ; If dataset_name is a group of datasets, show either a list of ; datasets (or subgroups) contained within or return ; the list. if n_data gt 0 then begin ; Concatenate base group name and its contents. ; Make sure that there are not too many "/"s sub_dataset_name = strarr(n_data) prefix = $ ( $ strmid(dataset_name, 0, 1) eq "/" $ ? dataset_name $ : "/" + dataset_name $ ) + $ ( $ strmid(dataset_name, strlen(dataset_name)-1, 1) eq "/" $ ? "" $ : "/" $ ) for i = 0, n_data - 1 do $ sub_dataset_name[i] = prefix + h5g_get_member_name(file_id, dataset_name, i) h5f_close, file_id ; If keyword LIST is set, just return the list of subgroups and datasets if keyword_set(list) then return, list eq 1 $ ? sub_dataset_name $ : strmid(sub_dataset_name, strpos(sub_dataset_name[0], "/", /reverse_search) + 1) $ else begin ; If keyword LIST is not set, show the list of subgroups and datasets, ; ask for selection and call read_hdf again with the chosen ; subgroup or dataset for i = 0, n_data - 1 do $ print, i, sub_dataset_name[i], format = '(i3, ": ", a)' s = "" read, s, prompt = "Read dataset number: " ; Recursively call read_hdf with new selected dataset name out = read_hdf( $ filename, $ sub_dataset_name[fix(s)], $ list = list, $ attributes = attributes, $ scale = scale, $ verbose = verbose $ ) return, out endelse endif ; ------------------------------------------------------------ -- ; Now that we have a valid datast name, we can go ahead ; and read it from the file if keyword_set(verbose) then print, "[INFO read_hdf:] Reading " + dataset_name dataset_id = h5d_open(file_id, dataset_name) dataset = h5d_read(dataset_id) ; Read the aattributes if (n_elements(attributes) ne 0) then tempvar = size(temporary(attributes)) n_att = h5a_get_num_attrs(dataset_id) if n_att gt 0 then begin for j = 0, n_att - 1 do begin att_id = h5a_open_idx(dataset_id, j) name = h5a_get_name(att_id) att_data = h5a_read(att_id) if strpos(strlowcase(name), "no_data") ne -1 then att_data = float(att_data) attributes = j ne 0 ? $ create_struct(attributes, name, att_data) : $ create_struct(name, att_data) h5a_close, att_id endfor endif h5d_close, dataset_id ; --- < Identify nodata value > --- ; --- Sometimes it's NODATA, sometimes NO_DATA --- if size(attributes, /type) eq 8 then begin tags = strlowcase(tag_names(attributes)) index = where(strpos(tags, "nodata") ne -1) if index[0] eq -1 then index = where(strpos(tags, "no_data") ne -1) attributes = create_struct(attributes, "nodata_id", index[0]) endif ; --- < /Identify nodata value > --- ; --- < Scale data > --- if keyword_set(scale) then begin gain = 1.0 offset = 0.0 if size(attributes, /type) eq 8 then begin tnames = strlowcase(tag_names(attributes)) i = where(tnames eq "gain", c) if c eq 1 then gain = attributes.(i[0]) i = where(tnames eq "scaling_factor", c) if c eq 1 then gain = attributes.(i[0]) i = where(tnames eq "intercept", c) if c eq 1 then offset = attributes.(i[0]) i = where(tnames eq "offset", c) if c eq 1 then offset = attributes.(i[0]) endif if gain ne 1 or offset ne 0 then begin if keyword_set(verbose) then print, "[INFO read_hdf:] Scaling " + dataset_name dataset = float(dataset) * gain + offset if attributes.nodata_id ne -1 then begin if keyword_set(verbose) then print, "[INFO read_hdf:] Scaling " + dataset_name + "." + tags[attributes.nodata_id] attributes.(attributes.nodata_id) = attributes.(attributes.nodata_id) * gain + offset endif endif endif ; --- < /Scale data > --- h5f_close, file_id return, dataset end ############################################################ ####### # map_image__define.pro ############################################################ ####### ;$Id: map_image__define.pro,v 1.9 2005/04/26 ; ; ;+ ; NAME: ; map_image ; ; PURPOSE: ; Colour coded display of georeferences data ; ; BASIS: ; This program is a mixture of Liam Gumley's IMAGEMAP, Rene Preusker's ; FUB_IMAGE and Chris O'Dell's WIS_IMAGE. They all use MAP_SET, ; CONVERT_COORD and and TV. What's new is that this one is object oriented, ; allowing the easy overlay of different images and the zooming into the images. ; Which doe not necessarily requires object oriented programs, but anyway. ; Uses COLORBAR (by David Fanning) to draw the colorbar. ; ; CALLING SEQUENCE: ; ; THERE ARE A LOT OF EXAMPLES AT THE END OF THE PROGRAM! ; CALL MAP_IMAGE_DEMO ; ; m = obj_new("map_image", image, lat, lon) ; obj_destroy, m ; ; OR, after compilation of map_image__define.pro ; ; map_image, image, lat, lon ; ; INPUTS: ; IMAGE: 1- or 2-dimensional data field. Or (x, y, 3) or (n, 3) ; data field for true colour images. ; LAT: ; Array or vector of latitudes ; LON: ; Array or vector of longitudes ; ; OPTIONAL KEYWORDS: ; Generally, all keywords related to MAP_SET, MAP_GRID, etc. cann be passed ; through via the _EXTRA keyword and are not all repeated here. ; Most important may be the projection keywords like STEREO, SINUSOIDAL ; and keywords like LONDEL, LATDEL, LIMIT or BOX_AXES. ; ; Keywords specific to MAP_IMAGE_DEFINE are: ; ; NO_COPY: ; Saves memory when set, but the original ; arrays are lost for the calling routine. ; ; NO_PROJECT: ; Surpresses image drawing directly when object is created. ; This is useful when you e.g. want to modify the colour table ; prior to projecting the data. ; ; NO_DRAW: ; Only surpresses drawing of the projected data. Useful when plotting ; several datasets subsequently in the same projection. ; ; NO_ERASE: ; If set, the plot region is not erased. Necessary when plotting ; several datasets subsequently in the same window. See examples below. ; ; NO_SCALE: ; If set, the original data is not bytescaled. ; ; P0LON, P0LAT, ROT: ; Must be set as keywords in order to get passed through to map_set. ; ; BOTTOM and NCOLOR: ; Are handled like in Dave Fannings COLORBAR.PRO ; ; WINDOW: ; Set plot window ; ; REVERT_BACKGROUND: ; If set, the back- and foreground colors are switched, ; typically useful for creatring black on white plots ; instead of the standard black background. ; ; CTABLE: ; Set color table ; ; RAINBOW: ; Use a predefined, colour enhanced rainbow colour table ; ; GREYSCALE: ; Use greyscale colour table. ; ; ELEVATION ; Use a 'elevation' type of color table, ranging from blue over green and ; brown to white ; Attention: Currently, only 79 colours are defined, so a call should be ; combined with ", ncolors=79" ; ; BWR ; Use a 'blue white red' colour table with negative values in blue, 0 in white ; and positve values in red ; ; DISCRETE: ; Use discrete colours. ; ; RED, GREEN, BLUE: ; Specify red, green and blue values for discrete colours. ; N.B.: The number of elements of the RED, GREEN and BLUE ; vectors must equal the number of elements of the DISCRETE ; vector minus 1 (Assume DISCRETE set to [0, 10, 20], then 2 colours ; are needed for the range 0 to 10 and 10 to 20). ; ; FLIP_COLOURS ; If set, the colour bar is reverted, eg. is ranging from red to blue instead ; of blue to red ; ; MAGNIFY: ; If plotting low resolution data, this keyword basically draws borders ; around each pixel in order to fill holes. Can be set to scalar values greater ; than zero. ; ; HORIZONTAL: ; If set, a horizontal colour bar is plotted. Default is vertical. ; ; DEBUG: ; Print some debug information. ; ;----------------------------------------------------------- -----------------------------; ; ; OBJECT METHODS: ; ; MAP_IMAGE::COLORS ; Set color table ; KEYWORDS: ; RAINBOW, GREYSCALE, CTABLE, DISCRETE ; ; MAP_IMAGE::WORLDVIEW ; Copied from the GEOV.PRO, provided by Andrey Savtchenko of the MODIS ; Data Support Team, this method shows the location of the data on a globe. ; ; MAP_IMAGE::project ; Performs the actual projection of the data on the desired grid. ; KEYWORDS: ; NO_DRAW, NO_ERASE ; ; MAP_IMAGE::display ; Displays the projected data. ; ; MAP_IMAGE::rdpix ; Executes RDPIX on the projected dataset ; ; MAP_IMAGE::colorbar ; This method calls David Fannings COLORBAR.pro. Different to COLORBAR, the ; default colour bar is vertical. You can set a keyword HORIZONTAL on order to get ; a horizontal one. ; ; MAP_IMAGE::draw_border ; Just draw a line around the draw area as this one sometimes gets lost ... ; ; MAP_IMAGE::zoom ; Call this method in order to interactively zoom into an image. ; ; MAP_IMAGE::init ; Creates a new object. ; KEYWORDS: ; NO_COPY, NO_PROJECT, DEBUG ; ; MAP_IMAGE::cleanup ; Dereference pointers etc. ; ; MAP_IMAGE::message ; Internal routine for debug information ; ; MAP_IMAGE::remove_tag ; Internal routine for the handling of structures ; ; MAP_IMAGE::set_range ; Internal routine which does a first guess for nicely looking values for the colour bar. ; ; MAP_IMAGE::set_limit ; Internal routine which sets a first-guess limit for the map projection. ; ; MAP_IMAGE::decomposed ; Sets DEVICE, decomposed = 0 for 8 bit images and decomposed = 1 for true colour images ; ; MAP_IMAGE::limit428 ; map_image expects the LIMIT keyword to MAP_SET to be an 8-element vector. This routines ; converts a given 4-element vector into an 8-element vector. ; ; MAP_IMAGE::set_one_var ; Internal method for structure handling. All keywords passed to map_image via _EXTRA ; are stored in an internal structure. ; ; MAP_IMAGE::set_var ; Internal method for structure handling. All keywords passed to map_image via _EXTRA ; are stored in an internal structure. ; ; MAP_IMAGE::get_var ; Return values of internal structure elements. ; ; MAP_IMAGE::debug ; Can be used for debug purposes in order to stop within a method (and get access to ; object variables). ;----------------------------------------------------------- -----------------------------; ; EXTERNAL SUBROUTINES ; win.pro ; colorbar.pro ; ;----------------------------------------------------------- -----------------------------; ; ; CAVEATS: ; MAP_IMAGE works on 8 bit colour mode, i.e. you should set ; DEVICE, PSEUDE=8 after you started IDL! ; ;----------------------------------------------------------- ----------------------------- ; - ; ; Modification history: ; Created by Peter Albert, 23.9.2004 ; ; 14.10.2004 PA map_image::colors now checking for keyword "discrete" after all other ; keywords -> discrete works together with "ctable" ; 01.12.2004 PA map_image::colors now sets self.var.draw.maxi if max(image) or the value ; provided via "MAXI= " is lower than max(discrete) ; It also checks ranges when the colour index (idx) is ; set. ; 13.12.2004 PA map_image::init and ::cleanup: now restoring original values ; of !p.color and !p.background ; ; The background keyword to dilate does not what I think, so ; magnifying the image plotting black on white (revert_background) ; did not work. I am now setting all pixels to 0 at the start (::project) ; and do not use the background keyword with dilate (::display). The ; remaining black pixels are set to !p.background afterwards. ; 28.12.2004 PA map_image::project: Consistent use of self.var.draw.xoffset instead of corners[0] should ; now avoid the fact that sometimes the map_set frame was erased by the image. ; Additionally, ceil and floor is now used when calculating size and offset ; from normalized device coordinates. ; ; 29.12.2004 PA colorbar / ex12: The new example ex12 shows how data with logarithmic scaling ; can be displayed. Attention: The subroutine colorbar.pro ; was modified! ; ; 05.01.2005 PA colorbar / ex12: The updated example 12 illustrates how logarithmic scaling can be achieved ; without modifications of colorbar.pro. ; map_image::project: There was a confusing offset of 2 at the numbers of used colors when ; bytescaling the data. Removed. ; ; 05.01.2005 PA/JS map_image::colors: Added colour table "maxcolors" ; ; 10.01.2005 PA map_image::colors: When keywords discrete and red are set, only set r[0] = 0 (and g, b) ; ; 11.01.2005 PA map_image::project / map_image::colors: Going back to saving 2 colours for background and void-color and 1 for ; the foreground colour. Using a maximum of 253 colours for the image. ; Avoid error messages when displaying byte data with a negative user-defined minimum value ; map_image:draw_border: Added routine ; ; 12.01.2005 PA map_image::project: Corected small bug with respect to size and offset of plot region ; map_image__define: Changed type of self.var.draw.color / background from byte to long ; map_image::project: No bytscl when displaying true colour images ; map_image::display: No setting of background colour when displaying true colour images ; ; 13.01.2005 PA Corrected some bugs with regard to true colour images, magnify and postscript output ; ; 24.01.2005 PA map_image::colors: Renamed keyword "maxcolors" to "extended_rainbow" in order to avoid ; conflicts with self.var.draw.maxi (When keyword "MAX) was set ; it was regognised for self.var.draw.maxi and for "maxcolors" ; ; 27.01.2005 PA map_image::colors: With discrete colours, the colours were sometimes defined such that ; colorbar.pro sometimes showed the steps at a slightly wrong position. ; Workaround with setting the number of colours. ; No revert_background for postscript ; ; 28.01.2005 PA map_image::set_one_var: Call set_range for new image only when mini and maxi are not yet set ; map_image::remove_tag: The new structure was built from *self.var.internal._extra instead from struct_in ; ; 01.02.2005 PA map_image::init, map_image::colorbar: When called with keyword "horizontal=0", a ; horizontal colourbar was drawn anyway. Fixed ; map_init::project: Added keyword no_scale (controlled by self.var.draw.no_scale) ; ; 11.02.2005 PA map_image::colors: Added keyword flip_colours ; ; 24.02.2005 PA map_init::colours: Added keyword "elevation" ; ; 01.03.2005 PA Added map_image::rdpix ; ; 03.03.2005 PA map_image::set_limit: Improved calculation of default limits ; for map_set ; map_image::colors: Is it a bug or is it a feature? When called for the first time ; TVLCT modifies !p.color. It will now be set back t the original ; value. ; 17.03.2005 PA map_init::colours: Added keyword "bwr" ; 15.04.2005 PA map_init::project: Check for NaNs from convert_coord before proceeding ; 26.04.2005 PA map_init::project: When drawing true colour images using /REVERT, areas where all three ; channels are equal to zero were drawn in white. Now, the minimum value ; for all three channels is set to one. ; 19.07.2005 PA map_init::project/display: Added _extra = _extra to calls to win. Allows setting XSIZE and ; YSIZE in calls to map_image ; 12.08.2005 PA map_image::colors: When keyword BWR was set together with NCOLORS and / or BOTTOM, ; white was not at the 0 value ; map_image::project: When the keyword REVERT_BACKGROUND was abbreviated to REV, it was interpreted ; as REVERSE by the first call to map_set, and the image was flipped around the ; x-axis ; 18.08.2005 PA map_image::display: deleted /horiz from call to map_grid ; 02.09.2005 PA map_image::worldview: Crashed when lon / lat data was ; 1D. Not fixed but no execution in ; this case. ; 25.10.2005 PA map_image::get_var: Added function ;----------------------------------------------------------- ----------------------------- ;----------------------------------------------------------- ----------------------------- ; ============================================================ ======================================== ; Set predefined colors ; ============================================================ ======================================== pro map_image::colors, $ rainbow = rainbow, $ extended_rainbow = extended_rainbow, $ greyscale = greyscale, $ elevation = elevation, $ bwr = bwr, $ ctable = ctable, $ discrete = discrete, $ red = red, green = green, blue = blue, $ revert_background = revert_background, $ flip_colours = flip_colours if keyword_set(rainbow) then begin self -> message, "[map_image::colors]: Rainbow color table" r = interpol([ 0, 0, 0, 0, 0, 200, 255, 255, 255, 180, 100], indgen(11), findgen(254) / 253 * 10) g = interpol([ 0, 0, 125, 255, 255, 255, 255, 125, 0, 0, 0], indgen(11), findgen(254) / 253 * 10) b = interpol([125, 255, 255, 255, 0, 0, 0, 0, 0, 0, 0], indgen(11), findgen(254) / 253 * 10) r = [0, r, 255] & g = [0, g, 255] & b = [0, b, 255] color_save = !p.color & tvlct, r, g, b & !p.color = color_save endif if keyword_set(extended_rainbow) then begin self -> message, "[map_image::colors]: Extended Rainbow color table" r = interpol([152, 217, 105, 0, 0, 7, 109, 210, 250, 250, 250, 221, 156, 100, 80], indgen(15), findgen(254) / 253 * 14) g = interpol([ 0, 0, 0, 125, 250, 245, 224, 250, 216, 159, 101, 46, 0, 0, 34], indgen(15), findgen(254) / 253 * 14) b = interpol([ 91, 230, 244, 250, 250, 140, 31, 0, 0, 0, 0, 46, 0, 0, 40], indgen(15), findgen(254) / 253 * 14) r = [0, r, 255] & g = [0, g, 255] & b = [0, b, 255] color_save = !p.color & tvlct, r, g, b & !p.color = color_save endif if keyword_set(elevation) then begin self -> message, "[map_image::colors]: 'Elevation' style: blue, green, yellow, brown, white" r = bytarr(256) g = bytarr(256) b = bytarr(256) r[2] = 0 & g[2] = 0 & b[2] = 150 r[3:20] = 0 & g[3:20] = bindgen(18) * (190.-120.) / 17. + 120 & b[3:20] = 0 r[21:40] = bindgen(20) * (220.) / 19. & g[21:40] = 190 & b[21:40] = 0 r[41:60] = 220 - bindgen(20) * (30) / 19. & g[41:60] = 190 - bindgen(20) * (80) / 19. & b[41:60] = 0 r[61:80] = bindgen(20) * 65 / 19. + 190 g[61:80] = bindgen(20) * 145 / 19. + 110 b[61:80] = bindgen(20) * 255 / 19. r[81] = 235 & g[81] = 235 & b[81] = 235 color_save = !p.color & tvlct, r, g, b & !p.color = color_save endif if keyword_set(bwr) then begin self -> message, "[map_image::colors]: Blue to red" saturation = fltarr(self.var.draw.ncolors) value = intarr(self.var.draw.ncolors) + 1 hue = intarr(self.var.draw.ncolors) white_index = round((self.var.draw.ncolors * (-1. * self.var.draw.mini ) / (self.var.draw.maxi - self.var.draw.mini))>0) hue [0 : white_index] = 240 saturation[0 : white_index] = (1. - findgen(white_index+1) / white_index) * self.var.draw.ncolors / 253. saturation[white_index + 1 : self.var.draw.ncolors - 1] = findgen(self.var.draw.ncolors-white_index-1) / (self.var.draw.ncolors-white_index-2) * self.var.draw.ncolors / 253. color_convert, hue, saturation, value, r, g, b, /hsv_rgb rr = intarr(254) & gg = intarr(254) & bb = intarr(254) rr[self.var.draw.bottom-1] = r gg[self.var.draw.bottom-1] = g bb[self.var.draw.bottom-1] = b r = [0, rr, 255] & g = [0, gg, 255] & b = [0, bb, 255] color_save = !p.color & tvlct, r, g, b & !p.color = color_save endif if keyword_set(ctable) then begin self -> message, "[map_image::colors]: Color table " + string(ctable, format = '(i2)') loadct, ctable, /silent endif if keyword_set(greyscale) then begin self -> message, "[map_image::colors]: Greyscale color table" loadct, 0, /silent endif if keyword_set(discrete) then begin self -> message, "[map_image::colors]: Discrete colors" if self.var.draw.maxi lt max(discrete) then begin self.var.draw.maxi = max(discrete) self -> message, "[map_image::colors]: Setting maxValue to " + string(self.var.draw.maxi) endif nnn = n_elements(discrete) - 1 ; Setting number of colors such that colorbar shows the steps at the correct position self.var.draw.ncolors = floor(float(self.var.draw.ncolors) / nnn) * nnn + 1 ; print, nnn, self.var.draw.ncolors ddd = bytscl([min(*self.var.data.image, max=maxi), discrete, maxi], $ min = self.var.draw.mini, $ max = self.var.draw.maxi, $ top = self.var.draw.ncolors - 1, /nan) + self.var.draw.bottom ddd = ddd[1:nnn+1] idx = intarr(256) - 1 for i = 0, nnn-1 do if ddd[i] lt ddd[i+1] then idx[ddd[i]:ddd[i+1]] = i if keyword_set(red) and keyword_set(blue) and keyword_set(green) then begin if n_elements(red) ne nnn or n_elements(green) ne nnn or n_elements(blue) ne nnn then begin print, nnn, format = '("[map_image::colors] ERROR: red, green and blue must have " , i3, " elements")' help, red, green, blue return endif endif else begin ;Pick the colours from within the chosen intervals tvlct, red, green, blue, /get index = round(findgen(nnn) / (nnn - 1.) * (self.var.draw.ncolors - 1) + self.var.draw.bottom) red = red [index] green = green[index] blue = blue [index] endelse r = red[idx>0] & g = green[idx>0] & b = blue[idx>0] bad = where(idx eq -1) if bad[0] ne -1 then begin r[bad] = self.var.draw.void_color g[bad] = self.var.draw.void_color b[bad] = self.var.draw.void_color endif color_save = !p.color & tvlct, r, g, b & !p.color = color_save endif if keyword_set(flip_colours) then begin tvlct, r, g, b, /get color_save = !p.color & tvlct, reverse(r), reverse(g), reverse(b) & !p.color = color_save endif ;Set white, black and void-color: tvlct, r, g, b, /get r[0] = 0 & g[0] = 0 & b[0] = 0 r[255] = 255 & g[255] = 255 & b[255] = 255 r[1] = self.var.draw.void_color & g[1] = self.var.draw.void_color & b[1] = self.var.draw.void_color color_save = !p.color & tvlct, r, g, b & !p.color = color_save if keyword_set(revert_background) then begin if strupcase(!d.name) eq 'PS' then begin self -> message, "[map_image::colors]:Don't revert colours when using postscript" endif else begin self -> message, "[map_image::colors]: Reverting background and plot colour:" background = !p.background color = !p.color self -> message, "[map_image::colors]: !p.background: " + $ strcompress(string(background), /rem) + " --> " + $ strcompress(string(color), /rem) self -> message, "[map_image::colors]: !p.color: " + $ strcompress(string(color), /rem) + " --> " + $ strcompress(string(background), /rem) !p.background = color !p.color = background endelse endif end ; ============================================================ ======================================== ; Debug messages ; ============================================================ ======================================== pro map_image::message, msg_string, time = time if keyword_set(time) then begin if self.var.internal.debug then print, string(systime(1)-self.var.internal.time, format = '(f5.1, " sec. : ")') + msg_string endif else begin if self.var.internal.debug then print, msg_string endelse end ; ============================================================ ======================================== ; Location on the globe ; ============================================================ ======================================== pro map_image::worldview ; This IDL subroutine was copied from geov.pro, develped at the ; Modis Data Support Team ; by Andrey Savtchenko,(asavtche@daac.gsfc.nasa.gov). tvlct, r, g, b, /get loadct,0, /silent dims=size(*self.var.data.latitude) if dims[0] eq 2 then begin lon0=(*self.var.data.longitude)[dims[1]/2,dims[2]/2] lat0=(*self.var.data.latitude)[dims[1]/2,dims[2]/2] ;latmin((*self.var.data.latitude))+(max((*self.var.data.lati tude))-min((*self.var.data.latitude)))/2. mpl=[(*self.var.data.latitude)[0,dims[2]-1],(*self.var.data. longitude)[0,dims[2]-1],$ (*self.var.data.latitude)[0,0],(*self.var.data.longitude)[0, 0],$ (*self.var.data.latitude)[dims[1]-1,0],(*self.var.data.longi tude)[dims[1]-1,0],$ (*self.var.data.latitude)[dims[1]-1,dims[2]-1],(*self.var.da ta.longitude)[dims[1]-1,dims[2]-1]] window,22,xsize=400,ysize=400,retain=2,title="granule geolocation" map_set,lat0,lon0,/ortho,/isotropic,/noborder,xmargin=0,$ ymargin=0,/horizon,/grid,color=100, londel = 15, latdel=15 map_continents,fill_continents=1,color=40 map_continents,/coasts,color=100 oplot,(*self.var.data.longitude)[*,0],(*self.var.data.latitu de)[*,0],psym=3,color=255 oplot,(*self.var.data.longitude)[*,dims[2]-1],(*self.var.dat a.latitude)[*,dims[2]-1],psym=3,color=255 oplot,(*self.var.data.longitude)[0,*],(*self.var.data.latitu de)[0,*],psym=3,color=255 oplot,(*self.var.data.longitude)[dims[1]-1,*],(*self.var.dat a.latitude)[dims[1]-1,*],psym=3,color=255 tvlct, r, g, b endif else begin print, "[map_image::worldview]: Latitude and Longitude data has to be 2D for this subroutine." endelse return end ; ============================================================ ======================================== ; Remove tag from structure ; ============================================================ ======================================== function map_image::remove_tag, struct_in, tag if size(struct_in, /type) eq 8 then begin st_names = strlowcase(tag_names(struct_in)) for i = 0, n_tags(struct_in) -1 do begin match = -1 for j = 0, n_elements(tag) - 1 do if strpos(st_names[i], strlowcase(tag[j])) ne -1 then match = 0 if match eq -1 then $ output = size(output, /type) eq 8 $ ; ? create_struct(output, st_names[i], (*self.var.internal._extra).(i)) $ ; : create_struct( st_names[i], (*self.var.internal._extra).(i)) ? create_struct(output, st_names[i], struct_in.(i)) $ : create_struct( st_names[i], struct_in.(i)) endfor endif else output = 0 return, output end ; ============================================================ ======================================== ; Project image ; ============================================================ ======================================== pro map_image::project, $ no_draw = no_draw, $ no_erase = no_erase, $ _extra = _extra self -> set_var, _extra = _extra self -> colors, _extra = _extra win, self.var.internal.window, _extra = _extra ; Define map projection if keyword no_erase is not set if not keyword_set(no_erase) or (size(*self.var.draw.draw_image))[0] eq 0 then begin ; First, see whether p0lon, p0lat or rot have user-defined values if size(*self.var.internal._extra, /type) eq 8 then begin p0lon = max(strpos(strlowcase(tag_names(*self.var.internal._extra)), 'p0lon')) eq -1 $ ? 0 $ : (*self.var.internal._extra).p0lon p0lat = max(strpos(strlowcase(tag_names(*self.var.internal._extra)), 'p0lat')) eq -1 $ ? 0 $ : (*self.var.internal._extra).p0lat rot = max(strpos(strlowcase(tag_names(*self.var.internal._extra)), 'rot')) eq -1 $ ? 0 $ : (*self.var.internal._extra).rot endif else begin p0lon = 0 & p0lat = 0 & rot = 0 endelse self -> message, /time, string(p0lon, p0lat, rot, format = '("[map_image::project]: p0lon, p0lat, rot: ", 3f6.1)') map_set,p0lat, p0lon, rot, $ _extra = self -> remove_tag(*self.var.internal._extra, ['tit', 'lond', 'latd', 'rev']) endif ; Corners of draw region: ; x, y, lower right, x, y, upper left in device coordinates ; don't use *self.var.proj.position as things can be confused using !p.multi corners = ([!d.x_size * !x.window, !d.y_size * !y.window])[[0, 2, 1, 3]] * self.var.draw.scalef self.var.draw.xoffset = ceil(corners[0]) self.var.draw.yoffset = ceil(corners[1]) ; Size of the draw region xsize = floor((!d.x_size * (!x.window[1] - !x.window[0])) * self.var.draw.scalef) ysize = floor((!d.y_size * (!y.window[1] - !y.window[0])) * self.var.draw.scalef) ; Create new image or add no data to existing image if not keyword_set(no_erase) or (size(*self.var.draw.draw_image))[0] eq 0 then begin self -> message, /time, string(xsize, ysize, format = '("[map_image::project]: Create new image (", i5, ", ", i5, ")")') ptr_free, self.var.draw.draw_image self.var.draw.draw_image = ptr_new(self.var.draw.true_color $ ? bytarr(xsize, ysize, 3) $ : bytarr(xsize, ysize)) endif else self -> message, /time, string(xsize, ysize, format = '("[map_image::project]: Add to existing image (", i5, ", ", i5, ")")') self -> message, /time, "[map_image::project]: Convert_coord ..." index = convert_coord(*self.var.data.longitude, *self.var.data.latitude, /data, /to_device) * self.var.draw.scalef self -> message, /time, "[map_image::project]: Convert_coord done" ; Check for pixels where convert_coord produces NaN: bad = where(finite(index[0, *]) eq 0 or finite(index[1, *]) eq 0, c) if c gt 0 then begin index[0, bad] = 0 index[1, bad] = 0 endif x = round(0 > (index[0, *] - self.var.draw.xoffset) < (xsize - 1)) y = round(0 > (index[1, *] - self.var.draw.yoffset) < (ysize - 1)) if c gt 0 then begin x[bad] = -1 y[bad] = -1 endif self -> message, /time, "[map_image::project]: Project image start" if self.var.draw.true_color then begin z = x*0 case (size(*self.var.data.image))[0] of 3: for i = 0, 2 do (*self.var.draw.draw_image)[x, y, z * 0 + i] = (*self.var.data.image)[*, *, i] > 1 2: for i = 0, 2 do (*self.var.draw.draw_image)[x, y, z * 0 + i] = (*self.var.data.image)[*, i] > 1 else: begin print, "[map_image::project]: I am confused about the size of the image... giving up" help, *self.var.data.image end endcase endif else begin ; When displaying discrete byte data it is sometimes useful to set self.var.draw.mini to -0.5 in order to ; have a nice colour bar. In these cases bytscl creates an error message when applied to the byte data. if self.var.draw.no_scale ne 0 then $ dum_data = *self.var.data.image $ else $ if size(*self.var.data.image, /type) eq 1 and self.var.draw.mini lt 0 then $ dum_data=bytscl(fix(*self.var.data.image), $ min = self.var.draw.mini, $ max = self.var.draw.maxi, $ top = self.var.draw.ncolors - 1, /nan) + self.var.draw.bottom $ else $ dum_data=bytscl(*self.var.data.image, $ min = self.var.draw.mini, $ max = self.var.draw.maxi, $ top = self.var.draw.ncolors - 1, /nan) + self.var.draw.bottom if (*self.var.draw.void_index)[0] ne -1 then begin self -> message, string(self.var.draw.void_color, format = '("[map_image::project]: Setting void pixels to ", i3)') dum_data[(*self.var.draw.void_index)] = 1 endif (*self.var.draw.draw_image)[x, y] = (dum_data) ptr_free, self.var.draw.draw_image_orig self.var.draw.draw_image_orig = ptr_new(fltarr((size(*self.var.draw.draw_image))[1:2])) (*self.var.draw.draw_image_orig)[x, y] = *self.var.data.image endelse self -> message, /time, "[map_image::project]: Project image done" if not keyword_set(no_draw) then self -> display, _extra = _extra ;, location = location end ; ============================================================ ======================================== ; Draw projected image ; ============================================================ ======================================== pro map_image::display, _extra = _extra self -> set_var, _extra = _extra ; Magnify image if self.var.draw.magnify ne 0 then $ structuring_element = replicate(1, 3, 3) if self.var.draw.magnify gt 0 then begin self -> message, "[map_image::project]: Magnify image" for mm = 1, fix(self.var.draw.magnify) do begin if strupcase(!d.name) eq 'PS' then begin idx = where(*self.var.draw.draw_image eq !p.background) if idx[0] ne -1 then (*self.var.draw.draw_image)[idx] = 0 endif *self.var.draw.draw_image = dilate(*self.var.draw.draw_image, structuring_element, /gray, /constrained) ;; You don't need the fill and loc stuff when the keyword "constrained" is set! ; fill = dilate( *self.var.draw.draw_image, structuring_element, /gray, background = !p.background) ; loc = strupcase(!d.name) ne 'PS' $ ; ? where( ( fill ge !p.background ) and ( *self.var.draw.draw_image eq !p.background ), count ) $ ; : where( ( fill le !p.background ) and ( *self.var.draw.draw_image eq !p.background ), count ) ; if count ge 1 then (*self.var.draw.draw_image)[loc] = fill[loc] ; print, count endfor endif if self.var.draw.magnify lt 0 then begin count = 1 while count gt 0 do begin self -> message, "[map_image::project]: Auto-Magnify image" *self.var.draw.draw_image = dilate( *self.var.draw.draw_image, structuring_element, /gray, background = (1 - self.var.draw.true_color) * !p.background, /constrained) ; fill = dilate( *self.var.draw.draw_image, structuring_element, /gray, background = !p.background) ; loc = strupcase(!d.name) ne 'PS' $ ; ? where( ( fill ge !p.background ) and ( *self.var.draw.draw_image eq !p.background ), count ) $ ; : where( ( fill le !p.background ) and ( *self.var.draw.draw_image eq !p.background ), count ) ; if count gt 0 then (*self.var.draw.draw_image)[loc] = fill[loc] endwhile endif ; The image was originally set to 0, as I am not able to let the BACKGROUND ; keyword of dilate do what I expect ... so here we set all black pixels to ; !p.background if self.var.draw.true_color and !p.color eq 0 then begin idx = where(total(*self.var.draw.draw_image, 3) eq 0) if idx[0] ne -1 then begin nnn = (size(*self.var.draw.draw_image))[5] / 3 (*self.var.draw.draw_image)idx, idx + nnn, idx + nnn * 2 = 255 endif endif else begin idx = where(*self.var.draw.draw_image eq 0) if idx[0] ne -1 then (*self.var.draw.draw_image)[idx] = !p.background endelse self.var.draw.magnify = 0 win, self.var.internal.window, _extra = _extra ss = size(*self.var.draw.draw_image) / self.var.draw.scalef tv, *self.var.draw.draw_image, $ self.var.draw.xoffset / self.var.draw.scalef, $ self.var.draw.yoffset / self.var.draw.scalef, $ xsize = ss[1], ysize = ss[2], $ true = self.var.draw.true_color * 3 map_continents, /hires, /continents, _extra = *self.var.internal._extra ;, col = !p.background map_grid, _extra = *self.var.internal._extra if not self.var.draw.no_color_bar then self -> colorbar self -> draw_border end ; ============================================================ ======================================== ; Use RDPIX ; ============================================================ ======================================== pro map_image::rdpix win, self.var.internal.window rdpix, *self.var.draw.draw_image_orig, $ self.var.draw.xoffset / self.var.draw.scalef, $ self.var.draw.yoffset / self.var.draw.scalef end ; ============================================================ ======================================== ; Draw color bar ; ============================================================ ======================================== pro map_image::colorbar, _extra = _extra self -> set_var, _extra = _extra ; self -> remove_tag: ; Create temporary _extra structure with all elements of *self.var.internal._extra but position. ; Position must not be passed to colorbar, as it controls the position of the mapped image ; The position of the colorbar is controlled by self.var.draw.cb_position. ; Now this is difficult :-) ; The default behaviour of 'colorbar' is to draw a horizontal colorbar. I want the default to ; be vertical, therefore I am checking whether a keyword like 'horizontal' is provided in ; _extra, otherwise, the keyword 'vertical' is passed to colorbar. dum = strpos(strlowcase(tag_names(*self.var.internal._extra)), 'horiz') if max(dum, index) eq 0 then vertical = (*self.var.internal._extra).(index) eq 0 else vertical = 1 colorbar, $ minrange = self.var.draw.mini, maxrange = self.var.draw.maxi, $ divisions = self.var.draw.n_lev, $ vertical = vertical, $ position = self.var.draw.cb_position, $ bot=self.var.draw.bottom, ncolors=self.var.draw.ncolors - 1, $ _extra = self -> remove_tag(*self.var.internal._extra, ['pos', 'iso']) end ; ============================================================ ======================================== ; Draw border ; ============================================================ ======================================== pro map_image::draw_border plots, !x.window[[0, 1, 1, 0, 0]], !y.window[[0, 0, 1, 1, 0]], /normal end ; ============================================================ ======================================== ; Select region of interest or zoom in and out by a choosen factor (in degree or percent) ; ============================================================ ======================================== pro map_image::zoom, _extra = _extra if strupcase(!d.name) eq 'PS' then begin print, "[map_image::zoom] ERROR: method zoom not defined for postscript mode!" return endif self -> message, "[map_image::zoom]: Zoom into image" self -> set_var, _extra = _extra first_corner = 0 second_corner = 0 selection_done = 0 device, get_graphics = old, set_graphics = 6 ;Set xor print, "Press left button to choose first corner." while not selection_done do begin cursor, xx, yy, wait = 2, /dev wait, .2 butn = !err if butn eq 1 then begin if first_corner eq 0 then begin if second_corner eq 1 then begin plots, [x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], /device plots, x0, y0, /device, psym = 4 endif second_corner = 0 x0 = xx & y0 = yy first_corner = 1 plots, x0, y0, /device, psym = 4 print, "Press second button to choose second corner." endif else begin x1 = xx & y1 = yy plots, [x0, x1, x1, x0, x0], [y0, y0, y1, y1, y0], /device second_corner = 1 & first_corner = 0 print, "Press right button to quit or left button for new selection." endelse endif if butn eq 4 and second_corner eq 1 then begin if x0 eq x1 or y0 eq y1 $ then print, "Please choose rectangular area." $ else selection_done = 1 endif endwhile device,set_graphics = old, cursor_standard=30 if x0 gt x1 then begin xx = x0 & x0 = x1 & x1 = xx endif if y0 gt y1 then begin yy = y0 & y0 = y1 & y1 = yy endif new_corners = ((convert_coord([x0, (x0+x1)/2, x1, (x0+x1)/2], [(y0+y1)/2, y1, (y0+y1)/2, y0], /device, /to_data))[0:1, *])[*] (*self.var.internal._extra).limit = new_corners[[1, 0, 3, 2, 5, 4, 7, 6]] self -> project end ; ============================================================ ======================================== ; Calculate range for color bar ; ============================================================ ======================================== pro map_image::set_range self.var.draw.mini = min(*self.var.data.image, max=maxi) self.var.draw.maxi = maxi if self.var.draw.mini eq self.var.draw.maxi then begin self.var.draw.mini = self.var.draw.mini - 1 self.var.draw.maxi = self.var.draw.maxi + 1 self.var.draw.n_lev = 2 endif else begin ;self.var.draw.maxi = percentiles((*self.var.data.image), value = [.95]) binsize = 10.^float(floor(alog10((self.var.draw.maxi-self.var.draw.min i)))) count=0 start_here: count=count+1 minout = floor(self.var.draw.mini / binsize) * binsize maxout = ceil(self.var.draw.maxi / binsize) * binsize nlev = round((maxout-minout) / binsize) if count gt 10 then stop if nlev lt 3 then begin binsize = binsize / 2. goto, start_here endif if nlev gt 8 then begin binsize = binsize * 2. goto, start_here endif self.var.draw.n_lev = nlev self.var.draw.mini = minout self.var.draw.maxi = maxout self -> message, string(minout, maxout, nlev, format = '("[map_image::set_range]: Mini, maxi, n_lev: ", 3f7.2)') endelse end ; ============================================================ ======================================== ; Set object variables ; ============================================================ ======================================== ;----------------------------------------------------------- ----------------------------- function map_image::set_limit d1=min(*self.var.data.longitude, i1) d2=max(*self.var.data.longitude, i2) d3=min(*self.var.data.latitude, i3) d4=max(*self.var.data.latitude, i4) i1 = where(*self.var.data.longitude eq d1) if n_elements(i1) gt 1 then i1 = i1[n_elements(i1)/2] i2 = where(*self.var.data.longitude eq d2) if n_elements(i2) gt 1 then i2 = i2[n_elements(i2)/2] i3 = where(*self.var.data.latitude eq d3) if n_elements(i3) gt 1 then i3 = i3[n_elements(i3)/2] i4 = where(*self.var.data.latitude eq d4) if n_elements(i4) gt 1 then i4 = i4[n_elements(i4)/2] limit = [$ (*self.var.data.latitude)[i1], d1, $ d4, (*self.var.data.longitude)[i4], $ (*self.var.data.latitude)[i2], d2, $ d3, (*self.var.data.longitude)[i3] $ ] return, limit end ;----------------------------------------------------------- ----------------------------- function map_image::decomposed, sss case (sss)[0] of 3:begin if (sss)[3] ne 3 then begin print, "[map_image::set_one_var]: For true color images, the image must be of dimension [xsize, ysize, 3]" stop endif self.var.draw.true_color=1 end 2: self.var.draw.true_color = (sss)[2] eq 3 1: self.var.draw.true_color = 0 else: begin print, "[map_image::set_one_var]:I don't know what to do with an image of these dimensions: " return, 1 end endcase if strupcase(!d.name) ne 'PS' then begin device, get_decomposed = decomp self.var.draw.decomposed = decomp self -> message, self.var.draw.true_color eq 1 $ ? "[map_image::set_one_var]: For true color images, I will set device, decomposed = 1" $ : "[map_image::set_one_var]: For 8-bit images, I will set device, decomposed = 0" device, decomp = self.var.draw.true_color endif if self.var.draw.true_color then self.var.draw.no_color_bar = 1 return, 0 end ;----------------------------------------------------------- ----------------------------- function map_image::limit428, limit ; Transform 4-element limit into 8-element vector ; ; If no map projection is defined yet, just assume (equidistant) ; cylindrical: x0 = limit[1] & y0 = limit[0] & x1 = limit[3] & y1 = limit[2] if !map.projection eq 0 then begin new_limit = [ $ (y0 + y1) / 2., x0, $ y1, (x0 + x1) / 2., $ (y0 + y1) / 2., x1, $ y0, (x0 + x1) / 2. $ ] endif else begin ; Normal coordinates of lower left and upper right corner: dum = convert_coord([x0, x1], [y0, y1], /data, /to_normal) ; Data coordinates of 4 points on each side: x0 = dum[0, 0] & x2 = dum[0, 1] & x1 = (x2 + x0) / 2. y0 = dum[1, 0] & y2 = dum[1, 1] & y1 = (y2 + y0) / 2. dum = convert_coord([x0, x1, x2, x1], [y1, y2, y1, y0], /normal, /to_data) ; Return correct limit vector new_limit = dum[[1, 0, 4, 3, 7, 6, 10, 9]] endelse self -> message, "[map_image::limit428]: " + $ "Making limit a 8-element vector:" self -> message, "[map_image::limit428]: " + $ "old: " + string(limit, format = '("[", 4(i4, ","), "]")') self -> message, "[map_image::limit428]: " + $ "new: " + string(new_limit, format = '("[", 8(i4, ","), "]")') return, new_limit end ;----------------------------------------------------------- ----------------------------- pro map_image::set_one_var, tag, value, $ no_copy = no_copy, $ error = error ; Look for the first matching tag name in self.var.(x) ; Some variables like limit can change their shape, so they are ; stored as pointers. ; If no tag in self.var.(x) matches, the tag and value are stored in ; self.var.internal._extra error = 0 for i = 0, n_tags(self.var) -1 do begin st_names = strlowcase(tag_names(self.var.(i))) for j = 0, n_tags(self.var.(i)) - 1 do $ if strpos(st_names[j], strlowcase(tag)) eq 0 then begin if st_names[j] eq 'image' then if self -> decomposed(size(value)) ne 0 then begin $ help, value error = 1 return endif if size(self.var.(i).(j), /type) eq 10 $ then if ptr_valid(self.var.(i).(j)) $ then *self.var.(i).(j) = value $ else self.var.(i).(j) = ptr_new(value, no_copy = no_copy) $ else self.var.(i).(j) = value if $ strpos(st_names[j] ,'ima') eq 0 and $ self.var.draw.mini eq 0 and $ self.var.draw.maxi eq 0 then self -> set_range ; Save two colours "at the bottom" for background and void color and one "at the top" for the foreground colour if st_names[j] eq 'bottom' then self.var.draw.bottom = self.var.draw.bottom > 2 if st_names[j] eq 'ncolors' or st_names[j] eq 'bottom' then self.var.draw.ncolors = self.var.draw.ncolors < (255-self.var.draw.bottom) return endif endfor if size(*self.var.internal._extra, /type) eq 8 then begin st_names = strlowcase(tag_names(*self.var.internal._extra)) for j = 0, n_tags(*self.var.internal._extra) - 1 do begin if strpos(st_names[j], strlowcase(tag)) eq 0 then begin ;self -> message, /time, "[map_image::set_one_var]: Setting self.var.internal._extra." + st_names[j] if strpos(strlowcase(tag), 'lim') eq 0 then begin if n_elements(value) eq 4 then value = self -> limit428(value) if value[0] eq -1 then value = self -> set_limit() endif (*self.var.internal._extra).(j) = value return endif endfor *self.var.internal._extra = create_struct(*self.var.internal._extra, tag, value) ;self -> message, /time, "[map_image::set_one_var]: Adding self.var.internal._extra." + tag;, ": ", value endif else begin ;self -> message, /time, "[map_image::set_one_var]: Creating self.var.internal._extra." + tag;, ": ", value *self.var.internal._extra = create_struct(tag, value) endelse end ;----------------------------------------------------------- ----------------------------- function map_image::get_var, tag ; Look for the first matching tag name in self.var.(x) ; Some variables like limit can change their shape, so they are ; stored as pointers. for i = 0, n_tags(self.var) -1 do begin st_names = strlowcase(tag_names(self.var.(i))) for j = 0, n_tags(self.var.(i)) - 1 do begin if strpos(st_names[j], strlowcase(tag)) eq 0 then begin if size(self.var.(i).(j), /type) eq 10 then begin if ptr_valid(self.var.(i).(j)) then $ return, *self.var.(i).(j) $ else $ return, -1 endif else begin return, self.var.(i).(j) endelse endif endfor endfor end ;----------------------------------------------------------- ----------------------------- pro map_image::set_var, $ error = error, $ _extra = _extra if size(_extra, /type) eq 8 then begin names = tag_names(_extra) for i = 0, n_tags(_extra) - 1 do self -> set_one_var, names[i], _extra.(i), error = error endif end ; ============================================================ ======================================== ; Internal subroutines ; ============================================================ ======================================== pro map_image::debug stop end ;----------------------------------------------------------- ----------------------------- function map_image::init, $ image, lat, lon, $ no_copy = no_copy, $ no_project = no_project, $ debug = debug, $ _extra = _extra if size(image, /type) * size(lat, /type) * size(lon, /type) eq 0 then begin print, "Image or lon or lat are not defined:" help, image, lat, lon return, 0 endif self.var.internal.debug = keyword_set(debug) self.var.internal.time = systime(1) if n_params() eq 3 then begin self -> set_one_var, 'image', image, no_copy = no_copy, error = error if error eq 1 then return, 0 self -> set_one_var, 'longitude', lon, no_copy = no_copy self -> set_one_var, 'latitude', lat, no_copy = no_copy endif else no_project = 1 ; First, fill all pointer variables with a default value just in case for i = 0, n_tags(self.var) -1 do $ for j = 0, n_tags(self.var.(i)) -1 do $ if size(self.var.(i).(j), /type) eq 10 then $ if not ptr_valid(self.var.(i).(j)) then begin ;self -> message, /time, "[map_image::init]: Setting " + (tag_names(self.var.(i)))[j] + " to " + string(i*10+j, format = '(i3)') self.var.(i).(j) = ptr_new(i*10+j) endif self -> set_one_var, 'void_index', -1 self -> set_one_var, 'void_color', 150 ; Reserve colour 0 for background and colour 1 for void-color self -> set_one_var, 'bottom', 2 self -> set_one_var, 'ncolors', 253 self.var.internal.window = !d.window ge 0 ? !d.window : 1 self.var.draw.scalef = strupcase(!d.name) eq 'PS' ? 0.03 : 1. ; Save latest color settinngs: tvlct, rold, gold, bold, /get self -> set_one_var, 'rold', rold self -> set_one_var, 'gold', gold self -> set_one_var, 'bold', bold self -> set_one_var, 'color', !p.color self -> set_one_var, 'background', !p.background if self.var.draw.no_color_bar then begin !p.position = [0.1,0.1,0.93,0.9] endif else begin vertical = 1 if size(_extra, /type) eq 8 then begin dum = strpos(strlowcase(tag_names(_extra)), 'horiz') if max(dum, index) eq 0 then vertical = (_extra).(index) eq 0 else vertical = 1 endif if vertical then begin ; Vertical colorbar !p.position = [0.24, 0.10, 0.93, 0.9] self.var.draw.cb_position = [0.13, 0.09, 0.16, 0.91] endif else begin ; Horizontal colorbar !p.position = [0.07, 0.3, 0.93, 0.9] self.var.draw.cb_position = [0.1, 0.08, 0.9, 0.14] endelse endelse self -> set_one_var, 'limit', self -> set_limit() if not keyword_set(no_project) then $ self -> project, _extra = _extra $ else $ self -> set_var, _extra = _extra return, 1 end ;----------------------------------------------------------- ----------------------------- pro map_image::cleanup, image, lat, lon if n_params() ge 1 then image = *self.var.data.image if n_params() ge 2 then lat = *self.var.data.latitude if n_params() ge 3 then lon = *self.var.data.longitude ; Find all pointers and dereference them for i = 0, n_tags(self.var) -1 do $ for j = 0, n_tags(self.var.(i)) -1 do $ if size(self.var.(i).(j), /type) eq 10 then ptr_free, self.var.(i).(j) if strupcase(!d.name) ne 'PS' then begin self -> message, string(self.var.draw.decomposed, format = '("[map_image::cleanup]: Setting device, decomposed = ", i1)') device, decomposed = self.var.draw.decomposed ; Restore original color settings: tvlct, self.var.draw.rold, self.var.draw.gold, self.var.draw.bold endif self -> message, "[map_image::colors]: !p.background: " + $ strcompress(string(!p.background), /rem) + " --> " + $ strcompress(string(self.var.draw.background), /rem) self -> message, "[map_image::colors]: !p.color: " + $ strcompress(string(!p.color), /rem) + " --> " + $ strcompress(string(self.var.draw.color), /rem) !p.color = self.var.draw.color !p.background = self.var.draw.background end ;----------------------------------------------------------- ----------------------------- pro map_image__define struct = {map_image, var: {var, $ data: {data, $ image: ptr_new(), latitude: ptr_new(), longitude: ptr_new() $ }, $ draw: {draw, $ image_title: "", legend: "", comment: "", $ mini: 0., maxi: 0., n_lev: 0, $ void_index:ptr_new(), void_color: 0, $ no_scale: 0, $ no_color_bar: 0, cb_position: [0., 0., 0., 0.], $ magnify: 0, $ draw_image: ptr_new(), xoffset: 0, yoffset: 0, scalef: 1., $ draw_image_orig: ptr_new(), $ decomposed: 0, true_color: 0, $ bottom: 0b, ncolors: 0, $ rold:bytarr(256), gold:bytarr(256), bold:bytarr(256), $ color:0l, background:0l $ }, $ internal: {internal, debug: 0, window: 0, _extra: ptr_new(), time: 0d} $ }} end ;----------------------------------------------------------- ----------------------------- ;----------------------------------------------------------- ----------------------------- pro map_image, image, lat, lon, _extra = _extra m = obj_new("map_image", image, lat, lon, _extra = _extra) obj_destroy, m end ;----------------------------------------------------------- ----------------------------- ;----------------------------------------------------------- ----------------------------- pro testdata, image, lon, lat, imgsize = imgsize imgsize = keyword_set(imgsize) ? imgsize : 500 image = dist(imgsize, imgsize) lon = fltarr(imgsize, imgsize) lat = fltarr(imgsize, imgsize) for i = 0, imgsize-1 do lon[*, i] = findgen(imgsize) / (imgsize-1) * 40. -10. for i = 0, imgsize-1 do lat[i, *] = findgen(imgsize) / (imgsize-1) * 40. + 30. end ;----------------------------------------------------------- ----------------------------- pro check_ptr q = ptr_valid(count = c) if c gt 0 then begin print, c, format = '("[example1]: Still ", i3, " valid heap variable(s)!")' q=ptr_valid(/cast) help, q, output = output for i = 0, n_elements(q) -1 do print, i, ": ", output ptr_free, ptr_valid(/cast) endif end ;----------------------------------------------------------- ----------------------------- ; EXAMPLES ;----------------------------------------------------------- ----------------------------- pro ex1, _extra = _extra ; Plain greyscale image testdata, image, lon, lat map_image, image, lat, lon, title = "Legend title / unit", _extra = _extra end pro ex2, _extra = _extra ; Low resolution data (looks bad!) testdata, image, lon, lat map_image, image, lat, lon, title = "Legend title / unit", $ limit=[48, -5, 52, 0], londel=1,latdel=1, _extra = _extra end pro ex3, _extra = _extra ; Magnify image testdata, image, lon, lat map_image, image, lat, lon, title = "Legend title / unit", $ limit=[48, -5, 52, 0], londel=1,latdel=1, $ magnify = 4, $ /rainbow, mini = 220, maxi=280, /box, _extra = _extra end pro ex4, _extra = _extra ; Different projection testdata, image, lon, lat map_image, image, lat, lon, title = "Legend title / unit", $ /ortho, p0lon = 0, p0lat = 50, /iso, /rainbow, $ limit = [20, -70, 70, 180, 20, 90, -20, 10], _extra = _extra end pro ex5, _extra = _extra ; Discrete colours and different colour table testdata, image, lon, lat, imgsize = 10 map_image, image, lat, lon, title = "Legend title / unit", $ magnify=20, mini=2, maxi=8, n_lev=6, discrete=[2, 3, 4, 5, 6, 7, 8], $ ctable = 15, bot = 20, ncol = 121, /box, _extra = _extra end pro ex6, _extra = _extra ; Discrete colors with explicitely defined colours !quiet=1 testdata, image, lon, lat, _extra = _extra, imgs=250 map_image, lon, lat, lon, title = "Longitude / degrees", $ red = [255, 0, 0, 255], green = [0, 255, 0, 255], blue = [0, 0, 255, 100], $ discrete=[0, 5, 15, 20, 25], $ /box, $ void_index = where(lon lt 0 or lon gt 25), $ min = 0, max=25, n_lev = 5, _extra = _extra end pro ex7, _extra = _extra ; Discrete colors, check for exact discrimination testdata, image, lon, lat map_image, lon, lat, lon, title="Longitude / degrees", $ mini = 4, maxi = 6, n_lev = 4, limit = [51, 4, 53, 6], londel = .5, latdel = .5, magnif = 2, $ format = '(f3.1)', /box, $ red = [122, 0, 255, 70], green = [0, 255, 70, 122], blue = [255, 70, 122, 0], $ discrete=[4, 4.5, 5, 5.5, 6.], _extra = _extra end pro ex8, _extra = _extra ; Zoom into image testdata, image, lon, lat m = obj_new("map_image", image, lat, lon, /no_copy, title = "Legend title / unit", /rainbow, _extra = _extra) m -> zoom, win = 2, magnify = 2 obj_destroy, m check_ptr !quiet=0 end pro ex9, _extra = _extra ; Multiple "overpasses" !quiet=1 testdata, image, lon, lat m = obj_new($ "map_image", image, lat, lon, $ title = "Legend title / unit", $ limit = [0, -65, 80, 10, 0, 35, -30, 10], $ /no_draw, $ /rainbow, $ _extra = _extra) m -> project, image = image, lon = lon, lat = lat - 50, /no_erase, /no_draw m -> project, image = image, lon = lon - 50, lat = lat - 50, /no_erase, /no_draw m -> project, image = image, lon = lon - 50, lat = lat, /no_erase, /no_draw m -> display obj_destroy, m check_ptr !quiet=0 end pro ex10, _extra = _extra ; Multiple images !quiet=1 testdata, image, lon, lat map_image, image, lat, lon, pos = [.05, .22, 0.45, .50], londel = 10, latdel = 10, /no_col, /rainbow, _extra = _extra map_image, image, lat, lon, /noerase, pos = [.05, .22, 0.45, .50] + [0.5, 0, 0.5, 0], londel = 20, latdel = 20, /sinu, /no_col, /rainbow, _extra = _extra map_image, image, lat, lon, /noerase, pos = [.05, .22, 0.45, .50] + [0, 0.45, 0, 0.45], /ortho, p0lon = 0, p0lat = 50, /iso, /rainbow, $ limit = [20, -70, 70, 180, 20, 90, -20, 10], /no_col, _extra = _extra m = obj_new("map_image", image, lat, lon, /noerase, pos = [.05, .22, 0.45, .50] + [0.5, 0.45, 0.5, 0.45], londel = 2, latdel = 2, limit = [55, -10, 60, -5, 55, 0, 50, -5], /no_col, magnif = 2, /rainbow, _extra = _extra) m -> colorbar, cb_position = [.04, .07, .96, .1], /horizontal, title = 'Legend title / unit' obj_destroy, m check_ptr !quiet=0 end pro ex11, _extra = _extra ; True color image testdata, image, lon, lat map_image, [bytscl(lon), bytscl(lat), bytscl(image)*0b], lat, lon, $ ; limit = [28, -15, 72, 35], $ _extra = _extra end pro ex12 ; Logarithmic colour scale testdata, image, lon, lat image = image / max(image) image = exp(image) image = (image-min(image)) / max(image) * 20 + .1 map_image, win = 1, image, lat, lon, /rainbow, min = 0, max = 20, n_lev = 4 map_image, win = 2, alog10(image), lat, lon, /rainbow, $ mini = alog10(0.1), maxi = alog10(100), $ /ylog, range = [0.1, 100], $ format = '(f5.1)', div = 3, minor = 10 end ;----------------------------------------------------------- ----------------------------- pro map_image_demo, _extra = _extra s = "" print, "[Example 1]: Just a plain greyscale image." & ex1, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex01.png", tvrd(true=1) print, "[Example 2]: Now, zoom into the image with the LIMIT keyword." & ex2, _extra = _extra ;write_png, "ex02.png", tvrd(true=1) print, "[Example 2]: This looks comparably bad, so let's try MAGNIFY" & read, s, prompt = " <ENTER> " print, "[Example 3]: Much better. We also chose a rainbow colour bar" print, "[Example 3]: and new min and max values. Plus the box around with the" print, "[Example 3]: lat / lon values." & ex3, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex03.png", tvrd(true=1) print, "[Example 4]: Now with a different projection:" & ex4, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex04.png", tvrd(true=1) print, "[Example 5]: A different colour table and discrete values:" & ex5, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex05.png", tvrd(true=1) print, "[Example 6]: Again discrete colours, but this time with explicitely defined colours:" & ex6, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex06.png", tvrd(true=1) print, "[Example 7]: And again, now check for exact discrimination of colours / values:" & ex7, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex07.png", tvrd(true=1) ;print, "[Example 8]: Interactively zoom into an image:" & ex8, _extra = _extra & read, s, prompt = " <ENTER> " print, "[Example 9]: Plot multiple datasets in one projection:" & ex9, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex09.png", tvrd(true=1) print, "[Example 10]: Plot multiple datasets in different projections, but with a common, horizontal colour bar:" & ex10, _extra = _extra & read, s, prompt = " <ENTER> " ;write_png, "ex10.png", tvrd(true=1) print, "[Example 11]: And finally, a true colour image:" & ex11, _extra = _extra ;write_png, "ex11.png", tvrd(true=1) print, "[Example 12]: Use logarithmic scaling:" & ex12, _extra = _extra ;write_png, "ex11.png", tvrd(true=1) end ############################################################ ####### # win.pro ############################################################ ####### pro win , i,xsize=xsize,ysize=ysize,_extra=_extra ;----------------------------------------------------------- -- ;+ ; NAME: ; win ; PURPOSE: ; win is a combination of wset and window ; if the window i already exists wset is used ; otherwise the window will be created. ; if device is postscript, nothing happens ; ; CATEGORY: ; Graphics ; ; CALLING SEQUENCE: ; win,1 ; ; KEYWORDS: ; all window/wset keywords ; ; MODIFICATION HISTORY: ; Written by R. Preusker, Jan, 2001. ; ; Copyright (C) 2001, Freie Universitaet Berlin ; This software may be used, copied, or redistributed as long ; as it is not sold and this copyright notice is reproduced on ; each copy made. This routine is provided as is without any ; express or implied warranties whatsoever. ; ; rene@amor.met.fu-berlin.de ;- if STRLOWCASE(!d.name) eq "ps" then return if not (keyword_set (xsize)) then xsize=400 if not (keyword_set (ysize)) then ysize=300 if (i eq 0) and (!d.window eq -1) then window,i,xsize=xsize,ysize=ysize,_extra=_extra catch, errorstat if errorstat ne 0 then window,i,xsize=xsize,ysize=ysize,_extra=_extra else wset,i end
Hi again, as promised, the files are now on a web server: http://wew.met.fu-berlin.de/idl/ I also added H5O_IS_GROUP, which is needed by READ_HDF. Cheers, Peter
I did a similar task a while back by georeferencing hundreds of hypercubes. Only difference was I had to create the array of lats/lons, which was a task in itself. This array of lats/lons (or any geocoordinates) is the equivalent of an Internal Geometry Model (IGM). Since you have an IGM for each file, yes, you can simply create a Geographic Lookup Table (GLT) using ENVI_GLT_DOIT, then warp the image using the GLT with ENVI_GEOREF_FROM_GLT_DOIT. There will be a lot of logistics around getting the data out of the HDF file, assigning FIDs and such, but what's above is the heart of what's required. Ask back if you need more help... Chris
Chris and Peter, Thanks a million for your help. I've got my code running now. I should note that I didn't use Liam's "imagemap" because of his strong caveat about the output maps not being quantitative; I have a lot of analysis to do on the mapped data. I have used imagemap (and will experiment with Peter's modification) for display and it is great. I appreciate your replies and thanks again. Barry
drdull wrote: > I've got a lot (thousands) of satellite images (HDF files) that I want > to georeference and project onto a standard grid. I want to automate > this using IDL/ENVI. Each input file includes the satellite sensor > data arrays as well as arrays with the earth location (latitude and > longitude)for each pixel. I'm interested in knowing if anyone has done > anything similar in batch mode and, if so, how you approached the > problem. Is it as simple as defining a projection, reading the data, > creating an "GLT" from the location data, and georeferencing to the new > GLT? Thanks. Can you tell us more about the satellite data? There may be a satellite-specific software package (non IDL) which would do the job more easily. Cheers, Liam. Practical IDL Programming http://www.gumley.com/