def test_nearest_remap(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(), self.area_def, 50000, segments=1)
     remap = kd_tree.resample_nearest(self.area_def, res.ravel(), swath_def, 5000, segments=1)
     cross_sum = remap.sum()
     expected = 22275.0
     self.assertEqual(cross_sum, expected, msg="Grid remapping nearest failed")

def main():
    conf = read_conf("config.txt")
    start = time.clock()
    # 观测数据
    hdf_L1B = glob.glob(str(conf['mod02-file']))
    L1B_obj = SD.SD(hdf_L1B[0], SD.SDC.READ)
    # 观测点坐标
    hdf_Geo = glob.glob(str(conf['mod03-file']))
    GEO_obj = SD.SD(hdf_Geo[0], SD.SDC.READ)
    
    swath_def = get_swath_def(GEO_obj)
    data = get_reflectance_data(L1B_obj)
    
    area_def = get_proj_area(zone=50, lon_min=115, lon_max=123, lat_min=37, lat_max=42)
    
    result = kd_tree.resample_nearest(swath_def, data, area_def, radius_of_influence=5000)
   
    
#    plt.axis()
#    plt.imshow(result)
#    plt.savefig(r'D:/mei/result.png')
    
    result = np.uint8(result * 255)
    img = Image.fromarray(result[:,:,4:], 'RGB')
    
    cw = ContourWriterAGG(conf['shp-path'])
    cw.add_coastlines(img, area_def, resolution='i', width=0.5)
    
    img.save(conf['img-output'])
        
    end = time.clock()
    print(end - start)

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()    

 def test_masked_full_multi(self):
     data = numpy.ones((50, 10))
     data[:, 5:] = 2
     mask1 = numpy.ones((50, 10))
     mask1[:, :5] = 0
     mask2 = numpy.ones((50, 10))
     mask2[:, 5:] = 0
     mask3 = numpy.ones((50, 10))
     mask3[:25, :] = 0
     data_multi = numpy.column_stack(
         (data.ravel(), data.ravel(), data.ravel()))
     mask_multi = numpy.column_stack(
         (mask1.ravel(), mask2.ravel(), mask3.ravel()))
     masked_data = numpy.ma.array(data_multi, mask=mask_multi)
     lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def,
                                    masked_data, self.area_def, 50000,
                                    fill_value=None, segments=1)
     expected_fill_mask = numpy.fromfile(os.path.join(os.path.dirname(__file__),
                                                      'test_files',
                                                      'mask_test_full_fill_multi.dat'),
                                         sep=' ').reshape((800, 800, 3))
     fill_mask = res.mask
     cross_sum = res.sum()
     expected = 357140.0
     self.assertAlmostEqual(cross_sum, expected,
                            msg='Failed to resample masked data')
     self.assertTrue(numpy.array_equal(fill_mask, expected_fill_mask),
                     msg='Failed to create fill mask on masked data')

 def test_nearest_empty_multi_masked(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     data_multi = numpy.column_stack((data.ravel(), data.ravel(), data.ravel()))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data_multi, self.area_def, 50000, segments=1, fill_value=None)
     self.assertEqual(res.shape, (800, 800, 3), msg="Swath resampling nearest empty multi masked failed")

 def test_orthoplot(self):
     area_def = utils.parse_area_file(os.path.join(os.path.dirname(__file__),
                                                   'test_files', 'areas.cfg'), 'ortho')[0]
     swath_def = geometry.SwathDefinition(self.lons, self.lats)
     result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
                                       radius_of_influence=20000,
                                       fill_value=None)
     plt = plot._get_quicklook(area_def, result)

 def test_nearest_mp(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(), self.area_def, 50000, nprocs=2, segments=1)
     cross_sum = res.sum()
     expected = 15874591.0
     self.assertEqual(cross_sum, expected, msg="Swath resampling mp nearest failed")

 def test_nearest_empty_masked(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(), self.area_def, 50000, segments=1, fill_value=None)
     cross_sum = res.mask.sum()
     expected = res.size
     self.assertTrue(cross_sum == expected, msg="Swath resampling nearest empty masked failed")

def data_reproject(data, area_def = None):
    '''
    将文件数据投影到指定区域。默认为数据本身的经纬度范围。
    返回一个掩膜数组，无效值被遮盖。
    '''
    if area_def is None:
        area_def = get_area_def(data)
    result = kd_tree.resample_nearest(get_swath_def(data), data.aot_550, area_def, radius_of_influence=5000)
    return np.ma.masked_equal(result, 0)

 def test_plate_carreeplot(self):
     area_def = utils.parse_area_file(os.path.join(os.path.dirname(__file__),
                                                   'test_files', 'areas.cfg'), 'pc_world')[0]
     swath_def = geometry.SwathDefinition(self.lons, self.lats)
     result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
                                       radius_of_influence=20000,
                                       fill_value=None)
     plt = plot._get_quicklook(area_def, result, num_meridians=0,
                               num_parallels=0)

 def test_easeplot(self):
     from pyresample import plot, kd_tree, geometry
     from pyresample import parse_area_file
     area_def = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_sh')[0]
     swath_def = geometry.SwathDefinition(self.lons, self.lats)
     result = kd_tree.resample_nearest(swath_def, self.tb37v, area_def,
                                       radius_of_influence=20000,
                                       fill_value=None)
     plot._get_quicklook(area_def, result)

 def test_nearest_segments(self):
     data = np.fromfunction(lambda y, x: y * x, (50, 10))
     lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(),
                                    self.area_def, 50000, segments=2)
     cross_sum = res.sum()
     expected = 15874591.0
     self.assertEqual(cross_sum, expected)

 def test_nearest_1d(self):
     data = numpy.fromfunction(lambda x, y: x * y, (800, 800))
     lons = numpy.fromfunction(lambda x: 3 + x / 100.0, (500,))
     lats = numpy.fromfunction(lambda x: 75 - x / 10.0, (500,))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(self.area_def, data.ravel(), swath_def, 50000, segments=1)
     cross_sum = res.sum()
     expected = 35821299.0
     self.assertEqual(res.shape, (500,), msg="Swath resampling nearest 1d failed")
     self.assertEqual(cross_sum, expected, msg="Swath resampling nearest 1d failed")

 def test_nearest_multi_unraveled(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     data_multi = numpy.dstack((data, data, data))
     res = kd_tree.resample_nearest(swath_def, data_multi, self.area_def, 50000, segments=1)
     cross_sum = res.sum()
     expected = 3 * 15874591.0
     self.assertEqual(cross_sum, expected, msg="Swath multi channel resampling nearest failed")

 def test_masked_nearest_1d(self):
     data = numpy.ones((800, 800))
     data[:400, :] = 2
     lons = numpy.fromfunction(lambda x: 3 + x / 100.0, (500,))
     lats = numpy.fromfunction(lambda x: 75 - x / 10.0, (500,))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     mask = numpy.ones((800, 800))
     mask[400:, :] = 0
     masked_data = numpy.ma.array(data, mask=mask)
     res = kd_tree.resample_nearest(self.area_def, masked_data.ravel(), swath_def, 50000, segments=1)
     self.assertEqual(res.mask.sum(), 108, msg="Swath resampling masked nearest 1d failed")

 def test_nearest_empty(self):
     data = numpy.fromfunction(lambda y, x: y * x, (50, 10))
     lons = numpy.fromfunction(lambda y, x: 165 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(),
                                    self.area_def, 50000, segments=1)
     cross_sum = res.sum()
     expected = 0
     self.assertEqual(cross_sum, expected,
                      msg='Swath resampling nearest empty failed')

 def test_masked_fill_int(self):
     data = numpy.ones((50, 10)).astype("int")
     lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(), self.area_def, 50000, fill_value=None, segments=1)
     expected_fill_mask = numpy.fromfile(
         os.path.join(os.path.dirname(__file__), "test_files", "mask_test_fill_value.dat"), sep=" "
     ).reshape((800, 800))
     fill_mask = res.mask
     self.assertTrue(
         numpy.array_equal(fill_mask, expected_fill_mask), msg="Failed to create fill mask on integer data"
     )

 def test_masked_fill_int(self):
     data = np.ones((50, 10)).astype('int')
     lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
     swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
     res = kd_tree.resample_nearest(swath_def, data.ravel(),
                                    self.area_def, 50000, fill_value=None, segments=1)
     expected_fill_mask = np.fromfile(os.path.join(os.path.dirname(__file__),
                                                   'test_files',
                                                   'mask_test_fill_value.dat'),
                                      sep=' ').reshape((800, 800))
     fill_mask = res.mask
     self.assertTrue(np.array_equal(fill_mask, expected_fill_mask))

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')    

 def test_nearest_masked_swath_target(self):
     """Test that a masked array works as a target."""
     data = np.fromfunction(lambda y, x: y * x, (50, 10))
     lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
     lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
     mask = np.ones_like(lons, dtype=np.bool)
     mask[::2, ::2] = False
     swath_def = geometry.SwathDefinition(
         lons=np.ma.masked_array(lons, mask=mask),
         lats=np.ma.masked_array(lats, mask=False)
     )
     res = kd_tree.resample_nearest(swath_def, data.ravel(),
                                    swath_def, 50000, segments=3)
     cross_sum = res.sum()
     # expected = 12716  # if masks aren't respected
     expected = 12000
     self.assertEqual(cross_sum, expected)