Source code for geophpy.operation.general

# -*- coding: utf-8 -*-
'''
   geophpy.operation.general
   -------------------------

   DataSet Object general operations routines.

   :copyright: Copyright 2014 Lionel Darras, Philippe Marty and contributors, see AUTHORS.
   :license: GNU GPL v3.

'''
### USER DEFINED PARAMETERS ##########################################

# list of "griddata" interpolation methods available for wumappy interface
gridding_interpolation_list = ['none', 'nearest', 'linear', 'cubic']

######################################################################



import numpy as np
from scipy.interpolate import griddata, RectBivariateSpline
from copy import deepcopy
from geophpy.misc.utils import *



def copy(dataset):
   '''
   cf. dataset.py

   '''
###ORIGINAL CODE having some flaws:
###- newdataset was a parameter instead of a "return"
###- values were correctly handled but the code was sequential
###- zimage was duplicated as a tuple instead of an np.array
###- fields and info were not duplicated but "referenced"
###- georef was missing ...
   #newdataset.data.values = []
   #for l in range(len(dataset.data.values)):
      ## for each line adds a value in the column
      #newdataset.data.values.append(dataset.data.values[l])
   #newdataset.data.values = np.array(newdataset.data.values)
   #newdataset.data.z_image = []
   #for l in range(len(dataset.data.z_image)):
      #newdataset.data.z_image.append(dataset.data.z_image[l])
   #newdataset.data.fields = dataset.data.fields
   #newdataset.info = dataset.info
###SECOND CODE somewhat better:
   ## initialisation of the duplicate DataSet Object
   #newdataset = DataSet._new()
   ## duplication of the raw values
   #newdataset.data.values = dataset.data.values.copy()
   ## duplication of the zimage
   #newdataset.data.z_image = dataset.data.z_image.copy()
   ## duplication of the fieldnames
   #newdataset.data.fields = dataset.data.fields.copy()
   ## duplication of the dataset meta data
   #newdataset.info = deepcopy(dataset.info)
   ## duplication of the georeferencing meta data
   #newdataset.georef.active = (dataset.georef.active & True)
   #newdataset.georef.points_list = dataset.georef.points_list.copy()
   #return newdataset
###THIRD CODE straightforward:
   return deepcopy(dataset)



def getgriddinginterpolationlist():
   """
   cf. dataset.py

   """
   return gridding_interpolation_list



def interpolate(dataset, interpolation="none", x_delta=None, y_delta=None, x_decimal_maxnb=2, y_decimal_maxnb=2, x_frame_factor=0., y_frame_factor=0.):
   '''
   cf. dataset.py

   '''
   T = dataset.data.values.T
   x_array = T[0]
   y_array = T[1]
   z_array = T[2]

   # Get distinct x values ###########################################
   x_list = np.unique(x_array)

   # Get the median delta between two distinct x values ##############
   if (x_delta == None):
      if (x_decimal_maxnb == None):
         # Undefined x rounding precision ############################
         # ...TBD... raise an error here !
         pass
      x_delta = np.median(np.around(np.diff(x_list), x_decimal_maxnb))
      # ...TBD... why not take the min diff value instead of the median ?
      #x_delta = np.nanmin(np.around(np.diff(x_list), x_decimal_maxnb))
   else:
      x_decimal_maxnb = getdecimalsnb(x_delta)

   # Determinate min and max x coordinates and number of x pixels ####
   xmin = x_array.min()
   xmax = x_array.max()
   xmin = (1.+x_frame_factor)*xmin - x_frame_factor*xmax
   xmax = (1.+x_frame_factor)*xmax - x_frame_factor*xmin
   xmin = round(xmin, x_decimal_maxnb)
   xmax = round(xmax, x_decimal_maxnb)
   # Don't forget the "+1" (éternel pb poteaux et intervalles !) #####
   nx   = np.around((xmax-xmin)/x_delta) + 1

   # Get the median delta between two distinct y values ##############
   if (y_delta == None):
      if (y_decimal_maxnb == None):
         # Undefined y rounding precision ############################
         # ...TBD... raise an error here !
         pass
      # Get median value for each profile ############################
      deltamedianlist = []
      for x in x_list:
         # Extract the y profile at x coordinate #####################
         profile = np.unique(y_array[np.where(x_array == x)])
         # Append the median of this y profile #######################
         deltamedianlist.append(np.median(np.around(np.diff(profile), y_decimal_maxnb)))
         # ...TBD... why not take the min diff value instead of the median ?
         #deltamedianlist.append(np.nanmin(np.around(np.diff(profile), y_decimal_maxnb)))
      # Get median value of all median values ########################
      y_delta = np.median(np.array(deltamedianlist))
      # ...TBD... why not take the min diff value instead of the median ?
      #y_delta = np.nanmin(np.array(deltamedianlist))
   else:
      y_decimal_maxnb = getdecimalsnb(y_delta)

   # Determinate min and max y coordinates and number of y pixels ####
   ymin = y_array.min()
   ymax = y_array.max()
   ymin = (1.+y_frame_factor)*ymin - y_frame_factor*ymax
   ymax = (1.+y_frame_factor)*ymax - y_frame_factor*ymin
   ymin = round(ymin, y_decimal_maxnb)
   ymax = round(ymax, y_decimal_maxnb)
   # Don't forget the "+1" (éternel pb poteaux et intervalles !) #####
   ny   = np.around((ymax - ymin)/y_delta) + 1

   # Build the (xi, yi) regular vector coordinates ###################
   xi = np.linspace(xmin, xmax, nx, endpoint=True)
   yi = np.linspace(ymin, ymax, ny, endpoint=True)

   # Build the (X, Y) regular grid coordinates #######################
   X, Y = np.meshgrid(xi, yi)

   # Build the Zimage regularly gridded matrix #######################
   if (interpolation == "none"):
      ## just project data into the grid
      ## if several data points fall into the same pixel, they are averaged
      ## don't forget to "peakfilt" the rawvalues beforehand to avoid averaging bad data points
      Z = X * 0.
      P = Z.copy()
      for x,y,z in dataset.data.values:
         indx = np.where(xi+x_delta/2. > x)
         indy = np.where(yi+y_delta/2. > y)
         Z[indy[0][0],indx[0][0]] += z
         P[indy[0][0],indx[0][0]] += 1
      Z = Z/P
   elif (interpolation in getgriddinginterpolationlist()) :
      ## perform data interpolation onto the grid
      ## the interpolation algorithm will deal with overlapping data points
      ## nevertheless don't forget to "peakfilt" the rawvalues beforehand to avoid interpolation being too much influenced by bad data points
      '''
      # Fill holes in each profiles with "nan" #######################
      ## this is to avoid filling holes with interpolated values
      nan_array = []
      for x in x_list:
         profile = np.unique(y_array[np.where(x_array == x)])
         nan_array = profile_completewithnan(x, profile, nan_array, y_delta, factor=2, ymin=ymin, ymax=ymax)
      if (len(nan_array) != 0):
         completed_array = np.append(dataset.data.values, np.array(nan_array), axis=0)
         T = completed_array.T
         x_array = T[0]
         y_array = T[1]
         z_array = T[2]
      '''
      # Interpolate the Zimage #######################################
      Z = griddata((x_array, y_array), z_array, (X, Y), method=interpolation)
   else:
      # Undefined interpolation method ###############################
      # ...TBD... raise an error here !
      pass

   # Fill the DataSet Object #########################################
   dataset.data.z_image = Z
   dataset.info.x_min = xmin
   dataset.info.x_max = xmax
   dataset.info.y_min = ymin
   dataset.info.y_max = ymax
   dataset.info.z_min = np.nanmin(Z)
   dataset.info.z_max = np.nanmax(Z)
   dataset.info.x_gridding_delta = x_delta
   dataset.info.y_gridding_delta = y_delta
   dataset.info.gridding_interpolation = interpolation



#def apodisation2d(val, apodisation_factor):
   '''
   2D apodisation, to reduce side effects

   Parameters :

   :val: 2-Dimension array

   :apodisation_factor: apodisation factor in percent (0-25)

   '''
#   if (apodisation_factor > 0):
#      # apodisation in the x direction
#      for profile in val.T:
#         _apodisation1d(profile, apodisation_factor)

      # apodisation in the y direction
#      for profile in val:
#         _apodisation1d(profile, apodisation_factor)



#def _apodisation1d(array1D, apodisation_factor):
   '''
   1D apodisation, to reduce side effects

   Parameters :

   :array1D: 1-Dimension array

   :apodisation_factor: apodisation factor in percent (0-25)

   '''
#   na = len(array1D)                                  # n is the number of array elements
#   napod = int(np.around((na * apodisation_factor)/100))     # napod is the number of array elements to treat
#   if (napod <= 1):                                   # one element at least must be treated
#      napod = 1
#   pi2 = np.pi/2.
#   for n in range(napod):                             # for napod first data
#      array1D[n] = array1D[n]*np.cos((napod-n)*pi2/napod)

#   for n in range(na-napod, na):                      # for napod last data
#      array1D[n] = array1D[n]*np.cos((n+1-na+napod)*pi2/napod)



[docs]def apodisation2d(val, apodisation_factor): ''' 2D apodisation, to reduce side effects Parameters : :val: 2-Dimension array :apodisation_factor: apodisation factor in percent (0-25) Returns : - apodisation pixels number in x direction - apodisation pixels number in y direction - enlarged array after apodisation ''' array2DTemp = [] array2D = [] if (apodisation_factor > 0): # apodisation in the x direction nx = len(val.T[0]) # n is the number of array elements napodx = int(np.around((nx * apodisation_factor)/100)) # napod is the number of array elements to treat if (napodx <= 1): # one element at least must be treated napodx = 1 for profile in val.T: array2DTemp.append(_apodisation1d(profile, napodx)) array2DTemp = (np.array(array2DTemp)).T # apodisation in the y direction ny = len(array2DTemp[0]) # n is the number of array elements napody = int(np.around((ny * apodisation_factor)/100)) # napod is the number of array elements to treat if (napody <= 1): # one element at least must be treated napody = 1 for profile in array2DTemp: array2D.append(_apodisation1d(profile, napody)) else: # apodisation factor = 0 array2D = val # return napodx, napody, np.array(array2D) return np.array(array2D)
def _apodisation1d(array1D, napod): ''' 1D apodisation, to reduce side effects Parameters : :array1D: 1-Dimension array :napod: apodisation pixels number Returns : 1-Dimension array of len(array1D) + napod elements ''' pi2 = np.pi/2. na = len(array1D) # n is the number of array elements nresult = na + 2*napod array1Dresult = [] for n in range(napod): array1Dresult.append(array1D[n]*np.cos((napod-n)*pi2/napod)) for n in range(na): array1Dresult.append(array1D[n]) for n in range(na-napod, na): # for napod last data array1Dresult.append(array1D[n]*np.cos((n+1-na+napod)*pi2/napod)) return array1Dresult def apodisation2Dreverse(val, valwithapod, napodx, napody): ''' To do the reverse apodisation ''' na = len(val) nb = len(val[0]) for n in range(na): for m in range(nb): val[n][m] = valwithapod[n+napody][m+napodx] def sample(dataset): ''' To rebuild the Values from a Z_image ''' X = zimage_xcoord(dataset) Y = zimage_ycoord(dataset) Z = dataset.data.z_image xi = dataset.data.values[:,0] yi = dataset.data.values[:,1] zi = dataset.data.values[:,2] zCubSpl = RectBivariateSpline(Y, X, Z, kx=3, ky=3) # why Y,X ? zj = zCubSpl.ev(xi, yi, dx=0, dy=0) zi *= 0. zi += zj