0

I want to use XML auxiliary files in Sentinel-1 files to remove noise of HV-pol images, I have extracted these noise GCP in noise XML files into a point shapefile using python( in below ). However, I don’t know how to rebuild the noise raster files?

I just know linear interpolation means, and I have tried as enter link description here.

However, the S-1 image size is 8992 * 10376 , script running error shows, MemoryError . Why popup the error message?

I posted my script includes to set noise field to those points and from point to raster in the end.

enter image description here

def Create_Points( root_dir, noise_dir, noise_value, ROWS, COLS, num_rows ,num_cols, rows, c ):

    env.overwriteOutput = True
    outpath = root_dir + noise_dir + r'measurement'
    newfc = 'newpoint.shp'
    ## Create points shapefiles
    arcpy.CreateFeatureclass_management( outpath, newfc, "Point" )
    newfc = os.path.join( outpath , newfc )
    cursor = arcpy.da.InsertCursor( newfc , [ "SHAPE@" ] )
    for i in range( len( rows ) ):
        for j in range( len( c ) ):
            cursor.insertRow( [ arcpy.Point( float( rows[ i ] ), float( c[ j ] ) ) ] )
    del cursor

    ## Add a new field to store noise values.
    source = ogr.Open( newfc, update=True )
    layer = source.GetLayer()
    noise_field = ogr.FieldDefn( 'noise', ogr.OFTReal )
    layer.CreateField( noise_field )
    #To set value to noise field.
    fields = [ 'FID', 'noise' ]
    with arcpy.da.UpdateCursor( newfc, fields ) as cursor:    
        for row in cursor:
            row [ 1 ] = float( noise_value[ row[ 0 ] ] )
            cursor.updateRow( row )
    # Close the Shapefile
    source = None

    ##Create noise tiff file based on point shapefile.
    # Reshape noise_value to a matrix of [ num_rows, num_cols ] size.
    noise_value = np.asarray( noise_value ).reshape( num_rows ,num_cols )
    position = np.chararray( ( num_rows , num_cols ), 10 )
    for i in range( num_rows ):
        for j in range( num_cols ):        
            position[ i ] [ j ] = str( rows[ i ] ) +',' + str( c[ j ] )    
    position = np.asarray( position ).ravel()
    location = np.zeros( ( num_rows * num_cols, 2 ) )
    for id in range( len( position ) ):
        location[ id ] [ 0 ] = int( ( position[ id ] ).split(",")[ 0 ] )
        location[ id ] [ 1 ] = int( ( position[ id ] ).split(",")[ 1 ] )
    #Convert noise_value to one dim array.
    rows_min, rows_max, c_min, c_max = [ min( rows ), max( rows ), min( c ), max( c ) ]
    n_rows = ( int( rows_max ) - int( rows_min ) + 1 )
    n_cols = ( int( c_max ) - int( c_min ) + 1 )
    xi = np.linspace( int( rows_min ), int( rows_max ), n_rows/20 )
    yi = np.linspace( int( c_min ), int( c_max ), n_cols/20 )

    grid_x, grid_y = np.meshgrid( xi, yi )
    noise_value = noise_value.reshape( num_rows * num_cols )
    noise_inter = griddata( location, noise_value, ( grid_x, grid_y ), method = 'linear' )
    noise_inter = np.kron( noise_inter, np.ones(( 20, 20 )))
    ##Create the output image.
    driver = gdal.GetDriverByName( 'GTiff' )
    outputRaster = os.path.join( outpath, 'noise.tif' )
    dst_ds = driver.Create( outputRaster, COLS, ROWS, 1, gdal.GDT_Int16 )
    band = dst_ds.GetRasterBand( 1 )
    band.WriteArray( noise_inter )
Falcon_Baikal
  • 121
  • 1
  • 11

1 Answers1

0

I have solved it.

For first question, I using another method, scipy.interpolate.griddata function(https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html) in python.

For second, Memory Error is because my spyder before is 32-bit and S-1 images size is too large. Then I installed 64-bit Spyder and ArcGIS for Desktop Background Geoprocessing(64-bit). In 64-bit Spyder environment, two trouble and solved strateyies is at the hyperlink, http://yangbaikal.blogspot.com/2018/10/import-gdal-and-arcpy-error-in-spyder.html .

Falcon_Baikal
  • 121
  • 1
  • 11