Fix colour composite

ERS_NRB/ancillary.py

@@ -141,80 +141,95 @@ def vrt_relpath(vrt):
     etree.indent(tree.getroot())
     tree.write(vrt, pretty_print=True, xml_declaration=False, encoding='utf-8')
 
-
 def create_rgb_vrt(outname, infiles, overviews, overview_resampling):
     """
-    Creates the RGB VRT file.
-    
+    Creation of the color composite VRT file.
     Parameters
     ----------
     outname: str
-        Full path to the output RGB VRT file.
+        Full path to the output VRT file.
     infiles: list[str]
-        A list of paths.
+        A list of paths pointing to the linear scaled measurement backscatter files.
     overviews: list[int]
         Internal overview levels to be defined for the created VRT file.
     overview_resampling: str
-        Resampling method for overview levels.
+        Resampling method applied to overview pyramids.
+    """    
+    # make sure order is right and co-polarization (VV or HH) is first
+    pols = [re.search('[hv]{2}', os.path.basename(f)).group() for f in infiles]
+    if pols[1] in ['vv', 'hh']:
+        infiles.reverse()
+        pols.reverse()
+
+    # format overview levels
+    ov = str(overviews)
+    for x in ['[', ']', ',']:
+        ov = ov.replace(x, '')
 
-    Returns
-    -------
-    None
-    """
-    # make sure order is right and VV polarization is first
-    paths_reorder = []
-    for i, f in enumerate(infiles):
-        pol = re.search('[hv]{2}', os.path.basename(f)).group()
-        if i == 1 and pol == 'vv':
-            paths_reorder.append(f)
-            paths_reorder.append(infiles[0])
-            infiles = paths_reorder
-
-    # create VRT and change its content
-    gdalbuildvrt(src=infiles, dst=outname, options={'separate': True})
+    # create VRT file and change its content
+    gdalbuildvrt(src=infiles, dst=outname,
+                 options={'separate': True})
 
     tree = etree.parse(outname)
     root = tree.getroot()
+    srs = tree.find('SRS').text
+    geotrans = tree.find('GeoTransform').text
     bands = tree.findall('VRTRasterBand')
+    vrt_nodata = bands[0].find('NoDataValue').text
+    complex_src = [band.find('ComplexSource') for band in bands]
+    for cs in complex_src:
+        cs.remove(cs.find('NODATA'))
 
     new_band = etree.SubElement(root, 'VRTRasterBand',
-                                attrib={'dataType': 'Float32', 'band': str(len(bands)), 'subClass': 'VRTDerivedRasterBand'})
-    vrt_nodata = etree.SubElement(new_band, 'NoDataValue')
-    vrt_nodata.text = 'nan'
-    for i, band in enumerate(bands, start=1):
-        new_band.insert(i, deepcopy(band.find('ComplexSource')))
-    # pxfun_type = etree.SubElement(new_band, 'PixelFunctionType')
-    # pxfun_type.text = 'diff'
-    pxfun_language = etree.SubElement(new_band, 'PixelFunctionLanguage')
-    pxfun_language.text = 'Python'
+                                attrib={'dataType': 'Float32', 'band': '3',
+                                        'subClass': 'VRTDerivedRasterBand'})
+    new_band_na = etree.SubElement(new_band, 'NoDataValue')
+    new_band_na.text = 'nan'
     pxfun_type = etree.SubElement(new_band, 'PixelFunctionType')
-    pxfun_type.text = 'div'
-    pxfun_code = etree.SubElement(new_band, 'PixelFunctionCode')
-    pxfun_code.text = etree.CDATA("""
-import numpy as np
-def div(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
-    np.divide(in_ar[0], in_ar[1], out=out_ar, where=in_ar[1]!=0, dtype='float32')
-    """)
+    pxfun_type.text = 'mul'
+    for cs in complex_src:
+        new_band.append(deepcopy(cs))
+
+    src = new_band.findall('ComplexSource')[1]
+    fname = src.find('SourceFilename')
+    fname_old = fname.text
+    src_attr = src.find('SourceProperties').attrib
+    fname.text = etree.CDATA("""
+    <VRTDataset rasterXSize="{rasterxsize}" rasterYSize="{rasterysize}">
+        <SRS dataAxisToSRSAxisMapping="1,2">{srs}</SRS>
+        <GeoTransform>{geotrans}</GeoTransform>
+        <VRTRasterBand dataType="{dtype}" band="1" subClass="VRTDerivedRasterBand">
+            <NoDataValue>{vrt_nodata}</NoDataValue>
+            <PixelFunctionType>{px_fun}</PixelFunctionType>
+            <ComplexSource>
+              <SourceFilename relativeToVRT="1">{fname}</SourceFilename>
+              <SourceBand>1</SourceBand>
+              <SourceProperties RasterXSize="{rasterxsize}" RasterYSize="{rasterysize}" DataType="{dtype}" BlockXSize="{blockxsize}" BlockYSize="{blockysize}"/>
+              <SrcRect xOff="0" yOff="0" xSize="{rasterxsize}" ySize="{rasterysize}"/>
+              <DstRect xOff="0" yOff="0" xSize="{rasterxsize}" ySize="{rasterysize}"/>
+            </ComplexSource>
+        </VRTRasterBand>
+        <OverviewList resampling="{ov_resampling}">{ov}</OverviewList>
+    </VRTDataset>
+    """.format(rasterxsize=src_attr['RasterXSize'], rasterysize=src_attr['RasterYSize'], srs=srs, geotrans=geotrans,
+               dtype=src_attr['DataType'], px_fun='inv', fname=fname_old, vrt_nodata=vrt_nodata,
+               blockxsize=src_attr['BlockXSize'], blockysize=src_attr['BlockYSize'],
+               ov_resampling=overview_resampling.lower(), ov=ov))
 
     bands = tree.findall('VRTRasterBand')
     for band, col in zip(bands, ['Red', 'Green', 'Blue']):
         color = etree.Element('ColorInterp')
         color.text = col
-        band.insert(1, color)
-        if band.attrib['band'] in ['1', '2']:
-            ndv = band.find('NoDataValue')
-            band.remove(ndv)
+        band.insert(0, color)
 
     ovr = etree.SubElement(root, 'OverviewList', attrib={'resampling': overview_resampling.lower()})
-    ov = str(overviews)
-    for x in ['[', ']', ',']:
-        ov = ov.replace(x, '')
     ovr.text = ov
 
     etree.indent(root)
     tree.write(outname, pretty_print=True, xml_declaration=False, encoding='utf-8')
 
 
+
 def generate_unique_id(encoded_str):
     """
     Returns a unique product identifier as a hexa-decimal string generated from the time of execution in isoformat.

docs/about/changelog.rst

@@ -1,12 +1,13 @@
 Changelog
 =========
 
-0.1.4 | 2022-10-07
+0.1.4 | 2022-10-10
 ------------------
 * Add compression as a config parameter
 * Compatibe with WSM
 * Use source range looks and azimuth looks
 * Write polarization code in final product
+* Fix colour composite
 `Full Changelog <https://github.com/SAR-ARD/ERS_NRB/compare/0.1.3...0.1.4>`_
 
 0.1.3 | 2022-09-22