PRO SubImage_Example ; Read the image file. file = 'IBCAO_V3_500m_SM.tif' image = Read_Tiff(file) image = Reverse(image, 2) ; Set up the map projection for the image. mapCoord = Obj_New('cgMap', 'Polar Stereographic', Ellipsoid='WGS 84', $ Center_Longitude=0, True_Scale_Latitude=75, $ XRange=[-2904000.d, 2904000.d], YRange=[-2904000.d, 2904000.d]) ; Display the image with map annotations. cgLoadCT, 4, /Brewer, /Reverse cgDisplay, 1000, 1000, Aspect=image cgImage, image, /Keep_Aspect, OPosition=opos, /Save, /Scale, $ MinValue=-4000, MaxValue=4000 mapCoord -> SetProperty, Position=opos cgMap_Continents, Map=mapCoord cgMap_Grid, Map=mapCoord, /Label, Lats=[Findgen(5)*10+40, 85], $ Lons=Findgen(19)*20.0, LonLab=75, LatLab = 0, Color='charcoal' ; Display the area we are interested in. area = mapCoord -> Forward([-80, -80, -30, -30, -80], $ [ 78, 64, 64, 78, 78], /NoForwardFix) cgPlotS, area, Color='dark red', Thick=thick ; Set up vectors that describe the image in XY projected meter space. scale = 500.D ; meters xmin = -2904000.D - (scale/2.0) xmax = 2904000.D + (scale/2.0) ymin = -2904000.D - (scale/2.0) ymax = 2904000.D + (scale/2.0) xvec = cgScaleVector(DIndgen(11617+1), xmin, xmax) yvec = cgScaleVector(DIndgen(11617+1), ymin, ymax) ; Find the pixel locations of the subimage. xpixel = Value_Locate(xvec, [Min(area[0,*]), Max(area[0,*])]) ypixel = Value_Locate(yvec, [Min(area[1,*]), Max(area[1,*])]) Print, 'XPixel Locations: ', xpixel Print, 'YPixel Locations: ', ypixel ; Draw the subimage area. x0 = Min(area[0,*]) x1 = Max(area[0,*]) y0 = Min(area[1,*]) y1 = Max(area[1,*]) cgPlots,[x0,x0,x1,x1,x0], [y0,y1,y1,y0,y0], Color='charcoal', SymSize=2, Thick=2 ; Extract and display the subimage. subImage = image[xpixel[0]:xpixel[1], ypixel[0]:ypixel[1]] cgDisplay, WID=1, 500, 500, Aspect=subimage cgImage, subimage, /Scale, MinValue=-4000, MaxValue=4000 subMapCoord = Obj_New('cgMap', 'Polar Stereographic', Ellipsoid='WGS 84', $ Center_Longitude=0, True_Scale_Latitude=75, Position=[0, 0, 1, 1], $ XRange=[xvec[xpixel]], YRange=[yvec[ypixel]]) cgMap_Continents, Map=subMapCoord cgMap_Grid, Map=subMapCoord, /Label, Lats=Findgen(12)*2+60, $ Lons=Findgen(16)*10-180, LonLab=74, LatLab = -140, Color='charcoal' ; Create a mask for the study area. xrange=[xvec[xpixel]] yrange=[yvec[ypixel]] dims = Size(subimage, /Dimensions) xvector = cgScaleVector(DIndgen(dims[0]), xrange[0], xrange[1]) yvector = cgScaleVector(DIndgen(dims[1]), yrange[0], yrange[1]) area = subMapCoord -> Forward([-80, -80, -30, -30, -80], $ [ 78, 64, 64, 78, 78], /NoForwardFix) mask = BytArr(dims[0], dims[1]) xpixels = Value_Locate(xvector, area[0,*]) ypixels = Value_Locate(yvector, area[1,*]) indices = PolyFillV(xpixels, ypixels, dims[0], dims[1]) mask[indices] = 1 ; Display the subimage with the mask applied. cgDisplay, WID=2, 500, 500, Aspect=subimage scaledImage = BytScl(subImage, Min=-4000, Max=4000, Top=254) + 1B cgLoadCT, 4, /Brewer, /Reverse, NColors=254, Bottom=1 TVLCT, cgColor('white', /Triple), 0 cgImage, scaledImage * mask, XRange=xrange, YRange=yrange cgMap_Continents, Map=subMapCoord cgMap_Grid, Map=subMapCoord, /Label, Lats=Findgen(12)*2+60, $ Lons=Findgen(16)*10-180, LonLab=74, LatLab = -140, Color='charcoal' END