|
| 1 | +#!/usr/bin/env python3 |
| 2 | +import socket |
| 3 | +import copy |
| 4 | +import numpy as np |
| 5 | +import matplotlib.pyplot as plt |
| 6 | +from cmaptools import getcolormap |
| 7 | +from processmesh import processmesh |
| 8 | +from processdata import processdata |
| 9 | +from InterpFromMeshToGrid import InterpFromMeshToGrid |
| 10 | +from InterpFromMeshToMesh2d import InterpFromMeshToMesh2d |
| 11 | +from applyoptions import applyoptions |
| 12 | + |
| 13 | +def plot_landsat(md,data,options,fig,axgrid,gridindex): |
| 14 | + """ |
| 15 | + Explain |
| 16 | + ------- |
| 17 | + This funtion loads Landsat Image Mosaic Antarctica (LIMA) for background image. |
| 18 | +
|
| 19 | + Usage |
| 20 | + ----- |
| 21 | + plot_landsat(md,data,options,plotlines,plotcols,i) |
| 22 | + """ |
| 23 | + |
| 24 | + #process mesh and data |
| 25 | + x2d, y2d, z2d, elements2d, is2d, isplanet=processmesh(md,[],options) |
| 26 | + data, datatype=processdata(md,data,options) |
| 27 | + |
| 28 | + ismask = options.exist('mask') |
| 29 | + if ismask: |
| 30 | + mask = options.getfieldvalue('mask') |
| 31 | + options2 = copy.deepcopy(options) |
| 32 | + options2.removefield('caxis',False) |
| 33 | + options2.removefield('log',False) |
| 34 | + mask, datatype=processdata(md,mask,options2) |
| 35 | + data[~mask] = np.nan |
| 36 | + |
| 37 | + #check is2d |
| 38 | + if not is2d: |
| 39 | + raise Exception('buildgridded error message: gridded not supported for 3d meshes, project on a layer') |
| 40 | + |
| 41 | + #Get options |
| 42 | + hostname = socket.gethostname().lower().replace('-','') |
| 43 | + transparency = options.getfieldvalue('transparency',.2) |
| 44 | + highres = options.getfieldvalue('highres',0) |
| 45 | + isunit = options.getfieldvalue('unit',1) |
| 46 | + |
| 47 | + #Get xlim, and ylim |
| 48 | + xlim=options.getfieldvalue('xlim',np.array([min(x2d),max(x2d)]))/isunit |
| 49 | + ylim=options.getfieldvalue('ylim',np.array([min(y2d),max(y2d)]))/isunit |
| 50 | + |
| 51 | + pwr = md.radaroverlay.pwr |
| 52 | + xm = md.radaroverlay.x |
| 53 | + ym = md.radaroverlay.y |
| 54 | + if md.mesh.epsg == 3031 & np.size(pwr)==0 | np.size(xm)==0 | np.size(ym) == 0: |
| 55 | + #Antarctica region |
| 56 | + if highres: |
| 57 | + print(' LIMA with geotiff') # {{{ |
| 58 | + |
| 59 | + # find merged mosaic landsat image |
| 60 | + limapath = {'simba00':'/drive/project_inwoo/issm/Data/LIMA/AntarcticaLandsat.tif'}; |
| 61 | + limapath = limapath[hostname] |
| 62 | + print(' LIMA path is %s'%(limapath)) |
| 63 | + |
| 64 | + # read image |
| 65 | + #im = imread(limapath); |
| 66 | + |
| 67 | + ## Region of LIMA data set |
| 68 | + #info = gdalinfo(limapath); # get geotiff info |
| 69 | + #xm = info.xmin + info.dx*np.arange(info.nx) |
| 70 | + #ym = info.ymax - info.dy*np.arange(info.ny) |
| 71 | + |
| 72 | + ## find region of model at LIMA |
| 73 | + #offset = 1e+4; |
| 74 | + #posx = np.where((xm > xlim[0]-offset)&(xm < xlim[1]+offset))[0] |
| 75 | + #posy = np.where((ym > ylim[0]-offset)&(ym < ylim[1]+offset))[0] |
| 76 | + # }}} |
| 77 | + else: |
| 78 | + print(' LIMA with reduced tiff') # {{{ |
| 79 | + #Find merged mosaic landsat image |
| 80 | + limapath = {'inwoob85md3h':'/drive/project_inwoo/issm/Data/LIMA/tiff_90pct/00000-20080319-092059124.tif', |
| 81 | + 'simba00':'/home/inwoo/data/LIMA/tiff_90pct/00000-20080319-092059124.tif'} |
| 82 | + if not hostname in limapath.keys(): |
| 83 | + raise Exception('Error: Landsat image at Antarctic region is downloaded at https://lima.usgs.gov/fullcontinent.php. Download geotiff image using "wget -c https://lima.usgs.gov/tiff_90pct.zip -O tiff_90pct.zip"'); |
| 84 | + |
| 85 | + limapath = limapath[hostname] |
| 86 | + print(' LIMA path is %s'%(limapath)) |
| 87 | + |
| 88 | + # read image |
| 89 | + #im = imread(limapath) |
| 90 | + |
| 91 | + ## Region of LIMA data set |
| 92 | + #info = gdalinfo(limapath) # get geotiff info |
| 93 | + #xm = info.xmin + info.dx*np.arange(info.nx) |
| 94 | + #ym = info.ymax - info.dy*np.arange(info.ny) |
| 95 | + |
| 96 | + ## find region of model at LIMA |
| 97 | + #offset = 1e+4 |
| 98 | + #posx = np.where((xm > xlim[0]-offset)&(xm < xlim[1]+offset))[0] |
| 99 | + #posy = np.where((ym > ylim[0]-offset)&(ym < ylim[1]+offset))[0] |
| 100 | + # }}} |
| 101 | + |
| 102 | + # Update region of radaroverlay |
| 103 | + md.radaroverlay.x = xm[posx] |
| 104 | + md.radaroverlay.y = ym[posy] |
| 105 | + md.radaroverlay.pwr = im[posy, posx,:] |
| 106 | + elif md.mesh.epsg == 3431 & np.size(pwr)==0 | np.size(xm)==0 | np.size(ym) == 0: |
| 107 | + #Greenland region |
| 108 | + raise Exception('Greenland region is not yet available.') |
| 109 | + |
| 110 | + #Check dataset. |
| 111 | + if (np.size(pwr)>0) & (np.size(xm)>0) & (np.size(ym)>0): |
| 112 | + #Existing radaroverlay |
| 113 | + if np.ndim(pwr) != 3: |
| 114 | + raise Exception('Error: Check np.ndim(md.radaroverlay.pwr) should be equal to 3.') |
| 115 | + |
| 116 | + #Check image size |
| 117 | + #shape of image should be (ny, nx, band) |
| 118 | + nx = len(xm) |
| 119 | + ny = len(ym) |
| 120 | + if (np.shape(pwr)[0]==nx) & (np.shape(pwr)[1]==ny): |
| 121 | + pwr = np.transpose(pwr,[1,0,2]) |
| 122 | + |
| 123 | + #Close-up to specific xlim |
| 124 | + posx = np.where((xlim[0]<=xm)&(xm<=xlim[1]))[0] |
| 125 | + posy = np.where((ylim[0]<=ym)&(ym<=ylim[1]))[0] |
| 126 | + xm = xm[posx] |
| 127 | + ym = ym[posy] |
| 128 | + pwr = pwr[posy[0]:posy[-1]+1,posx[0]:posx[-1]+1,:] |
| 129 | + else: |
| 130 | + raise Exception('Error: data array in md.radaroverlay is not implemented yet.') |
| 131 | + |
| 132 | + #Prepare grid |
| 133 | + if np.ndim(xm) == 1: |
| 134 | + data_grid=InterpFromMeshToGrid(elements2d+1,x2d/isunit,y2d/isunit,data,xm,ym,np.nan) |
| 135 | + else: |
| 136 | + data_grid=InterpFromMeshToMesh2d(elements2d+1,x2d,y2d,data,np.ravel(xm),np.ravel(ym),'default',np.nan) |
| 137 | + data_grid=np.reshape(data_grid,np.shape(xm)) |
| 138 | + |
| 139 | + data_nan=np.isnan(data_grid) |
| 140 | + if options.exist('caxis'): |
| 141 | + caxis_opt=options.getfieldvalue('caxis') |
| 142 | + data_grid[np.where(data_grid<caxis_opt[0])]=caxis_opt[0] |
| 143 | + data_grid[np.where(data_grid>caxis_opt[1])]=caxis_opt[1] |
| 144 | + data_min=caxis_opt[0]; |
| 145 | + data_max=caxis_opt[1]; |
| 146 | + else: |
| 147 | + data_min=np.nanmin(data_grid) |
| 148 | + data_max=np.nanmax(data_grid) |
| 149 | + |
| 150 | + options.addfielddefault('colormap',plt.cm.viridis) |
| 151 | + cmap = getcolormap(copy.deepcopy(options)) |
| 152 | + #TODO: Matlab version |
| 153 | + #image_rgb = ind2rgb(uint16((data_grid - data_min)*(length(colorm)/(data_max-data_min))),colorm); |
| 154 | + #NOTE: Python version for ind2rgb |
| 155 | + image_rgb = cmap((data_grid-data_min)/(data_max-data_min)) |
| 156 | + |
| 157 | + alpha=np.ones(np.shape(data_grid)) |
| 158 | + alpha[np.where(~data_nan)]=transparency |
| 159 | + alpha=np.repeat(alpha[:,:,np.newaxis],axis=2,repeats=3) |
| 160 | + |
| 161 | + final=alpha*(pwr/255)+(1-alpha)*image_rgb[:,:,:3] |
| 162 | + |
| 163 | + #Select plot area |
| 164 | + ax = axgrid[gridindex] |
| 165 | + |
| 166 | + xmin = min(xm)/isunit |
| 167 | + xmax = max(xm)/isunit |
| 168 | + ymin = min(ym)/isunit |
| 169 | + ymax = max(ym)/isunit |
| 170 | + #Draw RGB image |
| 171 | + h=ax.imshow(final, extent=[xmin, xmax, ymin, ymax], origin='lower') |
| 172 | + |
| 173 | + #last step: mesh gridded? |
| 174 | + if options.exist('edgecolor'): |
| 175 | + ax.triplot(x,y,triangles=elements, |
| 176 | + color=options.getfieldvalue('edgecolor'), |
| 177 | + linewdith=options.getfieldvalue('linewidth',1), |
| 178 | + ) |
| 179 | + |
| 180 | + #Apply options |
| 181 | + if ~np.isnan(data_min): |
| 182 | + options.changefieldvalue('caxis',[data_min, data_max]) # force caxis so that the colorbar is ready |
| 183 | + options.addfielddefault('axis','xy equal off') # default axis |
| 184 | + applyoptions(md,data,options,fig,axgrid,gridindex) |
0 commit comments