Mask/CLIP Raster images with different resolution in IDL

Question

I want to clip/mask raster image (500meters resolution) by another raster images (10 km resolution) using IDL Programming after the clip/mask process image should be in 500 meters resolution. I have 365 images pair and I want to a process by IDL programming. Can anybody write this code in IDL? Thanks in advance.

Answer 1

@MrFuppes

Initially, I write one code in IDL which able to mask/clip the same resolution images because those are in same dimension. But I cannot write the mask function for different resolution where is the problem in dimension and resolution.

In my code:

Band 1(500 meters) is my base images and I want to mask this image by cloud mask image(500 meters)(It works well) and then I want to do another mask by AOD_MASK images(10 km) resolution. But how to write MASK function for different dimension images I have no idea. If you have any suggestion regarding this issue then, please suggested me.

Here I attached my test DATA please check it.

Below is my IDL code:

 `Pro MASK

  COMPILE_OPT idl2
  ENVI, /RESTORE_BASE_SAVE_FILES
  ENVI_BATCH_INIT, /NO_STATUS_WINDOW
  print, '  ' + 'Batch Processing Start'
  print, '  ' + systime()

  ;******************************************************************************** 

  ;Input location of different file
  InputAOD = 'F:\DB_test_data\VAR2\TEST_DATA\Clip\'
  Input_cloud = 'F:\DB_test_data\VAR2\TEST_DATA\Clip\'
  Input_TOA_B1 = 'F:\DB_test_data\VAR2\TEST_DATA\Clip\'
  ;Output location
  Output = 'F:\DB_test_data\OUT_PUT\'

  ;  Search for inpiut data
  FileCloud = File_search(Input_cloud + "*Aerosol_Cldmask_Land_Ocean- 
  Aerosol_Cldmask_Land_Ocean.tif",COUNT=Count_in1)
  FileAOD = File_Search(InputAOD + '*Corrected_Optical_Depth_Land_2- 
  Corrected_Optical_Depth_Land.tif',COUNT=count_in)
  FileBAND1 = File_Search(Input_TOA_B1 + '*pssgrp_000501330352.EV_500_RefSB_1- 
  EV_500_RefSB.tif',COUNT=COUNT)

  nfiles = n_elements(count_in1)

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $$$$$$$$$$$$$$$$$
    FOR k=0, nfiles-1 DO BEGIN
    Print, 'Processing File =', k +1
    Print, 'Total Files=', count_in1
    string_count = strtrim(nfiles,2)

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    ;cloud(500_METER)
    ENVI_OPEN_FILE, FileCloud[k], r_FiD=FiD_cld
    IF (FiD_cld EQ -1) then return
    ENVI_FILE_QUERY, FiD_cld, dims=dims, nb=nb, ns=cld_ncols, nl=cld_nrows
    pos = [0]
    map_info = envi_get_map_info(FiD=FiD_cld)
    Cloud = Double(ENVI_GET_DATA(FiD=FiD_cld, dims=dims, pos=[0]))

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    ;BAND_1(500_METER)
    envi_open_file,FileBAND1[k],r_fid=fid_T_B1
    if (fid_T_B1 eq -1)then return
    envi_file_query,fid_T_B1,dims=dims,nb=nb, ns=toa_B1_ncols, nl=toa_B1_nrows
    pos=[0]
    map_info1=envi_get_map_info(fid=fid_T_B1)
    layer1 = fltarr(toa_B1_ncols,toa_B1_nrows,nb)
    FOR i=0,nb-1 DO BEGIN
      layer1[*,*,i]=ENVI_GET_DATA(fid=fid_T_B1,dims=dims,pos=pos[i])
    ENDFOR
    ;To replace 65535 with NAN
    i=WHERE(layer1 eq 65535,count)
    if (count gt 0)then layer1[i]=!VALUES.F_NAN
    TOA_B1=layer1*0.000053811826d

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    ;AOD(10KM_RESOLUTION)
    ENVI_OPEN_FILE, FileAOD[k], r_FiD=FiD_AOD
    IF (FiD_AOD EQ -1) then return
    ENVI_FILE_QUERY, FiD_AOD, dims=dims, nb=nb, ns=aod_ncols, nl=aod_nrows
    pos = [0]
    map_info = envi_get_map_info(FiD=FiD_AOD)
    layer2 = fltarr(aod_ncols,aod_nrows,nb)
    FOR i=0,nb-1 DO BEGIN
      layer2[*,*,i]=ENVI_GET_DATA(fid=FiD_AOD,dims=dims,pos=pos[i])
    ENDFOR
    ;To replace -9999 with NAN
    i=WHERE(layer2 eq -9999,count)
    if (count gt 0)then layer2[i]=!VALUES.F_NAN
    AOD=layer2*0.0010000000474974513d
    ;set the specific range for mask 
    AOD_MASK = 1.0*(AOD GE 0.0 and AOD LE 0.1) + (AOD GT 0.1)*0.0


    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    ; to get same dimensions for B1TOA and CLD data
    if (cld_ncols le toa_B1_ncols) then begin
      xsize = cld_ncols
    endif else begin
      xsize = toa_B1_ncols
    endelse
    if (cld_nrows le toa_B1_nrows) then begin
      ysize = cld_nrows
    endif else begin
      ysize = toa_B1_nrows
    endelse

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    ; to call parameters again with same dims
    Cloud = FSC_Resize_Image(Cloud,xsize, ysize)
    Band1 = FSC_Resize_Image(TOA_B1,xsize, ysize)

    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
    ;First masking processes with cloud mask images(500 meters)
    Cloud_masked = Cloud * Band1
    ;To replace 0 with NAN
    i=WHERE(Cloud_masked eq 0,count)
    if (count gt 0)then Cloud_masked[i]=!VALUES.F_NAN
    Cloud_masked_Band1=Cloud_masked
    ;Second masking processes with AOD_MASK images(10KM)




    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

    ;  Define output FileName
    FileName = STRJOIN(STRSPLIT(STRMID(FILE_BASENAME(FileBAND1[k], '.tif'), 0, 75), '.', /EXTRACT), '_')
    print, FileName
    ;  write output File 
    out_name_out = Output + FileName +'_masked_image.tif'
    map_info = envi_get_map_info(FiD=fid_T_B1)
    ENVI_WRITE_ENVI_FILE, Cloud_masked, out_name=out_name_out, map_info=map_info


    ;FileName = STRJOIN(STRSPLIT(STRMID(FILE_BASENAME(FileAOD[k], '.tif'), 0, 75), '.', /EXTRACT), '_')
    ;print, FileName
    ;  write output File
    ;out_name_out = Output + FileName +'_AOD.tif'
    ;map_info = envi_get_map_info(FiD=FiD_AOD)
    ;ENVI_WRITE_ENVI_FILE, AOD,out_name=out_name_out, map_info=map_info
    ;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$



    Print, FileName + ' ' + strtrim(k+1,2) + ' of ' + string_count + ' Processed'
  Endfor
  Print, '  ' + systime()
  Print, '  ' +  'Batch Processing End'
End
;$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
FUNCTION FSC_Resize_Image, image, cols, rows, $
  INTERPOLATE=interp, $
  MINUS_ONE=minus_one
  Compile_Opt idl2
  ;     Check parameters.
  IF N_Params() EQ 0 THEN BEGIN
    Print, 'USE SYNTAX: FSC_Resize_Image, image, cols, rows'
    RETURN, "FSC_Resize_Image Syntax Error"
  ENDIF
  ; Return to the caller on an error.
  On_Error, 2
  ;  Check image parameter for size. Only 2D and 3D images allowed.
  ndim = SIZE(image, /N_DIMENSIONS)
  dims = SIZE(image, /DIMENSIONS)
  IF ((ndim LT 2) OR (ndim GT 3)) THEN $
    Message, 'Input image must have two or three dimensions.'
  ;    Check for keywords. Default for minus_one is 0.
  interp = Keyword_Set(interp)
  IF N_Elements(minus_one) EQ 0 THEN minus_one = 0
  m1 = Keyword_Set(minus_one)
  ; 2D images are immediately passed to IDL's Congrid command, with
  ;  the CENTER keyword set.
  IF ndim EQ 2 THEN BEGIN
    RETURN, Congrid(image, cols, rows, CENTER=1, MINUS_ONE=minus_one, $
      INTERP=interp)
  ENDIF
  ;     24-bit images are handled differently. The "3" dimension of a 24-bit
  ;    image should not be interpolated.
  offset = 0.5 ; To center the pixel location (i.e., CENTER=1 for Congrid)
  index3 = Where(dims EQ 3)
  CASE index3 OF
    0: BEGIN
      srx = [0,1,2]
      sry = Float(dims[1] - m1) / (cols - m1) *(Findgen(cols) + offset) - offset
      srz = Float(dims[2] - m1) / (rows - m1) *(Findgen(rows) + offset) - offset
    END
    1: BEGIN
      srx = Float(dims[0] - m1) / (cols - m1) *(Findgen(cols) + offset) - offset
      sry = [0,1,2]
      srz = Float(dims[2] - m1) / (rows - m1) *(Findgen(rows) + offset) - offset
    END
    2: BEGIN
      srx = Float(dims[0] - m1) / (cols - m1) *(Findgen(cols) + offset) - offset
      sry = Float(dims[1] - m1) / (rows - m1) *(Findgen(rows) + offset) - offset
      srz = [0,1,2]
    END
  ENDCASE
  ;     Do the interpolation. Preserve nearest neighbor sampling, if required.
  IF interp THEN BEGIN
    RETURN, Interpolate(image, srx, sry, srz, /GRID)
  ENDIF ELSE BEGIN
    RETURN, Interpolate(image, Round(srx), Round(sry), Round(srz), /GRID)
  ENDELSE
 END`