;+ ; NAME: ; hip_map_sat ; ; PURPOSE: ; Map a Meteosat or GOES-E image in GIF-FORMAT from LAPETH ; archive, or a satellite image from ECMWF archive in GRIB ; format (METEOSAT (IR,VIS,WV) or GOES (IR) are available), or ; an image from another source such as Purdue University, ; Nottingham, CMC Canada, etc. ; ; A map projection is generated simulating either the view of the earth ; from the geostationary satellite or the projection of the ; remapped image. ; ; Ask Dominik for details about Satellite-images from ECMWF! ; Note: the colors for IR images from ECMWF have to be inverted! ; ; CATEGORY: ; map-plot, satellite images ; ; CALLING SEQUENCE: ; hip_map_sat,imagetype,satfile=satfile,title=title,img=img,colmin=colmin,$ ; colmax=colmax,_extra=e,pixmap=pixmap,continents=continents,contcol=contcol ; ; EXAMPLE: ; hip_map_sat,3,satfile='/home/dominik/idl/dataexamples/ly__0630.gif' ; or ; ; INPUTS: ; imagetype: (integer) satellite image type ; METEOSAT AND GOES IMAGES FROM Ruedi's LAPETH ; ARCHIVE, partly also available throught ; http://www.wetterzentrale.de ; 1 = Europe including eastern Atlantic, type 'd2_','e2_' (IR or WV) ; 2 = Earth (Tot-Satellite-Images), type 'dtot_','etot_' (IR or WV) ; 3 = North America (GOES-E image), type 'ly_' (IR) ; 4 = Central Europe, type 'c03_' (VIS) ; 5 = GOES Spezialausschnitt fuer Paper ; METEOSAT AND GOES (possibly also GM5) IMAGES FROM ECMWF MARS ARCHIVE ; (ask Dominik for details) ; 6 = Meteosat image in binary format produced from original GRIB-File ; which can be loaded from ECMWF data archive ; 7 = same as 6 but zoom on North Atlantic and Europe ; 8 = GOES image from ECMWF (whole image) ; 9 = same as 8 but zoom on eastern US/Canada ; IMAGES FROM CANADIAN METOFFICE CMC ; CMC images ; 10 = Canadian Meteoffice (CMC) GOES-8 image type g8ecan.jpg ; 11 = same as 8 but type g8enam.jpg (larger area) ; 12 = Canadian Metoffice (CMC) GOES-9 image type g9wcan.jpg ; 13 = same as 10 but type g9wnam.jpg (larger area) ; 14 = same as 11 but zoom on upper left image part ; REMAPPED METEOSAT IMAGES ; Only available at LAPETH ; 20 = Meteosat image, Europe ; COMPOSITE SATELLITE/RADAR IMAGE FROM PURDUE ; Purdue/Unisys images ; 23 = GOES8/9 satellite or satellite+radar composite (GIF image) ; 24 = GOES8/9 radar composite (GIF image) ; IMAGES FROM NOTTINGHAM FTP SERVER ; Nottingham images ; 25 = GOES 8 IR image type LY.JPG from NOTTINGHAM ; 26 = GOES 8 VIS image type LZ.JPG from ; IMAGES FROM NOAA/NCDC HISTORICAL (12hrly) GOES BROWSER ARCHIVE ; GOES archive ; 27 = GOES 8 (east) IR/VIS image type jpg, northern hemisphere ; 28 = GOES 9 (west) IR/VIS image type jpg ; 29 = Meteosat 5 image for INDOEX (Meteosat 5 at 63deg E) ; RADAR IMAGES FROM WSI Corporation ; 30 = US radar composite ; 31 = Cummins type 1 ; 32 = Cummins type 2 ; 33 = Insat image from the India Meteorological Department. Zoom on ; India ; ; OPTIONAL INPUT PARAMETERS: ; satfile: the name of the satellite (or radar) image (format ; depending on source). If no file is selected and no ; satellite image provided, then a file can be selected! ; img Use instead of satf+ile, if image is already loaded (prevents ; reading it from "satfile" again). ; ; title: map title ; extra: Additional keywords accepted by the IDL routine map_set ; ; KEYWORD INPUT PARAMETERS: ; pixmap Use this together with the CMC-GOES images to create ; a pixmap window instead of a normal on ; continents Draw the continental outlines ; ; OUTPUTS: ; img bytarr(n,m) the satellite image read from satfile ; contcol The color of the continental outlines ; ; COMMON BLOCKS: ; ; SIDE EFFECTS: ; If this procedure is called as a subroutine and the calling routine tries to plot ; the image to the postscript device then do activate the postcript device before ; you call this routine and do not pass the parameter ps_file (has no effect) ; The procedure hip_map_sat checks, if the ps-device is already active. If 'yes' it ; does not create a new postscript file but plots to the previously defined ps-file ; ; RESTRICTIONS: ; ; PROCEDURE: ; ; MODIFICATION HISTORY: ; first version Dec. 1996 by Dominik Brunner ; DB, Many additional image sources added and updated, Jan 2000 ;- pro hip_map_sat,imagetype,satfile=satfile,title=title,img=img,colmin=colmin,$ colmax=colmax,_extra=e,pixmap=pixmap,continents=continents,contcol=contcol contthick=contthick @hip_catch_error ; ON_ERROR,2 ; check for input IF n_params() NE 1 THEN message,'Usage: hip_map_sat,imagetype,satfile=satfile,'+$ 'title=title,img=img,colmin=colmin,colmax=colmax,_extra=e,'+$ '/pixmap,/continents,contcol=contcol' IF (n_elements(satfile) eq 0) AND (n_elements(img) EQ 0) then satfile=$ pickfile(title='Select the satellite image',filter='*.gif') IF (satfile eq '') AND (n_elements(img) EQ 0) THEN BEGIN print,'no satellite image chosen - nothing done' goto,endlabel ENDIF IF n_elements(title) EQ 0 THEN title='' IF n_elements(contthick) EQ 0 THEN contthick=1 ; define minimum and maximum color occupied by satellite image IF n_elements(colmin) EQ 0 THEN colmin=1 IF n_elements(colmax) EQ 0 THEN colmax=!d.table_size-2 IF colmax LT colmin THEN BEGIN savecol=colmin colmin=colmax colmax=savecol ENDIF IF colmax GT (!d.n_colors-2) THEN colmax=!d.table_size-2 IF colmax EQ colmin THEN BEGIN print,'Error: colmax must be greater than colmin' return ENDIF ncolors=colmax-colmin+1 ; number of colors ; define satellite image file type CASE imagetype OF 1: fileformat='GIF' 2: fileformat='GIF' 3: fileformat='GIF' 4: fileformat='GIF' 5: fileformat='GIF' 6: fileformat='GRB' 7: fileformat='GRB' 8: fileformat='GRB' 9: fileformat='GRB' 10:fileformat='JPG' 11:fileformat='JPG' 12:fileformat='JPG' 20:fileformat='GIF' 23:fileformat='GIF' 24:fileformat='GIF' 25:fileformat='JPG' 26:fileformat='JPG' 27:fileformat='JPG' 28:fileformat='JPG' 29:fileformat='JPG' 30:fileformat='GIF' 31:fileformat='GIF' 32:fileformat='GIF' 33:fileformat='GIF' ELSE: BEGIN print,'image type is not available/implemented' goto,endlabel END ENDCASE ;load satellite image and plot it IF n_elements(img) LE 1 THEN BEGIN ; check fileformat again IF strpos(strlowcase(satfile),'gif') NE -1 THEN fileformat='GIF' IF strpos(strlowcase(satfile),'jpg') NE -1 THEN fileformat='JPG' IF strpos(strlowcase(satfile),'jpeg') NE -1 THEN fileformat='JPG' tvlct,r0,g0,b0,/get ; load current colortable invert=0 CASE fileformat OF 'GIF': BEGIN read_gif,satfile,satimage,r,g,b minpix=min(satimage,max=maxpix) r=r[minpix:maxpix] & g=g[minpix:maxpix] & b=b[minpix:maxpix] ; load colortable provided by the gif-image into the ; colmin:colmax range (using linear interpolation) r0[colmin:colmax]=congrid(r,ncolors) g0[colmin:colmax]=congrid(g,ncolors) b0[colmin:colmax]=congrid(b,ncolors) tvlct,r0,g0,b0 ; load modified colortable END 'GRB': BEGIN ; ECMWF satellite image in GRIB format gribverbose,satfile,info=info,varlist=varlist,ok=ok IF NOT ok THEN return grib_read,info,varlist,0,filename=satfile,field=satimage satimage=reverse(satimage,0) print,'Note: Eventually the present colortable needs' print,'to be adapted (usually color index range 50-150)' IF strpos(strlowcase(satfile),'ir') NE -1 THEN satimage=255B-satimage END 'JPG': BEGIN ; JPEG image read_jpeg,satfile,satimage,ctable,colors=colmax-colmin+1,$ /TWO_PASS_QUANTIZE r=ctable[*,0] & g=ctable[*,1] & b=ctable[*,2] ; load colortable provided by the JPG-image into the ; colmin:colmax range (using linear interpolation) IF n_elements(r) LE (colmax-colmin+1) THEN BEGIN r0[colmin:colmin+n_elements(r)-1]=r g0[colmin:colmin+n_elements(g)-1]=g b0[colmin:colmin+n_elements(b)-1]=b ENDIF ELSE BEGIN r0[colmin:colmax]=congrid(r,ncolors) g0[colmin:colmax]=congrid(g,ncolors) b0[colmin:colmax]=congrid(b,ncolors) ENDELSE ;IF (imagetype LT 25) OR (imagetype GT 26) THEN tvlct,r0,g0,b0 tvlct,r0,g0,b0 END ELSE: exit ENDCASE ; rescale the bitmap to the color range minpix=min(satimage,max=maxpix) scale=float(colmax-colmin)/float(maxpix-minpix) satimage=byte(colmin+scale*(satimage-minpix)) img=satimage ENDIF ELSE satimage=img ; define special settings for individual image types ; such as IMAGE type, mapping area, etc. satp0lat=0 ; the latitude of the geostationary satellites CASE imagetype OF 1: BEGIN projtype='GEOSTAT' satP0lon=0. ; Meteosat ; Check if image has correct size IF n_elements(satimage) NE 640000L THEN BEGIN print,'Error: D2 image must be 800x800 pixels' return ENDIF ;clip lower image part (Africa) satimage=satimage[*,245:799] ; the full size is 800 x 800 pixels ;limits for map-projection [lats,lonw,latn,lone] limit=[28.4,-18.8,73.7,18.8] END 2: BEGIN projtype='GEOSTAT' satP0lon=0. ; Meteosat ; Check if image has correct size IF n_elements(satimage) NE 640000L THEN BEGIN print,'Error: DTOT image must be 800x800 pixels' return ENDIF ;clip lower image part (to 10degN and some of the right and left sides) satimage=satimage[34:766,499:799] limit=[13.5,-59,72.5,59] END 3: BEGIN ; equivalent to 25 projtype='GEOSTAT' satP0lon=-75. ; GOES8 ; Check if image has correct size IF n_elements(satimage) NE 640000L THEN BEGIN print,'Error: LY image must be 800x800 pixels' return ENDIF ;clip lower image part (North America) satimage=satimage[*,249:799] ;limits for map-projection [lats,lonw,latn,lone] limit=[13.25,-111.5,73,-53.8] END 4: BEGIN projtype='GEOSTAT' satP0lon=0. ; Meteosat subregion C03 ;limits for map-projection [lats,lonw,latn,lone] limit=[35.8,0.3,70.4,21] END 5: BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES8 ;clip lower image part (South America) ; Check if image has correct size IF n_elements(satimage) NE 640000L THEN BEGIN print,'Error: LY image must be 800x800 pixels' return ENDIF satimage=satimage[0:599,350:799] ;limits for map-projection [lats,lonw,latn,lone] limit=[20.8,-113,74,-67.7] END 6: BEGIN projtype='GEOSTAT' satP0lon=0. ; Meteosat (from MARS archive) IF n_elements(satimage) NE 6250000L THEN BEGIN print,'Error: Meteosat image from archive must be 2500x2500 pixels' return ENDIF ; entire view (only horizon is clipped) xclip=22 yclip=27 xsize=n_elements(satimage[*,0]) ysize=n_elements(satimage[0,*]) satimage=satimage(xclip:xsize-xclip-1,yclip:ysize-yclip-1) END 7: BEGIN projtype='GEOSTAT' satP0lon=0. ; Meteosat (from MARS archive) IF n_elements(satimage) NE 6250000L THEN BEGIN print,'Error: Meteosat image from archive must be 2500x2500 pixels' return ENDIF ; zoom on Europe/North Atlantic satimage=satimage[800:1610,1900:2450] ;limits for map-projection [lats,lonw,latn,lone] limit=[-40,40,40,120] END 8: BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES-E (8) IF n_elements(satimage) NE 1562500L THEN BEGIN print,'Error: GOES-E image from archive must be 1250x1250 pixels' return ENDIF ; entire view (only horizon is clipped) xclip=11 yclip=12 xsize=n_elements(satimage(*,0)) ysize=n_elements(satimage(0,*)) satimage=satimage[xclip:xsize-xclip-1,yclip:ysize-yclip-1] END 9: BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES-E (8) IF n_elements(satimage) NE 1562500L THEN BEGIN print,'Error: GOES-E image from archive must be 1250x1250 pixels' return ENDIF ; zoom on US/CANADA satimage=satimage[400:805,950:1225] ;limits for map-projection [lats,lonw,latn,lone] limit=[28.7,-21.3-75.,90,17.-75.] END 10:BEGIN ; CMC, GOES 8, type ecan projtype='STEREO1' ; North America limit=[29.25,-75-15.,61.8,-75+15.] p0lat=90. & p0lon=-75. END 11:BEGIN ; CMC, GOES 8, type enam projtype='STEREO1' ; North America limit=[13.7,-101.,58.2,-49.1] p0lat=90. & p0lon=-75. END 12:BEGIN projtype='STEREO1' ; North America limit=[29.25,-132.,61.5,-98] p0lat=90. & p0lon=-115. END 14:BEGIN ; CMC, GOES 8, type enam projtype='STEREO1' ; North America satimage=satimage[0:511,512:1023] limit=[36.2,-117.75,61.,-75.1] p0lat=90. & p0lon=-75. END 20:BEGIN projtype='STEREO3' ; Meteosat image from SMA limit=[28.,-23.5,67.85,18.1] p0lat=90. & p0lon=0. END 23:BEGIN projtype='STEREO2' ;limits for map-projection [lats,lonw,latn,lone] ;limit=[21.,-110.,55.,-85.] limit=[23.,-116.6,54.7,-77.4] p0lat=90. & p0lon=-97. END 24:BEGIN projtype='STEREO3' ;limits for map-projection [lats,lonw,latn,lone] ;limit=[21.,-110.,55.,-85.] ;satimage=satimage[0:898,26:623] limit=[24.3,-118.2,52.7,-75.9] p0lat=90. & p0lon=-97. END 25:BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES8 ;clip lower image part (South America) satimage=satimage[*,249:799] ;limits for map-projection [lats,lonw,latn,lone] limit=[13.25,-111.5,73,-53.8] END 26:BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES8 ;limits for map-projection [lats,lonw,latn,lone] limit=[17.1,-87.5,73.,-53.3] END 27:BEGIN projtype='GEOSTAT' satP0lon=-75. ; GOES8 ;limits for map-projection [lats,lonw,latn,lone] limit=[2.3,-110.8,59.,-39.] END 28:BEGIN projtype='GEOSTAT' satP0lon=-135. ; GOES9 satimage=satimage[64:n_elements(satimage[*,0])-1,*] ;limits for map-projection [lats,lonw,latn,lone] limit=[2.3,-107.4-60.,59.,-42.6-60.] END 29: BEGIN projtype='GEOSTAT' satP0lon=63. ; METEOSAT 5 (new position over Indian Ocean) ; entire view (only horizon is clipped) ;limits for map-projection [lats,lonw,latn,lone] limit=[-34.,31.,34.,108.5] END 30: BEGIN projtype='MERC1' ;limits for map-projection [lats,lonw,latn,lone] limit=[20.6,-121.4,55.7,-69.4] p0lat=36. & p0lon=-95. END 31:BEGIN projtype='MERC1' ;limits for map-projection [lats,lonw,latn,lone] ;limit=[21.,-110.,55.,-85.] p0lat=40. p0lon=-96.5 ; limit=[19.2,-117.88333,54.0333,-60.28333] limit=[22.,-118.25,56.5,-74.78] satimage=satimage[0:550,*] END 32:BEGIN projtype='STEREO2' ;limits for map-projection [lats,lonw,latn,lone] ;limit=[21.,-110.,55.,-85.] limit=[20.52,-117.37,55.,-36.54] p0lat=90. & p0lon=-97. END 33: BEGIN projtype='CYL1' ; Insat image in cylindrical projection ; Check if image has correct size IF n_elements(satimage) NE 281600L THEN BEGIN print,'Error: Insat image must have 512x550 pixels!' return ENDIF ;limits for map-projection [lats,lonw,latn,lone] limit=[-9.45,45.,48.1,104.9] END ENDCASE ;check if postscript device is active ps_active=!d.flags AND 1 IF NOT ps_active THEN BEGIN ; output to screen IF (projtype EQ 'STEREO1') THEN BEGIN ; create quadratic window size=max([!d.x_size,!d.y_size]) window,xsize=size+total(!x.omargin+1)*!d.x_ch_size,$ ysize=size+total(!y.omargin+1)*!d.y_ch_size,pixmap=pixmap ENDIF IF projtype EQ 'STEREO2' THEN BEGIN ; create window with aspect ratio 640:512 size=!d.y_size window,xsize=size*640/512+total(!x.omargin+1)*!d.x_ch_size,$ ysize=size+total(!y.omargin+1)*!d.y_ch_size,pixmap=pixmap ENDIF IF projtype EQ 'MERC1' THEN BEGIN window,xsize=n_elements(satimage[*,0]),ysize=n_elements(satimage[0,*]) ENDIF ; define image size in pixels xsize=round(!d.x_size-total(!x.omargin+1)*!d.x_ch_size) ysize=round(!d.y_size-total(!y.omargin+1)*!d.y_ch_size) satimage=congrid(satimage,xsize,ysize) tv,satimage,(!x.omargin[0]+1)*!d.x_ch_size,$ (!y.omargin[0]+1)*!d.y_ch_size ENDIF ELSE BEGIN tv,satimage,(!x.omargin[0]+1)*!d.x_ch_size,(!y.omargin[0]+1)*!d.y_ch_size,$ xsize=!d.x_size-total(!x.omargin+1)*!d.x_ch_size,$ ysize=!d.y_size-total(!y.omargin+1)*!d.y_ch_size ENDELSE CASE projtype OF 'GEOSTAT': map_set,satP0lat,satP0lon,/satel,sat_p=[6.61,0,0],limit=limit,$ title=title,/noerase,xmargin=1,ymargin=1,_extra=e 'STEREO1': map_set,p0lat,p0lon,limit=limit,/noerase,/ster,/iso,$ /grid,londel=10.,latdel=10.,xmargin=1,ymargin=1,/nob,_extra=e 'STEREO2': map_set,p0lat,p0lon,limit=limit,/noerase,/ster,/iso,/nob,_extra=e 'STEREO3': map_set,p0lat,p0lon,limit=limit,/noerase,/ster,/nob 'MERC1': map_set,p0lat,p0lon,limit=limit,/noerase,/mercator,_extra=e,$ xmargin=1,ymargin=1 'CYL1': map_set,limit=limit,/noerase,/mercator,_extra=e,$ xmargin=1,ymargin=1 ELSE:return ENDCASE endlabel: IF keyword_set(continents) THEN BEGIN IF n_elements(contcol) EQ 1 THEN $ map_continents,/hires,color=contcol,thick=contthick $ ELSE map_continents,thick=contthick,/hires ENDIF end