在OpenTopoMap上貼上NCL等值圖

    44    larg=len(sys.argv)
    45    if larg<6:
    46        fname=sys.argv[1]
    47        with open(fname,'r') as f:
    48            ll=[l for l in f]
    49        ll=ll[8:]
    50        X,Y=([float(l.split()[i]) for l in ll] for i in [0,1])
    51        mnX,mnY,mxX,mxY=(min(X),min(Y),max(X),max(Y))   
    52    else:
    53        mnX,nx,dx,mnY,ny,dy=(float(sys.argv[i]) for i in range(larg-6,larg))
    54        mxX,mxY=mnX+nx*dx,mnY+ny*dy
    55    lat_min,lon_min=TWD2LL(mnX,mnY)
    56    lat_max,lon_max=TWD2LL(mxX,mxY)

    33    zoomi = 22 #initial zoomming level
...
  137    for zoom in range(zoomi,5,-1):
  138      x_min, x_max, y_min, y_max = bbox_to_xyz(
  139        lon_min, lon_max, lat_min, lat_max, zoom)
  140      ntiles=(x_max - x_min + 1) * (y_max - y_min + 1)
  141      if 150 > ntiles > 20:break
  142    print(f"Downloading {(x_max - x_min + 1) * (y_max - y_min + 1)} tiles")

    28    store_dir = os.path.join(os.path.dirname(__file__), '../pngs')
    29    temp_dir = ‘temp'
...
  101    def download_tile(x, y, z, tile_server):
  102        url = tile_server.replace(
  103            "{x}", str(x)).replace(
  104            "{y}", str(y)).replace(
  105            "{z}", str(z))
  106        store = f'{store_dir}/{x}_{y}_{z}.png'
  107        path = f'{temp_dir}/{x}_{y}_{z}.png'
  108        if os.path.exists(store):
  109            os.system('ln -sf '+store+' '+path)   
  110            return(path)
  111        if "street" in tile_server:
  112            os.system(wget+url+' -O '+path)
  113        else:
  114            urllib.request.urlretrieve(url, path)
  115        os.system('cp '+path+' '+store)   
  116        return(path)

    41    ;---Create simple contour plots 
    42      wks_type = "png"
    43      wks_type@wkWidth = 2000
    44      wks_type@wkHeight = 2000
    45      wks = gsn_open_wks(wks_type,"tmp_cn")

#kuang@114-32-164-198 /Users/Data/GIS/OSM_20210318/merged_GeoTIFF
#$ cat -n NCLonOTM.py
    1    #!/opt/anaconda3/envs/env_name/bin/python
...
    11    image=cv2.imread('tmp_cn.png')
    12    ny,nx,nz=image.shape
    13   
    14    #locate the frame of contour plot
    15    nwidth=np.array([len(np.where(image[j,:,:]!=255)[0]) for j in range(ny)])
    16    nheight=np.array([len(np.where(image[:,i,:]!=255)[0]) for i in range(nx)])
    17    wh={'j':'nwidth','i':'nheight'}
    18    nxy={'j':'ny','i':'nx'}
    19    for ij in wh:
    20      n=nxy[ij]
    21      exec('id1=int(np.mean(np.where('+wh[ij]+'>=max('+wh[ij]+')-2)[0]))')
    22      exec('id2=int(np.mean(np.where('+wh[ij]+'==max(['+wh[ij]+'[i] for i in range('+n+') if abs(i-id1) > 2]))))')
    23      if id1==id2:
    24        exec('id1=np.min(np.where('+wh[ij]+'>=max('+wh[ij]+')-2)[0])') 
    25        exec('id2=np.max(np.where('+wh[ij]+'>=max('+wh[ij]+')-2)[0])') 
    26      exec(ij+'0=min([id1,id2])')
    27      exec(ij+'1=max([id1,id2])')
    28    size = (j1-j0,i1-i0)

    30    #darken the lines inside of the frame
    31    image3=image[:]
    32    boo=(image[:,:,0]==image[:,:,1]) & (image[:,:,1]==image[:,:,2])& (image[:,:,2]!=255)
    33    idx=np.where(boo)
    34    idx0=np.array([idx[0][i] for i in range(len(idx[0])) if idx[0][i]>j0 and idx[0][i]<j1 and idx[1][i]>i0 and idx[1][i]<i1])
    35    idx1=np.array([idx[1][i] for i in range(len(idx[0])) if idx[0][i]>j0 and idx[0][i]<j1 and idx[1][i]>i0 and idx[1][i]<i1])
    36    if len(idx0)>0 and len(idx1)>0:image3[idx0,idx1,:]=0

    38    #resize (alignment the resolutions)
    39    image2=cv2.imread('fitted.png')
    40    ny2,nx2,nz=image2.shape
    41    if ny2<size[0] and size[0]/ny2>1.3 :sys.exit('make a more detail basemap')
    42    rate=(ny2/size[0],nx2/size[1])
    43    size2=(int(rate[0]*ny),int(rate[1]*nx))
    44    im_resized = cv2.resize(image3,size2, interpolation=cv2.INTER_CUBIC)
    45    image3=np.array(im_resized)
    46    j0,i0=int(j0*rate[0]),int(i0*rate[1])

    48    #replace the nodata location with OSM plot
    49    for j in range(j0,j0+ny2):
    50      for i in range(i0,i0+nx2):
    51        if (image3[j,i,:]==0).all():continue 
    52    #    if sum(image3[j,i,:])>=200+255:
    53        image3[j,i,:]=image2[j-j0,i-i0,:]
    54    cv2.imwrite("NCLonOTM.png",image3)

    26  def PLT_parser(fname):
    27    with open(fname,'r') as f:
    28      l=[line.strip('\n') for line in f]
    29    if l[0][0]=='*':l=l[8:]
    30    X,Y=([float(l[i].split()[j]) for i in range(len(l))] for j in range(2))
    31    sX ,sY=list(set(X)),list(set(Y))
    32    sX.sort()
    33    sY.sort()
    34    nx,ny=len(sX),len(sY)
    35    dx,dy=[round(sX[i+1]-sX[i],3) for i in range(nx-1)],[round(sY[i+1]-sY[i],3) for i in range(ny-1)]
    36    if len(set(dx))!=1 or len(set(dy))!=1:
    37      fname=pth+fn
    38      print """not a regular RE GRIDCART system, sorry! your input is:
    39      <a data-auto-download href="%s">%s</a>
    40          """  % (fname.replace(WEB,'../../../'),fname.split('/')[-1])
    41      sys.exit('not a regular RE GRIDCART system')
    42    return ' %d %d %d %d %d %d ' % (int(min(X)),nx,int(dx[0]),int(min(Y)),ny,int(dy[0]))

    15  FIT='/Users/Data/GIS/OSM_20210318/merged_GeoTIFF/tiles_to_tiffFit.py '
...
    48  form = cgi.FieldStorage()
    49  STR = str(form.getvalue("iscinp"))
    50  os.system('echo "'+STR+'"'+OUT)
    51  if len(STR)>=4: #in case of input a string
    52    cmd ='cd '+pth+';'
    53    cmd+= FIT+STR+OUT+';'
    54    os.system('echo "'+cmd+'"'+OUT)
    55    r=os.system(cmd+OUT)
    56    if r!=0:sys.exit('error in ncl')
    57    fnames=['fitted.png']

    58  else:
    59    fileitem = form['filename']
    60    if fileitem.filename:
    61      fn = os.path.basename(fileitem.filename)
    62      open(pth+fn, 'wb').write(fileitem.file.read())
    63      with open(pth+fn, 'r') as f:
    64        ll=[l.strip('\n') for l in f]
    65      if ll[0]=='DSAA':
    66        x, y, c, (ny, nx) = load_surfer(pth+fn)
    67        with open(pth+'tmp.PLT','w') as f:
    68          for j in range(ny):
    69            for i in range(nx):
    70              f.write('%f %f %f\n' % (x[j,i],y[j,i],c[j,i]))

    71      elif ll[0].split()[1] in ['AERMOD','ISCST3']:
    72        cmd ='cd '+pth+';'
    73        if pth+fn!=pth+'userinp.PLT': cmd+='cp '+pth+fn+' '+pth+'userinp.PLT;'
    74        cmd+= SED+'userinp.PLT'+CUT+';'
    75        os.system('echo "'+cmd+'"'+OUT)
    76        r=os.system(cmd+OUT)
    77      else:
    78        print 'wrong format! '+ll[0]
    79        sys.exit('wrong format')
    80      cmd ='cd '+pth+';'
    81      cmd+= FIT+' tmp.PLT '+OUT+';'
    82      os.system('echo "'+cmd+'"'+OUT)
    83      r=os.system(cmd+OUT)

    16  OVL='/Users/Data/GIS/OSM_20210318/merged_GeoTIFF/NCLonOTM.py'
...
    85      STR=PLT_parser(pth+'tmp.PLT')
    86      cmd ='cd '+pth+';'
    87      cmd+='echo "'+STR+'">param.txt;'
    88      cmd+='echo "'+fn+'">title.txt;'
    89      cmd+= NCL+OUT+';'
    90      cmd+= OVL
    91      os.system('echo "'+cmd+'"'+OUT)
    92      r=os.system(cmd+OUT)
    93      if r!=0:sys.exit('error in ncl')

    94      fnames=['fitted.png', 'tmp_cn.png', "NCLonOTM.png"]
    95      descip={ 'fitted.png':'OpenTopoMap: ',
    96          'tmp_cn.png':'CLN_contour: ',
    97          'NCLonOTM.png':'contour post on OTM: '}
    98
    99  for fn in fnames:
   100    fname=pth+fn
   101    print """\
   102    %s<a data-auto-download href="%s">%s</a></br>
   103    """  % (descip[fn],fname.replace(WEB,'../../../'),fname.split('/')[-1])
   104  print """\
   105    </body>
   106    </html>
   107    """