Python Basemap.projtran方法代码示例

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
def ss_plot():
    #Pandas method of importing data frame and getting extents
    conn = psycopg2.connect("dbname='reach_4a' user='root'  host='localhost' port='9000'")
    df = pd.read_sql_query('select easting, northing, texture, sidescan_intensity from mosaic_2014_09', con=conn)
    minE = min(df['easting'])
    maxE = max(df['easting'])
    minN = min(df['northing'])
    maxN = max(df['northing'])
    conn.close()
    print 'Done Importing Data from Database'
    
    #Create grid for countourf plot
    res = 0.25
    grid_x, grid_y = np.meshgrid( np.arange(np.floor(minE), np.ceil(maxE), res), np.arange(np.floor(minN), np.ceil(maxN), res))
    grid_lon, grid_lat = trans(grid_x,grid_y,inverse=True)
    
    #Re-sampling procedure
    m_lon, m_lat = trans(df['easting'].values.flatten(), df['northing'].values.flatten(), inverse=True)
    orig_def = geometry.SwathDefinition(lons=m_lon, lats=m_lat)
    target_def = geometry.SwathDefinition(lons=grid_lon.flatten(), lats=grid_lat.flatten())
    print 'Now Resampling...'
    result = kd_tree.resample_nearest(orig_def, df['sidescan_intensity'].values.flatten(), target_def, radius_of_influence=1, fill_value=None, nprocs = cpu_count())
    print 'Done Resampling!!!' 
    
    #format side scan intensities grid for plotting
    gridded_result = np.reshape(result,np.shape(grid_lon))
    gridded_result = np.squeeze(gridded_result)
    gridded_result[np.isinf(gridded_result)] = np.nan
    gridded_result[gridded_result<=0] = np.nan
    grid2plot = np.ma.masked_invalid(gridded_result)
    
#    x = df['easting']
#    y = df['northing']
#    z = df['sidescan_intensity']
#    
#    xi = df['easting']
#    yi = df['northing']
#    
#    X,Y= np.meshgrid(xi,yi)
#    grid_lon, grid_lat = trans(X,Y,inverse=True)
#    Z = griddata((x, y), z, (X, Y),method='nearest')
#    print 'Done Gridding Data'
    
    
    print 'Now mapping...'
    #Create Figure
    fig = plt.figure(frameon=True)
    ax = plt.subplot(1,1,1)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], llcrnrlon=np.min(grid_lon)-0.0009, llcrnrlat=np.min(grid_lat)-0.0009,urcrnrlon=np.max(grid_lon)+0.0009, urcrnrlat=np.max(grid_lat)+0.0009)
    gx,gy = map.projtran(grid_lon,grid_lat)
    map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=1200)
    im = map.pcolormesh(gx, gy, grid2plot, cmap='gray',vmin=0.1, vmax=30)
        
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbr = plt.colorbar(im, cax=cax)
    cbr.set_label('Sidescan Intensity [dBw]', size=8)
    for t in cbr.ax.get_yticklabels():
        t.set_fontsize(8)
    plt.show()    

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
def ss_plot():
    #Pandas method of importing data frame and getting extents
    db_connect="dbname='reach_4a' user='root'  host='localhost' port='9000'"
    conn = psycopg2.connect(db_connect)
    df = pd.read_sql_query('SELECT * from mb_may_2012_1m tt inner join (	SELECT s.easting, s.northing, s.texture, s.sidescan_intensity  FROM ss_2012_05 s) ss on tt.easting=ss.easting and tt.northing=ss.northing;', con=conn)
    minE = df['easting'].min()[0]
    maxE = df['easting'].max()[0]
    minN = df['northing'].min()[0]
    maxN = df['northing'].max()[0]
    conn.close()
    print 'Done Importing Data from Database'
    
    #Create grid for countourf plot
    res = 1
    grid_x, grid_y = np.meshgrid( np.arange(np.floor(minE), np.ceil(maxE), res), np.arange(np.floor(minN), np.ceil(maxN), res))
    grid_lon, grid_lat = trans(grid_x,grid_y,inverse=True)
    
    #Re-sampling procedure
    m_lon, m_lat = trans(df['easting'].values.flatten(), df['northing'].values.flatten(), inverse=True)
    orig_def = geometry.SwathDefinition(lons=m_lon, lats=m_lat)
    target_def = geometry.SwathDefinition(lons=grid_lon.flatten(), lats=grid_lat.flatten())
    print 'Now Resampling...'
    result = kd_tree.resample_nearest(orig_def, df['sidescan_intensity'].values.flatten(), target_def, radius_of_influence=1, fill_value=None, nprocs = cpu_count())
    print 'Done Resampling!!!' 
    
    #format side scan intensities grid for plotting
    gridded_result = np.reshape(result,np.shape(grid_lon))
    gridded_result = np.squeeze(gridded_result)
    gridded_result[np.isinf(gridded_result)] = np.nan
    gridded_result[gridded_result<=0] = np.nan
    grid2plot = np.ma.masked_invalid(gridded_result)
       
    
    print 'Now mapping...'
    #Create Figure
    fig = plt.figure(frameon=True)
    ax = plt.subplot(1,1,1)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], llcrnrlon=np.min(grid_lon)-0.0009, llcrnrlat=np.min(grid_lat)-0.0009,urcrnrlon=np.max(grid_lon)+0.0009, urcrnrlat=np.max(grid_lat)+0.0009)
    gx,gy = map.projtran(grid_lon,grid_lat)
    map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=1200)
    im = map.pcolormesh(gx, gy, grid2plot, cmap='gray',vmin=0.1, vmax=30)
        
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.1)
    cbr = plt.colorbar(im, cax=cax)
    cbr.set_label('Sidescan Intensity [dBw]', size=8)
    for t in cbr.ax.get_yticklabels():
        t.set_fontsize(8)
    plt.savefig(r'C:\workspace\Texture_Classification\output\May2012_1m_sidescan_intensity.png')    

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
def print_contour_map(cs2cs_args, humlon, humlat, glon, glat, datm, sonpath, p, vmin, vmax): #merge, 

    #levels = [0,0.25,0.5,0.75,1.25,1.5,1.75,2,3,5]
          
    print("drawing and printing map ...")
    fig = plt.figure(frameon=False)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #26949,
     resolution = 'i', #h #f
     llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
     urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

    #if dogrid==1:
    gx,gy = map.projtran(glon, glat)

    ax = plt.Axes(fig, [0., 0., 1., 1.], )
    ax.set_axis_off()
    fig.add_axes(ax)

    try:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
    except:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1000, ypixels=None, dpi=300)
    #finally:
    #   print "error: map could not be created..."
             
    #if dogrid==1:
    if 2>1:
       if datm.size > 25000000:
          print("matrix size > 25,000,000 - decimating by factor of 5 for display")
          #map.contourf(gx[::5,::5], gy[::5,::5], datm[::5,::5], levels, cmap='YlOrRd')
          map.pcolormesh(gx, gy, datm, cmap='YlOrRd', vmin=vmin, vmax=vmax)
       else:
          #map.contourf(gx, gy, datm, levels, cmap='YlOrRd')
          map.pcolormesh(gx[::5,::5], gy[::5,::5], datm[::5,::5], cmap='pink', vmin=vmin, vmax=vmax)
                   
    custom_save2(sonpath,'class_map_imagery'+str(p))
    del fig 

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]

#.........这里部分代码省略.........
                plt.xlim(1,8); plt.ylim(1,8)
                #plt.show()
                custom_save(sonpath,'e1e2_kmeans'+str(p))
                del fig

             except:
                print "plot could not be produced"

          if doplot==1:
             try:

                fig = plt.figure()
                s=plt.scatter(Zes[labels==0],Zns[labels==0],marker='o',c='k', s=10, linewidth=0, vmin=0, vmax=8);
                s=plt.scatter(Zes[labels==1],Zns[labels==1],marker='o',c='r', s=10, linewidth=0, vmin=0, vmax=8);
                s=plt.scatter(Zes[labels==2],Zns[labels==2],marker='o',c='b', s=10, linewidth=0, vmin=0, vmax=8);
                custom_save(sonpath,'rgh_hard_kmeans'+str(p))
                del fig

             except:
                print "plot could not be produced"

          if doplot==1:
             try:

                print "drawing and printing map ..."
                fig = plt.figure(frameon=False)
                #fig.subplots_adjust(wspace = 0.4, hspace=0.4)
                map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #epsg=26949,
                   resolution = 'i', #h #f
                   llcrnrlon=np.min(Zlon)-0.0001, llcrnrlat=np.min(Zlat)-0.0001,
                   urcrnrlon=np.max(Zlon)+0.0001, urcrnrlat=np.max(Zlat)+0.0001)

                # draw point cloud
                x,y = map.projtran(Zlon, Zlat)

                cs = map.scatter(x.flatten(), y.flatten(), 1, rough.flatten(), linewidth=0, vmin=0, vmax=8)

                map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
  
                cbar = map.colorbar(cs,location='bottom',pad="5%")
                cbar.set_label('E1')
                cbar.set_ticks([0,2,4,6,8])

                custom_save(sonpath,'map_rgh'+str(p))
                del fig 

             except:
                print "plot could not be produced"

          if doplot==1:
             try:
                fig = plt.figure()
                #fig.subplots_adjust(wspace = 0.4, hspace=0.4)
                map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1],
                   resolution = 'i', #h #f
                   llcrnrlon=np.min(Zlon)-0.0001, llcrnrlat=np.min(Zlat)-0.0001,
                   urcrnrlon=np.max(Zlon)+0.0001, urcrnrlat=np.max(Zlat)+0.0001)

                # draw point cloud
                x,y = map.projtran(Zlon, Zlat)

                cs = map.scatter(x.flatten(), y.flatten(), 1, hard.flatten(), linewidth=0, vmin=0, vmax=8)

                map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)

                cbar = map.colorbar(cs,location='bottom',pad="5%")

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]

#.........这里部分代码省略.........
      geotransform = (xmin, xres, 0, ymax, 0, -yres)

      ds.SetGeoTransform(geotransform)
      ss_band = ds.GetRasterBand(1)
      ss_band.WriteArray(np.flipud(datout)) #datout)
      ss_band.SetNoDataValue(-99)
      ss_band.FlushCache()
      ss_band.ComputeStatistics(False)
      del ds   
   
   except:
      print("error: geotiff could not be created... check your gdal/ogr install")


   try:

      # =========================================================
      print("creating kmz file ...")
      ## new way to create kml file
      pixels = 1024 * 10

      fig, ax = humutils.gearth_fig(llcrnrlon=glon.min(),
                     llcrnrlat=glat.min(),
                     urcrnrlon=glon.max(),
                     urcrnrlat=glat.max(),
                     pixels=pixels)
      cs = ax.pcolormesh(glon, glat, datm, vmax=scalemax, cmap='gray')
      ax.set_axis_off()
      fig.savefig(os.path.normpath(os.path.join(sonpath,'map'+str(p)+'.png')), transparent=True, format='png')
      del fig, ax

      # =========================================================
      fig = plt.figure(figsize=(1.0, 4.0), facecolor=None, frameon=False)
      ax = fig.add_axes([0.0, 0.05, 0.2, 0.9])
      cb = fig.colorbar(cs, cax=ax)
      cb.set_label('Intensity [dB W]', rotation=-90, color='k', labelpad=20)
      fig.savefig(os.path.normpath(os.path.join(sonpath,'legend'+str(p)+'.png')), transparent=False, format='png')
      del fig, ax, cs, cb

      # =========================================================
      humutils.make_kml(llcrnrlon=glon.min(), llcrnrlat=glat.min(),
         urcrnrlon=glon.max(), urcrnrlat=glat.max(),
         figs=[os.path.normpath(os.path.join(sonpath,'map'+str(p)+'.png'))],
         colorbar=os.path.normpath(os.path.join(sonpath,'legend'+str(p)+'.png')),
         kmzfile=os.path.normpath(os.path.join(sonpath,'GroundOverlay'+str(p)+'.kmz')),
         name='Sidescan Intensity')

   except:
      print("error: map could not be created...")


   #y1 = np.min(glat)-0.001
   #x1 = np.min(glon)-0.001
   #y2 = np.max(glat)+0.001
   #x2 = np.max(glon)+0.001

   print("drawing and printing map ...")
   fig = plt.figure(frameon=False)
   map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1],
    resolution = 'i', #h #f
    llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(glat)-0.001,
    urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(glat)+0.001)

   try:
      map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
   except:
      map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1000, ypixels=None, dpi=300)
   #finally:
   #   print "error: map could not be created..."

   #if dogrid==1:
   gx,gy = map.projtran(glon, glat)

   ax = plt.Axes(fig, [0., 0., 1., 1.], )
   ax.set_axis_off()
   fig.add_axes(ax)

   #if dogrid==1:
   if 2>1:
      if datm.size > 25000000:
         print("matrix size > 25,000,000 - decimating by factor of 5 for display")
         map.pcolormesh(gx[::5,::5], gy[::5,::5], datm[::5,::5], cmap='gray', vmin=np.nanmin(datm), vmax=scalemax) #vmax=np.nanmax(datm)
      else:
         map.pcolormesh(gx, gy, datm, cmap='gray', vmin=np.nanmin(datm), vmax=scalemax) #vmax=np.nanmax(datm)
      del datm, dat
   else:
      ## draw point cloud
      x,y = map.projtran(humlon, humlat)
      map.scatter(x.flatten(), y.flatten(), 0.5, merge.flatten(), cmap='gray', linewidth = '0')

   #map.drawmapscale(x1+0.001, y1+0.001, x1, y1, 200., units='m', barstyle='fancy', labelstyle='simple', fontcolor='k') #'#F8F8FF')
   #map.drawparallels(np.arange(y1-0.001, y2+0.001, 0.005),labels=[1,0,0,1], linewidth=0.0, rotation=30, fontsize=8)
   #map.drawmeridians(np.arange(x1, x2, 0.002),labels=[1,0,0,1], linewidth=0.0, rotation=30, fontsize=8)

   custom_save2(sonpath,'map_imagery'+str(p))
   del fig


   del humlat, humlon
   return res #return the new resolution

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]

#.........这里部分代码省略.........
    else:
       if weight < 4:
          Sdat_g = (np.nansum(w * S.flatten()[inds], axis=1) / np.nansum(w, axis=1)).reshape(shape)
       else:
          Sdat_g = (np.nansum(S.flatten()[inds], axis=1)).reshape(shape)
       del w
       dist = np.nanmean(dist,axis=1).reshape(shape)

    del S

    Sdat_g[dist>1] = np.nan
    Sdat_g[Sdat_g<noisefloor] = np.nan

    dat = Sdat_g.copy()
    dat[dist>1] = 0
    dat2 = replace_nans.RN(dat.astype('float64'),1000,0.01,2,'localmean').getdata()
    dat2[dat==0] = np.nan
    del dat

    dat2[dat2<noisefloor] = np.nan

    Sdat_g = dat2.copy()
    del dat2
   
    Sdat_g[Sdat_g==0] = np.nan
    Sdat_g[np.isinf(Sdat_g)] = np.nan
    Sdat_gm = np.ma.masked_invalid(Sdat_g)
    del Sdat_g

    glon, glat = trans(grid_x, grid_y, inverse=True)
    del grid_x, grid_y
    
    # =========================================================
    print("creating kmz file ...")
    ## new way to create kml file  
    pixels = 1024 * 10
 
    fig, ax = humutils.gearth_fig(llcrnrlon=glon.min(),
                     llcrnrlat=glat.min(),
                     urcrnrlon=glon.max(),
                     urcrnrlat=glat.max(),
                     pixels=pixels)
    cs = ax.pcolormesh(glon, glat, Sdat_gm)
    ax.set_axis_off()
    fig.savefig(os.path.normpath(os.path.join(sonpath,'class_overlay1.png')), transparent=True, format='png')    
    

    fig = plt.figure(figsize=(1.0, 4.0), facecolor=None, frameon=False)
    ax = fig.add_axes([0.0, 0.05, 0.2, 0.9])
    cb = fig.colorbar(cs, cax=ax)
    cb.set_label('Texture lengthscale [m]', rotation=-90, color='k', labelpad=20)
    fig.savefig(os.path.normpath(os.path.join(sonpath,'class_legend.png')), transparent=False, format='png')  


    humutils.make_kml(llcrnrlon=glon.min(), llcrnrlat=glat.min(),
         urcrnrlon=glon.max(), urcrnrlat=glat.max(),
         figs=[os.path.normpath(os.path.join(sonpath,'class_overlay1.png'))], 
         colorbar=os.path.normpath(os.path.join(sonpath,'class_legend.png')),
         kmzfile=os.path.normpath(os.path.join(sonpath,'class_GroundOverlay.kmz')), 
         name='Sidescan Intensity')


    # =========================================================
    print("drawing and printing map ...")
    fig = plt.figure(frameon=False)
    map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], 
     resolution = 'i', #h #f
     llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
     urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

    gx,gy = map.projtran(glon, glat)
       
    try:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='ESRI_Imagery_World_2D', xpixels=1000, ypixels=None, dpi=300)
    except:
       map.arcgisimage(server='http://server.arcgisonline.com/ArcGIS', service='World_Imagery', xpixels=1000, ypixels=None, dpi=300)
    #finally:
    #   print "error: map could not be created..."
      
    ax = plt.Axes(fig, [0., 0., 1., 1.], )
    ax.set_axis_off()
    fig.add_axes(ax)

    if Sdat_gm.size > 25000000:
       print("matrix size > 25,000,000 - decimating by factor of 5 for display")
       map.pcolormesh(gx[::5,::5], gy[::5,::5], Sdat_gm[::5,::5], vmin=np.nanmin(Sdat_gm), vmax=np.nanmax(Sdat_gm))
    else:
       map.pcolormesh(gx, gy, Sdat_gm, vmin=np.nanmin(Sdat_gm), vmax=np.nanmax(Sdat_gm))

    custom_save2(sonpath,'class_map_imagery')
    del fig 

   
    if os.name=='posix': # true if linux/mac
       elapsed = (time.time() - start)
    else: # windows
       elapsed = (time.clock() - start)
    print("Processing took "+str(elapsed)+"seconds to analyse")

    print("Done!")

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
tex_levels = list(np.arange(0,135,5))

ss_level=[0,2.5,5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35]
       
fig = plt.figure(figsize=(15,6))
ax = fig.add_subplot(1,5,1)
ax.set_title('50 square pixel')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.min(glon-0.0009), 
            llcrnrlat=np.min(glat-0.0006),
            urcrnrlon=np.max(glon+0.0009), 
            urcrnrlat=np.max(glat+0.0009))
m.wmsimage(server='http://grandcanyon.usgs.gov/arcgis/services/Imagery/ColoradoRiverImageryExplorer/MapServer/WmsServer?', layers=['0'], xpixels=1000)
x,y = m.projtran(glon, glat)
im = m.contourf(x,y,ss_data_50.T, cmap='Greys_r',levels=ss_level)
im2 = m.contourf(x, y, tex_data_50.T, alpha=0.4, cmap='YlOrRd', levels=tex_levels)#levels=tex_levels
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cax2 = divider.append_axes("right", size="5%", pad=0.3)
cbr = plt.colorbar(im, cax=cax)
cbr2 = plt.colorbar(im2,cax=cax2)



ax1 = fig.add_subplot(1,5,2)
ax1.set_title('70 square pixel')
m1 = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.min(glon-0.0009), 

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
glon, glat = trans(xx, yy, inverse=True)
#ortho_lon, ortho_lat = trans(ortho_x, ortho_y, inverse=True)
cs2cs_args = "epsg:26949"
ss_level=[0,2.5,5,7.5,10,12.5,15,17.5,20,22.5,25,27.5,30,32.5,35]

fig = plt.figure(figsize=(15,12))
ax = plt.subplot2grid((5,2),(0, 0),rowspan=4)
ax.set_title('September 2014 \n R01767')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.min(glon) - 0.0002, 
            llcrnrlat=np.min(glat) - 0.0006,
            urcrnrlon=np.max(glon) + 0.0002, 
            urcrnrlat=np.max(glat) + 0.0006)
m.wmsimage(server='http://grandcanyon.usgs.gov/arcgis/services/Imagery/ColoradoRiverImageryExplorer/MapServer/WmsServer?', layers=['3'], xpixels=1000)
x,y = m.projtran(glon, glat)
im = m.contourf(x,y,data.T, cmap='Greys_r',levels=ss_level)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
cbr = plt.colorbar(im, cax=cax)
cbr.set_label('Sidescan Intensity [dBW]', size=10)

#read shapefile and create polygon collections
##NOTE: Shapefile has to be in WGS84
m.readshapefile( r"C:\workspace\Merged_SS\window_analysis\shapefiles\tex_seg_2014_09_67_geo","layer",drawbounds = False)

#sand, sand/gravel, gravel/sand, ledge, gravel, gravel/boulders, boulders, boulder
s_patch, sg_patch, gs_patch, g_patch, gb_patch, b_patch  =[],[],[],[],[],[]

for info, shape in zip(m.layer_info, m.layer):
    if info['substrate'] == 'sand':

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
print 'Now plotting R02028 Acoutic sediment classifications...'
#Begin the plot
cs2cs_args = "epsg:26949"
fig = plt.figure(figsize=(15,12))
ax = plt.subplot2grid((5,2),(0, 0),rowspan=4)
ax.set_title('R02028 \n May 2014 Acousic Sediment Classifications')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.nanmin(r28_lon)-0.0009, 
            llcrnrlat=np.nanmin(r28_lat)-0.0006,
            urcrnrlon=np.nanmax(r28_lon)+0.0009, 
            urcrnrlat=np.nanmax(r28_lat)+0.0006)
#m.wmsimage(server='http://grandcanyon.usgs.gov/arcgis/services/Imagery/ColoradoRiverImageryExplorer/MapServer/WmsServer?', layers=['3'], xpixels=1000)

x,y = m.projtran(r28_lon, r28_lat)
im = m.contourf(x,y,R02028.T, cmap='Greys_r',levels=ss_level)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbr = plt.colorbar(im, cax=cax)
m.readshapefile(r"C:\workspace\Reach_4a\Multibeam\mb_sed_class\output\shapefiles\may2014_3m_buff_geo","layer",drawbounds = False)

#sand, sand/gravel, gravel, sand/rock, rock
s_patch, sg_patch, g_patch, sr_patch, r_patch, = [],[],[],[],[]

bound = max(stat['count'] for stat in stats_28)/2
for info, shape in zip(m.layer_info, m.layer):
    if info['substrate'] == 'sand' and info['count_28'] > bound:
        s_patch.append(Polygon(np.asarray(shape),True))
    if info['substrate'] == 'sand/gravel' and info['count_28'] > bound:
        sg_patch.append(Polygon(np.asarray(shape),True))        

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
circ5 = Line2D([0], [0], linestyle="none", marker="o", markersize=10, markerfacecolor=colors[4],alpha=a_val)

print 'Now plotting R02028 Acoutic sediment classifications...'
#Begin the plot
cs2cs_args = "epsg:26949"
fig = plt.figure(figsize=(15,12))
ax = plt.subplot2grid((5,2),(0, 0),rowspan=4)
ax.set_title('August 2013 Acousic \n Sediment Classifications')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.nanmin(glon)-0.0009, 
            llcrnrlat=np.nanmin(glat)-0.0006,
            urcrnrlon=np.nanmax(glon)+0.0009, 
            urcrnrlat=np.nanmax(glat)+0.0006)
m.wmsimage(server='http://grandcanyon.usgs.gov/arcgis/services/Imagery/ColoradoRiverImageryExplorer/MapServer/WmsServer?', layers=['3'], xpixels=1000)
x,y = m.projtran(aug_13_lon, aug_13_lat)
im = m.contourf(x,y,aug_sed_class.T, cmap='coolwarm', levels=[0,1,2,3,4,5])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbr = plt.colorbar(im, cax=cax)
ax.legend((circ1, circ2, circ3,circ4,circ5),('rock','sand/rock','Gravel','Sand/Gravel','sand'),numpoints=1, loc='best')

print 'Now plotting May 2014 Acoustic Sediment Classifications...'
ax = plt.subplot2grid((5,2),(0, 1),rowspan=4)
ax.set_title('May 2014 Acousic \n Sediment Classifications')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.nanmin(may_lon)-0.0009, 
            llcrnrlat=np.nanmin(may_lat)-0.0006,
            urcrnrlon=np.nanmax(may_lon)+0.0009, 
            urcrnrlat=np.nanmax(may_lat)+0.0006)

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
def make_map(e, n, t, d, dat_port, dat_star, pix_m, res, cs2cs_args, sonpath, p, dogrid):
   
   trans =  pyproj.Proj(init=cs2cs_args)   
   
   merge = np.vstack((dat_port,dat_star))
   #merge = np.vstack((np.flipud(port_fp[p]),star_fp[p]))
   
   merge[np.isnan(merge)] = 0

   merge = merge[:,:len(n)]

   # get number pixels in scan line
   extent = int(np.shape(merge)[0]/2)

   yvec = np.linspace(pix_m,extent*pix_m,extent)

   print "getting point cloud ..."
   # get the points by rotating the [x,y] vector so it lines up with boat heading, assumed to be the same as the curvature of the [e,n] trace
   X=[]; Y=[];
   for k in range(len(n)): 
      x = np.concatenate((np.tile(e[k],extent) , np.tile(e[k],extent)))
      #y = np.concatenate((n[k]+yvec, n[k]-yvec))
      rangedist = np.sqrt(np.power(yvec, 2.0) - np.power(d[k], 2.0))
      y = np.concatenate((n[k]+rangedist, n[k]-rangedist))
      # Rotate line around center point
      xx = e[k] - ((x - e[k]) * np.cos(t[k])) - ((y - n[k]) * np.sin(t[k]))
      yy = n[k] - ((x - e[k]) * np.sin(t[k])) + ((y - n[k]) * np.cos(t[k]))
      xx, yy = calc_beam_pos(d[k], t[k], xx, yy)
      X.append(xx)
      Y.append(yy) 

   del e, n, t, x, y #, X, Y

   # merge flatten and stack
   X = np.asarray(X,'float').T
   X = X.flatten()

   # merge flatten and stack
   Y = np.asarray(Y,'float').T
   Y = Y.flatten()

   X = X[np.where(np.logical_not(np.isnan(Y)))]
   merge = merge.flatten()[np.where(np.logical_not(np.isnan(Y)))]
   Y = Y[np.where(np.logical_not(np.isnan(Y)))]

   Y = Y[np.where(np.logical_not(np.isnan(X)))]
   merge = merge.flatten()[np.where(np.logical_not(np.isnan(X)))]
   X = X[np.where(np.logical_not(np.isnan(X)))]


   X = X[np.where(np.logical_not(np.isnan(merge)))]
   Y = Y[np.where(np.logical_not(np.isnan(merge)))]
   merge = merge[np.where(np.logical_not(np.isnan(merge)))]

   # write raw bs to file
   outfile = sonpath+'x_y_ss_raw'+str(p)+'.asc' 
   with open(outfile, 'w') as f:
      np.savetxt(f, np.hstack((humutils.ascol(X.flatten()),humutils.ascol(Y.flatten()), humutils.ascol(merge.flatten()))), delimiter=' ', fmt="%8.6f %8.6f %8.6f")

   humlon, humlat = trans(X, Y, inverse=True)

   if dogrid==1:
      grid_x, grid_y = np.meshgrid( np.arange(np.min(X), np.max(X), res), np.arange(np.min(Y), np.max(Y), res) )  

      dat = griddata(np.c_[X.flatten(),Y.flatten()], merge.flatten(), (grid_x, grid_y), method='nearest') 

      ## create mask for where the data is not
      tree = KDTree(np.c_[X.flatten(),Y.flatten()])
      dist, _ = tree.query(np.c_[grid_x.ravel(), grid_y.ravel()], k=1)
      dist = dist.reshape(grid_x.shape)

   del X, Y #, bearing #, pix_m, yvec

   if dogrid==1:
      ## mask
      dat[dist> np.floor(np.sqrt(1/res))-1 ] = np.nan #np.floor(np.sqrt(1/res))-1 ] = np.nan
      del dist, tree

      dat[dat==0] = np.nan
      dat[np.isinf(dat)] = np.nan
      datm = np.ma.masked_invalid(dat)

      glon, glat = trans(grid_x, grid_y, inverse=True)
      del grid_x, grid_y

   print "drawing and printing map ..."
   fig = plt.figure(frameon=False)
   map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #26949,
    resolution = 'i', #h #f
    llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
    urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

   if dogrid==1:
      gx,gy = map.projtran(glon, glat)

   ax = plt.Axes(fig, [0., 0., 1., 1.], )
   ax.set_axis_off()
   fig.add_axes(ax)

   if dogrid==1:
#.........这里部分代码省略.........

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
#Get grid of coordinates from raster
xx_1, yy_1 = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres]
trans =  pyproj.Proj(init="epsg:26949") 
glon_1, glat_1 = trans(xx_1, yy_1, inverse=True)

fig = plt.figure(figsize=(12,24))
plt.suptitle('May 2014')
ax = fig.add_subplot(1,3,1)
ax.set_title('R01359')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.min(glon_1 - 0.0007), 
            llcrnrlat=np.min(glat_1 - 0.0006),
            urcrnrlon=np.max(glon_1 + 0.0005), 
            urcrnrlat=np.max(glat_1 + 0.0006))
x,y = m.projtran(glon, glat)
m.wmsimage(server='http://grandcanyon.usgs.gov/arcgis/services/Imagery/ColoradoRiverImageryExplorer/MapServer/WmsServer?', 
              layers=['0'], xpixels=1000)
im = m.contourf(x, y, data_59.T, cmap='Greys_r',levels=ss_level)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.1)
cbr = plt.colorbar(im,cax=cax)

ax_2 = fig.add_subplot(1,3,2)
ax_2.set_title('R01360')
m = Basemap(projection='merc', 
            epsg=cs2cs_args.split(':')[1], 
            llcrnrlon=np.min(glon_1 - 0.0007), 
            llcrnrlat=np.min(glat_1 - 0.0006),
            urcrnrlon=np.max(glon_1 + 0.0005), 
            urcrnrlat=np.max(glat_1 + 0.0006))

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]

#.........这里部分代码省略.........

          dat = griddata(np.c_[X.flatten(),Y.flatten()], merge.flatten(), (grid_x, grid_y), method='nearest')

          ## create mask for where the data is not
          tree = KDTree(np.c_[X.flatten(),Y.flatten()])
          dist, _ = tree.query(np.c_[grid_x.ravel(), grid_y.ravel()], k=1)
          dist = dist.reshape(grid_x.shape)

       del X, Y

       if dogrid==1:
          ## mask
          dat[dist> 1 ] = np.nan 

          del dist, tree

          dat[dat==0] = np.nan
          dat[np.isinf(dat)] = np.nan

          datm = np.ma.masked_invalid(dat)

          glon, glat = trans(grid_x, grid_y, inverse=True)
          del grid_x, grid_y

       try:
          print "drawing and printing map ..."
          fig = plt.figure(frameon=False)
          map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #26949,
           resolution = 'i', #h #f
           llcrnrlon=np.min(humlon)-0.001, llcrnrlat=np.min(humlat)-0.001,
           urcrnrlon=np.max(humlon)+0.001, urcrnrlat=np.max(humlat)+0.001)

          if dogrid==1:
             gx,gy = map.projtran(glon, glat)

          ax = plt.Axes(fig, [0., 0., 1., 1.], )
          ax.set_axis_off()
          fig.add_axes(ax)

          if dogrid==1:
             map.pcolormesh(gx, gy, datm, cmap='YlOrRd', vmin=0.5, vmax=2)
             del dat
          else: 
             ## draw point cloud
             x,y = map.projtran(humlon, humlat)
             map.scatter(x.flatten(), y.flatten(), 0.5, merge.flatten(), cmap='YlOrRd', linewidth = '0')

          custom_save(sonpath,'class_map'+str(p))
          del fig 

       except:
          print "error: map could not be created..."

       kml = simplekml.Kml()
       ground = kml.newgroundoverlay(name='GroundOverlay')
       ground.icon.href = 'class_map'+str(p)+'.png'
       ground.latlonbox.north = np.min(humlat)-0.001
       ground.latlonbox.south = np.max(humlat)+0.001
       ground.latlonbox.east =  np.max(humlon)+0.001
       ground.latlonbox.west =  np.min(humlon)-0.001
       ground.latlonbox.rotation = 0

       kml.save(sonpath+'class_GroundOverlay'+str(p)+'.kml')

    if dowrite==1:

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]
ortnMapE=Basemap(projection='ortho',lat_0=latCtr,lon_0=lonCtr,resolution='c',
                 ax=ax[1])
ortnMapE.drawcoastlines(linewidth=0.5)
ortnMapE.drawcountries(linewidth=0.25)
ortnMapE.drawmeridians(numpy.arange(0,360,30))
ortnMapE.drawparallels(numpy.arange(-90,90,30))
ax[1].set_title(r'${}\ over\ Europe$'.format(shapeRdr0.records()[0][4]),
    fontsize=titleFontSize)
plotPrefecture(shp=shape,colour='gold',lwdth=2,bMap=ortnMapE,axes=ax[1],
               latOff=dLatJ-latJpn,longOff=dLonJ-lonJpn)

fig.show()

' Japan and Kitakyushu overlaid on Europe.'
fig,ax=matplotlib.pyplot.subplots(1,1,figsize=(16,8))
mercMapE=Basemap(projection='merc',llcrnrlat=30,urcrnrlat=75,llcrnrlon=-25,
                urcrnrlon=40,lat_ts=10,ax=ax,resolution='l')
mercMapE.drawcoastlines(linewidth=0.5)
mercMapE.drawcountries(linewidth=0.25)
mercMapE.drawparallels(numpy.arange(mercMapE.latmin,mercMapE.latmax,10.))
mercMapE.drawmeridians(numpy.arange(mercMapE.lonmin,mercMapE.lonmax,15.))
ax.set_title(r'$Europe,\ true\ lat.$',fontsize=titleFontSize)
plotPrefecture(shp=shape,colour='gold',lwdth=2,bMap=mercMapE,axes=ax,
               latOff=0,longOff=dLonJ-lonJpn)
# Show annotation at the true latitude.
xKIT,yKIT=mercMapE.projtran(130.834730+dLonJ-lonJpn,33.8924837)
xTXT,yTXT=mercMapE.projtran(110.834730+dLonJ-lonJpn,45.8924837)
ax.scatter([xKIT],[yKIT],s=50,c='crimson')
ax.annotate('Here', xy=(xKIT,yKIT),xytext=(xTXT,yTXT),color='crimson',
	arrowprops=dict(facecolor='crimson', shrink=0.05))
fig.show()

# 需要导入模块: from mpl_toolkits.basemap import Basemap [as 别名]
# 或者: from mpl_toolkits.basemap.Basemap import projtran [as 别名]

#.........这里部分代码省略.........

            targ_def = pyresample.geometry.SwathDefinition(lons=longrid.flatten(), lats=latgrid.flatten())
            del longrid, latgrid

            orig_def = pyresample.geometry.SwathDefinition(lons=humlon.flatten(), lats=humlat.flatten()) 

            dat, stdev, counts = pyresample.kd_tree.resample_gauss(orig_def, merge.flatten(), targ_def, radius_of_influence=influence, neighbours=nn, sigmas=sigmas, fill_value=None, with_uncert = np.nan, nprocs = cpu_count(), epsilon = eps)

      del X, Y

      dat = dat.reshape(shape)

      if mode>1:
         stdev = stdev.reshape(shape)
         counts = counts.reshape(shape)

      mask = dat.mask.copy()

      dat[mask==1] = 0

      if mode>1:
         dat[(stdev>numstdevs) & (mask!=0)] = np.nan
         dat[(counts<nn) & (counts>0)] = np.nan

      dat2 = replace_nans.RN(dat.astype('float64'),1000,0.01,2,'localmean').getdata()
      dat2[dat==0] = np.nan

      # get a new mask
      mask = np.isnan(dat2)

      mask = ~binary_dilation(binary_erosion(~mask,structure=np.ones((15,15))), structure=np.ones((15,15)))
      #mask = binary_fill_holes(mask, structure=np.ones((15,15)))
      #mask = ~binary_fill_holes(~mask, structure=np.ones((15,15)))

      dat2[mask==1] = np.nan
      dat2[dat2<1] = np.nan

      del dat
      dat = dat2
      del dat2


   if dogrid==1:
      ### mask
      #if np.floor(np.sqrt(1/res))-1 > 0.0:
      #   dat[dist> np.floor(np.sqrt(1/res))-1 ] = np.nan #np.floor(np.sqrt(1/res))-1 ] = np.nan
      #else:
      #   dat[dist> np.sqrt(1/res) ] = np.nan #np.floor(np.sqrt(1/res))-1 ] = np.nan

      #del dist, tree

      dat[dat==0] = np.nan
      dat[np.isinf(dat)] = np.nan
      datm = np.ma.masked_invalid(dat)

      glon, glat = trans(grid_x, grid_y, inverse=True)
      del grid_x, grid_y

   try:
      print "drawing and printing map ..."
      fig = plt.figure(frameon=False)
      map = Basemap(projection='merc', epsg=cs2cs_args.split(':')[1], #26949,
       resolution = 'i', #h #f
       llcrnrlon=np.min(humlon)-0.00001, llcrnrlat=np.min(humlat)-0.00001,
       urcrnrlon=np.max(humlon)+0.00001, urcrnrlat=np.max(humlat)+0.00001)

      if dogrid==1:
         gx,gy = map.projtran(glon, glat)

      ax = plt.Axes(fig, [0., 0., 1., 1.], )
      ax.set_axis_off()
      fig.add_axes(ax)

      if dogrid==1:
         map.pcolormesh(gx, gy, datm, cmap='gray', vmin=np.nanmin(datm), vmax=np.nanmax(datm))
         del datm, dat
      else: 
         ## draw point cloud
         x,y = map.projtran(humlon, humlat)
         map.scatter(x.flatten(), y.flatten(), 0.5, merge.flatten(), cmap='gray', linewidth = '0')

      custom_save(sonpath,'map'+str(p))
      del fig 

   except:
      print "error: map could not be created..."

   kml = simplekml.Kml()
   ground = kml.newgroundoverlay(name='GroundOverlay')
   ground.icon.href = 'map'+str(p)+'.png'
   ground.latlonbox.north = np.min(humlat)-0.00001
   ground.latlonbox.south = np.max(humlat)+0.00001
   ground.latlonbox.east =  np.max(humlon)+0.00001
   ground.latlonbox.west =  np.min(humlon)-0.00001
   ground.latlonbox.rotation = 0

   #kml.save(sonpath+'GroundOverlay'+str(p)+'.kml')
   kml.save(os.path.normpath(os.path.join(sonpath,'GroundOverlay'+str(p)+'.kml')))

   del humlat, humlon