Add i.saocom and i.sar toolsets

src/imagery/i.saocom/Makefile

@@ -0,0 +1,13 @@
+MODULE_TOPDIR =../..
+
+PGM = i.saocom
+
+SUBDIRS = i.saocom.geocode \
+        i.saocom.import \
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: parsubdirs htmldir
+
+install: installsubdirs
+	$(INSTALL_DATA) $(PGM).html $(INST_DIR)/docs/html/

src/imagery/i.saocom/i.saocom.geocode/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.saocom.geocode
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

src/imagery/i.saocom/i.saocom.geocode/i.saocom.geocode.html

@@ -0,0 +1,52 @@
+<h2>DESCRIPTION</h2>
+
+The <em>i.saocom.geocode</em> module allows the projection of SAOCOM-1 L1A raw or derived products from radar coordinates to a cartographic coordinate system (CS). This can be applied to bands imported by <em>i.saocom.import</em>, any raster map calculated from these inputs, or directly to the SAOCOM original files. The tool does not currently support the geocoding of image subsets, so the whole extent must be processed.
+
+<p>
+<center>
+<img src="saocomGeocodeWF.png" width="1100" height="600">
+<p><em>Graphical description of the workflow followed by <em>i.saocom.geocode</em></em>
+</center>
+
+<h2>NOTES</h2>
+
+When the tool is directly run on the original SAOCOM files, it will internally call <em>i.saocom.import</em> to import the real and imaginary bands to a temporary XY unprojected location. In this case the <b>data</b> option must be used. If the user wants to geocode already imported files, the <b>map</b> option must be specified. The tool will not allow the use of both options, either one of them must be selected.
+
+This module makes use of the <b>GCPS.csv</b> file generated by <em>i.saocom.import</em> and associated to a specific image. This GCP file is associated to the image basename specified in the option <b>basename</b>.
+
+The tool will write the output to a Geotiff external file in the directory specified by the <b>dbase</b> option. The output files will be projected to the CS of the target location indicated in the option <b>location</b>. If the specified location does not exist, GRASS will create it with the CS EPSG:4326. The output map name will be the same as the input map, with the suffix <em>_geo</em> added at the end.
+
+
+<h2>EXAMPLES</h2>
+
+<h3>Within an XY unprojected location where the SAOCOM-1 bands have already been imported and processed, create an amplitude image from the real and imaginary bands, geocode it and then import it into a location called <em>saocom_geo</em></h3>
+
+<div class="code"><pre>
+#Calculate the amplitude image 
+i.sar.amplitude real=SAO1B_20211222_hh_real imag=SAO1B_20211222_hh_imag output=SAO1B_20211222_hh_amp
+
+#Geocode it and import it to new location. The external temporary raster will be saved as GeoTiff in the $HOME directory
+i.saocom.geocode map=SAO1B_20211222_hh_amp basename=SAO1B_20211222 dbase=$HOME location=saocom_geo
+</pre></div>
+
+<h3>Geocode the real and imaginary bands of a SAOCOM L1A image and import them to a location called <em>geo_location</em></h3>
+
+<div class="code"><pre>
+zip_file = 'S1B_OPER_SAR_EOSSP__CORE_L1A_OLF_20211225T165228.zip'
+i.saocom.geocode data=zip_file pols=['hh'] multilook=[8,4] basename='SAO1B_20211222' location='geo_location' dbase ='$HOME' )
+</pre></div>
+
+The previous example will run <em>i.saocom.import</em> to generate the real and imaginary bands for the polarimetric channels specified in the <b>pols</b> option and will apply
+the multilook factor indicated in the <b>multilook</b> option. After that it will geocode the real and imaginary bands and import them to the target location. 
+
+
+<h2>SEE ALSO</h2>
+
+<a href="i.saocom.import.html">i.saocom.import</a>,
+<a href="i.sar.amplitude.html">i.sar.amplitude</a>,
+<a href="i.sar.speckle.html">i.sar.speckle</a>,
+<a href="i.sar.pauliRGB.html">i.sar.pauliRGB</a>
+
+<h2>AUTHORS</h2>
+
+Santiago Seppi, CONAE, Argentina.

src/imagery/i.saocom/i.saocom.geocode/i.saocom.geocode.py

@@ -0,0 +1,294 @@
+#!/usr/bin/env python3
+
+############################################################################
+#
+# MODULE:    i.saocom.geocode
+#
+# AUTHOR:   Santiago Seppi
+#
+# PURPOSE:    Geocode SAOCOM SLC derived products using information from Ground Control Points (GCPs). The program writes the geocoded grid to a temporary external file, and then imports it to a new location and mapset.
+#
+# COPYRIGHT: (C) 2002-2023 by the GRASS Development Team
+#
+#         This program is free software under the GNU General Public
+#         License (>=v2). Read the file COPYING that comes with GRASS
+#         for details.
+#
+#
+#############################################################################
+
+# %Module
+# % description: Geocode SAOCOM SLC derived products
+# % keyword: imagery
+# % keyword: saocom
+# % keyword: sar
+# % keyword: radar
+# % overwrite: yes
+# %End
+# %option G_OPT_R_INPUT
+# % key: map
+# % required: no
+# % description: Map to be geocoded, in radar coordinates
+# %end
+# %option
+# % key: data
+# % type: string
+# % required: no
+# % multiple: no
+# % description: Path to data directory (ZIP or folder)
+# %end
+# %option
+# % key: is_zip
+# % type: string
+# % required: no
+# % multiple: no
+# % answer: yes
+# % options: yes,no
+# % description: Whether the data directory is zipped or not. Only required if geocoding is to ble applied on the original files
+# %end
+# %option
+# % key: polarizations
+# % type: string
+# % required: no
+# % multiple: yes
+# % answer: ['hh','hv','vh','vv']
+# % description: Polarizations to process . Only required if geocoding is to ble applied on the original files. (Default:['hh','hv','vh','vv'])
+# %end
+# %option
+# % key: multilook
+# % type: integer
+# % required: no
+# % multiple: yes
+# % answer: 1,1
+# % description: Azimuth and range factors to apply (eg: [4,2])). Only required if geocoding is to be applied on the original files
+# %end
+# %option
+# % key: basename
+# % type: string
+# % required: yes
+# % multiple: no
+# % description: Basename folder contaning the GCPS info. This folder is generated by i.saocom.import
+# %end
+# %option
+# % key: location
+# % type: string
+# % multiple: no
+# % required: yes
+# % description: Location where the geocoded map will be stored
+# %end
+# %option G_OPT_M_DIR
+# % key: dbase
+# % multiple: no
+# % required: yes
+# % description: GRASS GIS database directory to store the external geocoded map
+# %end
+
+
+import os
+import numpy as np
+import grass.script as gs
+import rasterio
+import numpy as np
+from osgeo import gdal
+import pandas as pd
+import shutil
+
+
+def save_GeoTiff(fn, crs, transform, mat, meta=None, nodata=None, bandnames=[]):
+    if len(mat.shape) == 2:
+        count = 1
+    else:
+        count = mat.shape[0]
+
+    if not meta:
+        meta = {}
+
+    meta["driver"] = "GTiff"
+    meta["height"] = mat.shape[-2]
+    meta["width"] = mat.shape[-1]
+    meta["count"] = count
+    meta["crs"] = crs
+    meta["transform"] = transform
+
+    if "dtype" not in meta:  # if no datatype is specified, use float32
+        meta["dtype"] = np.float32
+
+    if nodata is None:
+        pass
+    else:
+        meta["nodata"] = nodata
+
+    with rasterio.open(fn, "w", **meta) as dst:
+        if count == 1:  # Mono-band raster
+            dst.write(mat.astype(meta["dtype"]), 1)
+            if bandnames:
+                dst.set_band_description(1, bandnames[0])
+        else:  # Multi-band raster
+            for b in range(count):
+                dst.write(mat[b].astype(meta["dtype"]), b + 1)
+            for b, bandname in enumerate(bandnames):
+                dst.set_band_description(b + 1, bandname)
+
+
+def read_gcps(gcp_base_folder):
+    df_gcps = pd.read_csv(os.path.join(gcp_base_folder, "GCPS.csv"))
+    gcps = []
+    for i, d in df_gcps.iterrows():
+        gcp = rasterio.control.GroundControlPoint(
+            d.row, d.col, d.x, d.y, d.z, d.id, d.info
+        )
+        gcps.append(gcp)
+    return gcps
+
+
+def geocode_file(input_map, basename, outdir, output_map, format="GTiff"):
+    # Write the map to Geotiff file
+    gs.run_command(
+        "r.out.gdal",
+        input=input_map,
+        output=os.path.join(outdir, output_map),
+        format="GTiff",
+    )
+    # Read as rasterio dataset
+    ds = rasterio.open(os.path.join(outdir, output_map))
+    ds = ds.read()[0]
+    # Read the GCPS dataframe
+    env = gs.gisenv()
+    gcp_base_folder = os.path.join(
+        env["GISDBASE"], env["LOCATION_NAME"], env["MAPSET"], "cell_misc", basename
+    )
+    gcps = read_gcps(gcp_base_folder)
+    # Save the geocoded raster (it will replace the previous intermediate file)
+    geoTs = rasterio.transform.from_gcps(gcps)
+    crs = rasterio.crs.CRS.from_epsg(4326)
+    save_GeoTiff(fn=os.path.join(outdir, output_map), crs=crs, transform=geoTs, mat=ds)
+
+
+def export_to_location(outdir, location, input_map, int_map, env):
+    # Run gdalWarp. This is made to avoid ERROR: Input map is rotated - cannot import
+    output_warp = f"gdalwarp_{int_map}"
+    os.system(
+        f"gdalwarp {os.path.join(outdir,int_map)} {os.path.join(outdir,output_warp)}"
+    )
+
+    gs.warning(_("Switching location"))
+
+    # Create the new location with EPSG:4326, in case it does not exist
+    location_folder = env["GISDBASE"]
+    out_location = os.path.join(location_folder, location)
+    if not os.path.exists(out_location):
+        gs.create_location(
+            env["GISDBASE"],
+            location,
+            epsg=4326,
+            desc="Location created by i.saocom.geocode",
+        )
+
+    gs.run_command("g.mapset", mapset="PERMANENT", location=location)
+    gs.run_command(
+        "r.import", input=os.path.join(outdir, output_warp), output=f"{input_map}"
+    )
+    os.remove(os.path.join(outdir, output_warp))
+
+
+def main():
+    input_map = options["map"]
+    data = options["data"]
+    zip_v = options["is_zip"]
+    pols = options["polarizations"]
+    multilook = options["multilook"]
+    basename = options["basename"]
+    location = options["location"]
+    outdir = options["dbase"]
+    env = gs.gisenv()
+
+    if not input_map and not data:
+        gs.fatal(_("Either one of input map or data folder/zip must be specified"))
+
+    if input_map and data:
+        gs.fatal(
+            _("Either one of input map or data folder/zip must be specified, not both")
+        )
+
+    if not input_map and data:
+        # Import real and imaginary bands to a temporary location and geocode them to external file
+        gs.message(_("Running i.saocom.import"))
+        gs.create_location(
+            env["GISDBASE"], f"{basename}_XY_tempLocation", overwrite=1
+        )
+        gs.run_command(
+            "g.mapset", mapset="PERMANENT", location=f"{basename}_XY_tempLocation"
+        )
+        gs.run_command(
+            "i.saocom.import",
+            data=data,
+            is_zip=zip_v,
+            pols=pols,
+            multilook=multilook,
+            basename=basename,
+        )
+        # Get the list of maps to be geocoded
+        map_list = gs.list_grouped(type="raster", pattern=f"{basename}*")[
+            "PERMANENT"
+        ]
+        for m in map_list:
+            gs.run_command("g.region", raster=m)
+            geocode_file(
+                input_map=m,
+                basename=basename,
+                outdir=outdir,
+                output_map=f"{m}.tif",
+                format="GTiff",
+            )
+            # Import the geocoded bands into the target location
+            export_to_location(
+                outdir=outdir,
+                location=location,
+                input_map=f"{m}_geo",
+                int_map=f"{m}.tif",
+                env=env,
+            )
+            # Remove intermediate files
+            os.remove(os.path.join(outdir, f"{m}.tif"))
+            # Go back to temporary location to continue exporting
+            gs.run_command(
+                "g.mapset", mapset="PERMANENT", location=f"{basename}_XY_tempLocation"
+            )
+        # Go back to original location
+        gs.run_command(
+            "g.mapset", mapset=env["MAPSET"], location=env["LOCATION_NAME"]
+        )
+        shutil.rmtree(os.path.join(env["GISDBASE"], f"{basename}_XY_tempLocation"))
+
+    if input_map and not data:
+        # Check if this an XY location
+        proj = gs.read_command("g.proj", flags="g").split("=")[1].split("\n")[0]
+        if proj != "xy_location_unprojected":
+            raise ValueError(
+                "Current location is not XY unprojected (radar coordinates)"
+            )
+        geocode_file(
+            input_map=input_map,
+            basename=basename,
+            outdir=outdir,
+            output_map=f"{input_map}.tif",
+            format="GTiff",
+        )
+        export_to_location(
+            outdir=outdir,
+            location=location,
+            input_map=f"{input_map}_geo",
+            int_map=f"{input_map}.tif",
+            env=env,
+        )
+        # Remove intermediate files
+        os.remove(os.path.join(outdir, f"{input_map}.tif"))
+        # Go back to original location
+        gs.run_command(
+            "g.mapset", mapset=env["MAPSET"], location=env["LOCATION_NAME"]
+        )
+
+
+if __name__ == "__main__":
+    options, flags = gs.parser()
+    main()