ENVI二次开发的内部（猜测）示例代码分析

;--------------------------------------------------------+

; NAME

;

;   GEOREF_RSAT

;

; TYPE

;

;   ENVI User Function

;

; CALLING SEQUENCE

;

;   GEOREF_RSAT, event

;

; PURPOSE

;

;   An ENVI User Function to automatically georeference

;   (register) radarsat format image data using the embedded

;   ground control points.

;

; MODIFICATION HISTORY

;

;   April 1999     -- written under ENVI 3.0, IDL 5.03

;   September 2000 -- added About Button and a check for missing GCPs

;

; DISTRIBUTION INFO

;

;   Distributed from the RSI Web Site, see Tech Tip

;   ENVI249 "Automatically Registering Radarsat

;   Imagery from Embedded GCPs" at:

;   http://www.rsinc.com/services/search.cfm

;

; AUTHOR INFO

;

;   David Gorodetzky

;   RSI Profession Services Group

;   goro@rsinc.com

;

; PROGRAM NOTES

;

;   According to CCRS documentation available at the time this

;   routine was written, the lat/lon of the first, middle and

;   last pixel in each line of the radarsat image is embedded

;   into the line prefix in bytes 132 to 156.  The geocoded

;   data are 4 byte integers arranged as follows:

;

;   1st_lat, middle_lat, last_lat, 1st_lon, middle_lon, last_lon

;

;   The units for the lat/lon data are millionths of a

;   degree.  Note that the accuracy of the registration

;   is limited by the accuracy of the embedded GCPs, which

;   are estimated to be about 400 m.  Also, the datum of

;   the lat/lon coordinates is unknown and assumed to be

;   WGS84.  The use of the embedded geocoded data to georerference

;   the image appears to be a good first step in getting the geometry

;   correct, although the absolute positional accuracy of the

;   registered result is often in error by as much as 1500 m.

;   Translation errors like this can be corrected if an additional

;   source od ground control is available, such as an accurate

;   vector layer.

;

;   This file contains the following routines in this order:

;

;   GEOREF_RSAT_EVENT

;

;     The event handler for the first GUI.  It stores all user

;     defined data in a pointer called info.

;

;   GEOREF_RSAT

;

;     The processing routine.  It also makes the second GUI

;     which is auto-managed by ENVI.

;

;---------------------------------------------------------------+



FORWARD_FUNCTION ENVI_FILE_TYPE, WIDGET_AUTO_BASE, WIDGET_PMENU, $

  AUTO_WID_MNG, WIDGET_PARAM, WIDGET_PROJ, WIDGET_OUTF, WIDGET_MENU, $

  WIDGET_OUTFM, ENVI_TRANSLATE_PROJECTION_UNITS, MAGIC_MEM_CHECK, $

  ENVI_REAL_MAIN_BASE

 

;----------------------------------------------------------------

PRO GEOREF_RSAT_EVENT, event

  ;

  ; Event handler for the first GUI.

  ;

  ;---------------------------------------------------------------+



  WIDGET_CONTROL, event.top, get_uvalue=info

  WIDGET_CONTROL, event.id, get_uvalue=what

 

  CASE what OF

 

    'OK':     BEGIN

      ;

      ; transfer user inupted data into info pointer

      ;

      WIDGET_CONTROL, (*info).num_gcp, get_value=skip

      (*info).skip = skip

      WIDGET_CONTROL, (*info).proj, get_value=sel_proj

      (*info).sel_proj = sel_proj

      WIDGET_CONTROL, (*info).outf, get_value=outfile

      (*info).outfile = outfile

      WIDGET_CONTROL, (*info).pts_only, get_value=get_pts_only

      (*info).get_pts_only = get_pts_only

      ;

      ; destroy the GUI and return

      ;

      WIDGET_CONTROL, event.top, /destroy

      RETURN

    END

    'CANCEL': BEGIN

      (*info).cancel = 1

      WIDGET_CONTROL, event.top, /destroy

      RETURN

    END

    'ABOUT':  BEGIN

      ;

      ; display about message text

      ;

      msg = STRARR(28)

      msg[0] = "GEOREF_RSAT User Function"

      msg[1] = "----------------------------------"

      msg[2] = ""

      msg[3] = "Written by David Gorodetzky, goro@rsinc.com"

      msg[4] = 'Research Systems, Inc. (RSI)'

      msg[5] = "Professional Services Group (PSG)"

      msg[6] = "http://www.rsinc.com"

      msg[7] = ""

      msg[8] = "For additional information on this function please"

      msg[9] = "see the Tech Tip article 'Automatically Registering"

      msg[10] ="Radarsat Imagery from Embedded GCPs' on the RSI web"

      msg[11] ="site at http://www.rsinc.com/services/search.cfm"

      msg[12] =""

      msg[13] ="For more information about consulting services"

      msg[14] ="related to ENVI, IDL or remote sensing contact"

      msg[15] ="RSI at (303) 786-9900 and ask for the"

      msg[16] ="Professional Services Group, or send an e-mail"

      msg[17] ="to << consult@rsinc.com >>"

      msg[18] = ""

      msg[19] ='DISCLAIMER'

      msg[20] = ""

      msg[21] ='  This additional functionality is provided free'

      msg[22] ='  of charge as a service for our ENVI users.'

      msg[23] ='  However, this procedure is not supported by'

      msg[24] ='  Research Systems Technical Support.  The'

      msg[25] ='  procedure GEOREF_RSAT has been tested'

      msg[26] ='  and we believe that it works correctly, but it'

      msg[27] ='  is in no way guaranteed.'

     

      ENVI_INFO_WID, msg, title='About RSI User Functions'

     

    END

    ELSE:

  ENDCASE

 

END



;----------------------------------------------------------------

PRO GEOREF_RSAT, event

  ;

  ; The processing routine.  This also makes the second GUI

  ; which is auto-managed by ENVI.

  ;

  ;---------------------------------------------------------------+



  ; use CATCH to prevent unexpected errors from crashing ENVI

  ;

  !error_state.code = 0

  CATCH, error

  IF (error NE 0) THEN BEGIN

    IF (N_ELEMENTS(rep_base) NE 0) THEN ENVI_REPORT_INIT, base=rep_base, /finish

    ok = DIALOG_MESSAGE(!error_state.msg, /error, /cancel)

    IF (STRUPCASE(ok) EQ 'CANCEL') THEN BEGIN

      IF (N_ELEMENTS(lun) NE 0) THEN FREE_LUN, lun

      RETURN

    ENDIF

  ENDIF

 

  ; select the input image, disable spatial and spectral

  ; subsetting and force the user to select by band -- also

  ; require that the input image is a RADARSAT file type

  ;

  ft = ENVI_FILE_TYPE("RADARSAT")

  ENVI_SELECT, title='Select RADARSAT Input File', $

    fid=fid, /no_dims, /no_spec, /band_only,$; file_type=ft, $

    dims=dims, pos=pos

  IF (fid EQ -1) THEN RETURN

 

  ; get some basic info about the file

  ;

  ENVI_FILE_QUERY, fid, ns=ns, nl=nl, offset=offset, $

    fname=fname, data_type=dt, sname=sname

  OPENR, lun, fname, /Get_LUN

  fsize = FSTAT(lun)

 

  ; determine the number of bytes per image pixel (bpip)

  ; based on the difference between the image size and

  ; the file size, accounting for the line prefix of 192 bytes

  ;

  CASE 1 OF

    fsize.size EQ ( (ns+192)*(nl) + offset):       bpip = 1

    fsize.size EQ ( ((ns*2)+192)*nl + offset):       bpip = 2

    ELSE:        BEGIN

      extra_tlb = WIDGET_AUTO_BASE(title='specify data type')

      txt1 = 'Unable to compute the data type of the imagery.'

      txt2 = 'Select the number of bytes per image pixel.  If'

      txt3 = 'the data file was scaled to byte when it was'

      txt4 = 'opened remember to choose the original data type.'

      sbase0 = WIDGET_BASE(extra_tlb, /col, xsize=255, /frame)

      msg = WIDGET_LABEL(sbase0, value=txt1, /align_left)

      msg = WIDGET_LABEL(sbase0, value=txt2, /align_left)

      msg = WIDGET_LABEL(sbase0, value=txt3, /align_left)

      msg = WIDGET_LABEL(sbase0, value=txt4, /align_left)

      sbase1 = WIDGET_BASE(sbase0, /row)

      choose = WIDGET_PMENU(sbase1, /auto, uvalue='bpip', default=0, $

        list = ['1  (byte)', '2  (integer)', '4  (long integer/float)'], $

        prompt='bytes per image pixel')

      lookup = [1,2,4]

      extra_result = AUTO_WID_MNG(extra_tlb)

      IF (extra_result.accept EQ 0) THEN RETURN

      bpip = lookup(extra_result.bpip)

    END

  ENDCASE

 

  ; create the GUI

  ;

  ENVI_CENTER, x, y

  TLB = WIDGET_BASE(/col, xoffset=x, yoffset=y-100, /modal, $

    group_leader=envi_real_main_base(), title='Extract RADARSAT GCP Parameters')

   

  sbase1 = WIDGET_BASE(TLB, /frame, /col)

  line_skip_txt1 = 'Each image line yields 3 GCPs (first, middle, and last pixel).'

  line_skip_txt2 = 'Set the number of lines to skip between extracting GCPs.'

  label2 = WIDGET_LABEL(sbase1, value=line_skip_txt1, /align_left)

  label3 = WIDGET_LABEL(sbase1, value=line_skip_txt2, /align_left)

 

  rbase1 = WIDGET_BASE(sbase1, /row)

  num_gcp = WIDGET_PARAM(rbase1, prompt='line skip factor ', $

    uvalue = 'skip', xs=3, dt = 2, default=10, floor=0, $

    ceil=(nl/2), /all_events)

   

  proj = WIDGET_PROJ(TLB, proj=proj, uvalue='proj', $

    prompt='select output projection for registeration', /frame)

   

  sbase2 = WIDGET_BASE(TLB, /frame, /col)

  def_fname = STRMID(sname, 0, RSTRPOS(sname, '.'))+'.pts'

  outfile = WIDGET_OUTF(sbase2, prompt='output .pts file name', $

    uvalue='outfile', default=def_fname)

   

  sbase3 = WIDGET_BASE(sbase2, /row)

  pts_only = WIDGET_MENU(sbase2, uvalue='pts_only', $

    list=["don't georeference, just write points file"])

   

  row_base = WIDGET_BASE(TLB, /row)

  ok = WIDGET_BUTTON(row_base, value='  OK  ', uvalue='OK')

  can = WIDGET_BUTTON(row_base, value='Cancel', uvalue='CANCEL')

  about = WIDGET_BUTTON(row_base, value='About This Function', uvalue='ABOUT')

 

  WIDGET_CONTROL, TLB, /realize

  info_data = { $

    num_gcp:num_gcp, $         ; widget ID for line skip factor widget_param

    skip:10, $              ; default line skip factor

    proj:proj, $            ; widget ID for widget_proj

    sel_proj:{envi_proj_struct}, $; projection selection

    outf:outfile, $          ; widget ID for the OUTF widget for pts file

    outfile:def_fname, $       ; default .pts filename

    pts_only:pts_only, $       ; widget ID for the checkbox for doing pts only

    get_pts_only:0, $          ; 0 = georef, 1=only get the .pts

    CANCEL:0 $              ; cancel instructions

    }

   

  info = PTR_NEW(info_data, /no_copy)

  WIDGET_CONTROL, TLB, set_uvalue=info

  XMANAGER, 'GEOREF_RSAT', TLB, /no_block

 

  ; return if the user cancels

  ;

  IF (*info).cancel EQ 1 THEN RETURN

 

  IF ( (*info).get_pts_only NE 1) THEN BEGIN  ; if doing the full registration

 

    TLB = WIDGET_AUTO_BASE(title='Registration Parameters', /ybig)

   

    sbase0 = WIDGET_BASE(TLB, /col, /frame)

   

    sbase1 = WIDGET_BASE(sbase0, /row)

    model = WIDGET_PMENU(sbase1, /auto, default=2, prompt='warp model', $

      list=['RST', 'Polynomial', 'Triangulation'], uvalue='model')

     

    sbase2 = WIDGET_BASE(sbase0, /row)

    resample = WIDGET_PMENU(sbase2, /auto, default=0, prompt='resampling method', $

      list=['Nearest Neighbor', 'Bilinear Interpolation', 'Cubic Convolution'], $

      uvalue='resample')

     

    sbase4 = WIDGET_BASE(sbase0, /row)

    back = WIDGET_PARAM(sbase4, prompt='background value', $

      uvalue = 'back', xs=5, dt = 4, field=1, default=0, /auto)

     

    sbase5 = WIDGET_BASE(sbase0, /row)

    pixX = WIDGET_PARAM(sbase5, /auto, prompt='Output Pixel Sizes: X ', $

      uvalue='pixX', xs=5, dt=4, field=1, default=100)

    pixY = WIDGET_PARAM(sbase5, /auto, prompt=' Y ', $

      uvalue='pixY', xs=5, dt=4, field=1, default=100)

     

    warpfile = WIDGET_OUTFM(TLB, /auto, uvalue='warpfile', $

      /frame, prompt='Registered Image Filename')

     

    result2 = AUTO_WID_MNG(TLB)

    IF result2.accept EQ 0 THEN RETURN

   

  ENDIF

 

  ; convert the widget info back into the correct data type

  ;

  (*info).skip = LONG((*info).skip)

 

  ; set up varaibles and arrays needed to extract the GCPs

  ;

  num_gcp_lines = LONG( (nl-1) / ((*info).skip+1) + 1 )

  full_line = (ns*bpip)+192L

  temp = LONARR(6)

  lat = LONARR(3,num_gcp_lines)

  lon = LONARR(3,num_gcp_lines)

  pix_y = LONARR(3,num_gcp_lines)

 

  ; set up the status reporting widget

  ;

  rstr = ['Input File: '+fname, 'Output File: '+(*info).outfile]

  ENVI_REPORT_INIT, rstr, base=rep_base, /interupt, $

    title='extracting GCPs from RADARSAT'

  ENVI_REPORT_INC, rep_base, num_gcp_lines+3

  ENVI_REPORT_STAT, rep_base, 1, num_gcp_lines+3

 

  ; read in the first image line's GCPs

  ;

  POINT_LUN, lun, offset+132

  READU, lun, temp

  IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)

  lat(*,0) = temp(0:2)

  lon(*,0) = temp(3:5)

  pix_y(*,0) = 1

 

  ; read in the rest of the GCPs, except for the last line

  ;

  line = (*info).skip + 2

  file_pointer = ( (line-1)*full_line ) + 132 + offset

  FOR array_row=1, num_gcp_lines-2 DO BEGIN

    ;

    ; update the status reporting widget and

    ; check for the function being cancelled

    ;

    ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel

    IF (cancel) THEN BEGIN

      FREE_LUN, lun

      ENVI_REPORT_INIT, base=rep_base, /finish

      RETURN

    ENDIF

    ;

    ; extract the rest of the GCPs except for the last line

    ;

    POINT_LUN, lun, file_pointer

    READU, lun, temp

    IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)

    lat(*,array_row) = temp(0:2)

    lon(*,array_row) = temp(3:5)

    pix_y(*,array_row) = line

    line = line + (*info).skip + 1

    file_pointer = file_pointer + ((*info).skip+1)*full_line

  ENDFOR

 

  ; read in the last image line's GCPs

  ;

  ; first, update the status report widget

  ;

  ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel

  IF (cancel) THEN BEGIN

    FREE_LUN, lun

    ENVI_REPORT_INIT, base=rep_base, /finish

    RETURN

  ENDIF

  ;

  lastline = offset + (nl-1)*full_line + 132

  POINT_LUN, lun, lastline

  READU, lun, temp

  IF ( BYTE(255,1) EQ 0 ) THEN temp = SWAP_ENDIAN(temp)

  lat(*,array_row) = temp(0:2)

  lon(*,array_row) = temp(3:5)

  pix_y(*,array_row) = nl

 

  ; don't forget to free the file pointer

  ;

  FREE_LUN, lun

 

  ; make the array that will become the .pts file

  ;

  ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel

  IF (cancel) THEN BEGIN

    ENVI_REPORT_INIT, base=rep_base, /finish

    RETURN

  ENDIF

  GCPs = DBLARR(4, num_gcp_lines*3)

  ;

  ; define the x-pixel positions (1st, middle, last)

  ;

  pix_x = LONARR(3,num_gcp_lines)

  pix_x(0,0) = 1L

  pix_x(1,0) = LONG(ns/2)

  pix_x(2,0) = ns

  pix_x(0,*) = pix_x(0,0)

  pix_x(1,*) = pix_x(1,0)

  pix_x(2,*) = pix_x(2,0)

  ;

  ; convert the lat/lons into the selected output projection,

  ; remember to convert them from millionths to degrees first

  ;

  ENVI_PROJ_CONVERT, map_x, map_y, REFORM( (lat/1e6), num_gcp_lines*3), $

    REFORM( (lon/1e6), num_gcp_lines*3), projection=(*info).sel_proj, /to_map

  ;

  ; if the file is missing the embedded GCPs then they usually

  ; end up being zeros, so do a quick check to see if the data

  ; are OK -- if more than 50% of the GCPs are zero then

  ; warn the user

  ;

  w1 = WHERE(map_x EQ 0, c1)

  w2 = WHERE(map_y EQ 0, c2)

  IF ( c1 GT (0.5*3*num_gcp_lines) OR c2 GT (0.5*3*num_gcp_lines) ) THEN BEGIN

    warn = DIALOG_MESSAGE("Most of the GCPs are zeros!  You're Radarsat image "+$

      'is probably missing the GCPs.  Do you want to Continue?', /question)

    IF STRUPCASE(warn) EQ 'NO' THEN BEGIN

      FREE_LUN, lun

      ENVI_REPORT_INIT, base=rep_base, /finish

      RETURN

    ENDIF

  ENDIF

  ;

  ; put the GCP data into the .pts array

  ;

  GCPs(0,*) = map_x

  GCPs(1,*) = map_y

  GCPs(2,*) = REFORM(pix_x, num_gcp_lines*3)

  GCPs(3,*) = REFORM(pix_y, num_gcp_lines*3)

 

  ; write the .pts file to disk

  ;

  OPENW, lun, (*info).outfile, /Get_LUN

  ENVI_REPORT_STAT, rep_base, array_row+1, num_gcp_lines+3, cancel=cancel

  IF (cancel) THEN BEGIN

    FREE_LUN, lun

    ENVI_REPORT_INIT, base=rep_base, /finish

    RETURN

  ENDIF

  PRINTF, lun, 'ENVI Registration GCP File'

  PRINTF, lun, GCPs, format="(32000(2(F20.2), 2(I20)/))"

  FREE_LUN, lun

 

  ENVI_REPORT_INIT, base=rep_base, /finish

 

  ; as long as the user doesn't choose to write

  ; the points file only, do the warp

  ;

  IF ((*info).get_pts_only NE 1) THEN BEGIN

 

    ; x0 and y0 are the map coordinates of the upper left

    ; pixel of the output warped image -- these should be

    ; set to the Westernmost (X) and Northernmost (Y) GCPs

    ;

    x0 = MIN( GCPs(0,*) )

    y0 = MAX( GCPs(1,*) )

   

    ; xsize and ysize are the size of the ouput image in map coordinates

    ;

    xsize = MAX(GCPs(0,*)) - MIN(GCPs(0,*))

    ysize = MAX(GCPs(1,*)) - MIN(GCPs(1,*))

   

    ; define the map information structure

    ;

    map_info = {envi_map_struct}

    map_info.proj = (*info).sel_proj

    map_info.ps = [result2.pixX, result2.pixY]

    map_info.mc = [0, 0, x0, y0]

    CASE map_info.proj.type OF

      1:      map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('degrees')    ; Geographic Lat/Lon gets degrees

      8:      map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('feet')        ; State Plane gets feet

      ELSE:   map_info.proj.units = ENVI_TRANSLATE_PROJECTION_UNITS('meters')     ; everything else gets meters

    ENDCASE

   

    ; if the user wants to output the result

    ; to memory, make sure it will not exceed

    ; ENVI's cache limit

    ;

    mem_dims = [-1L, 0, LONG(xsize/result2.pixX), $

      0, LONG(ysize/result2.pixY) ]

    mem_result = MAGIC_MEM_CHECK(dims=dims, fid=fid, nb=N_ELEMENTS(pos), $

      in_memory=result2.warpfile.in_memory, out_dt=dt, $

      out_name=result2.warpfile.name)

    IF (mem_result.cancel EQ 1) THEN RETURN

   

    ; set the warp method

    ;

    CASE result2.model OF

      0:  BEGIN

        CASE result2.resample OF

          0:  method=0     ; RST w/nearest neighbor

          1:  method=1     ; RST w/bilinear

          2:  method=2     ; RST w/cubic convolution

          ELSE:

        ENDCASE

      END

      1:  BEGIN

        CASE result2.resample OF

          0:  method=3     ; Polynomial w/nearest neighbor

          1:  method=4     ; Polynomial w/bilinear

          2:  method=5     ; Polynomial w/cubic convolution

          ELSE:

        ENDCASE

      END

      2:  BEGIN

        CASE result2.resample OF

          0:  method=6     ; Triangulation w/nearest neighbor

          1:  method=7     ; Triangulation w/bilinear

          2:  method=8     ; Triangulation w/cubic convolution

          ELSE:

        ENDCASE

      END

      ELSE:

    ENDCASE

   

    ; make sure the correct save file is restored

    ;

    ENVI_CHECK_SAVE, /registration

   

    ; do the registration warp, note that if the user chooses

    ; to do a Polynomial warp, they always get 1st order

    ;

    REG_WARP_DOIT, fid=fid, pos=pos, dims=dims, map_info=map_info, $

      out_name=mem_result.out_name, method=method, degree=1, $

      in_memory=mem_result.in_memory, pts=GCPs, background=result2.back, $

      to_image=0, xsize=xsize, ysize=ysize, xstart=0, ystart=0, $

      x0=x0, y0=y0, pixel_size=map_info.ps     

  ENDIF

 

END