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