pixelFunctionTest.py

logscale.vrt

<VRTDataset rasterXSize="9002" rasterYSize="5857">
  <SRS dataAxisToSRSAxisMapping="1,2">PROJCS["USA_Contiguous_Lambert_Conformal_Conic",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",39],PARAMETER["central_meridian",-96],PARAMETER["standard_parallel_1",33],PARAMETER["standard_parallel_2",45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]</SRS>
  <GeoTransform> -1.4748244013199583e+05,  2.5000000000000000e+02,  0.0000000000000000e+00,  7.4684201824673638e+05,  0.0000000000000000e+00, -2.5000000000000000e+02</GeoTransform>
  <VRTRasterBand dataType="Float32" band="1" subClass="VRTDerivedRasterBand">
    <NoDataValue>-32768</NoDataValue>
    <ColorInterp>Gray</ColorInterp>
    <ComplexSource>
      <SourceFilename relativeToVRT="1">../sources/CT_f32_gmted2010_srtm_merge_102004_1000m_avg_250m_cubicspline.tif</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="9002" RasterYSize="5857" DataType="Float32" BlockXSize="256" BlockYSize="256" />
      <SrcRect xOff="0" yOff="0" xSize="9001.78684747202" ySize="5856.50952181305" />
      <DstRect xOff="0.213152527980506" yOff="0.490478186946362" xSize="9001.78684747202" ySize="5856.50952181305" />
      <NODATA>-32768</NODATA>
    </ComplexSource>
    <PixelFunctionLanguage>Python</PixelFunctionLanguage>
<PixelFunctionType>pixelFunctionTest.runOperation</PixelFunctionType>
<PixelFunctionArguments op="globalLogScaleLandA" />
  </VRTRasterBand>
</VRTDataset>

pixelFunctionTest.py

import numpy as np

import math

def globalLogScaleLandA(a):
  if a > 0:
    return a+0.1 #math.log(float(a))
  elif a == 0:
    return 0
  elif a > -32768:
    return -1 * math.log(float(abs(a)))
  else:
    return a
  
def runOperation(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, radius, gt, **kwargs):
  op = kwargs['op'].decode('utf-8')
  vGlobalLogScaleLandA = np.vectorize(globalLogScaleLandA)
  np.round_(vGlobalLogScaleLandA(in_ar[0]), out=out_ar)