-
Notifications
You must be signed in to change notification settings - Fork 148
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add i.saocom and i.sar toolsets #949
base: grass8
Are you sure you want to change the base?
Changes from 13 commits
745f27c
6d6f61a
d3d1b7c
43a8502
508ce86
17b21ec
a5dea51
30a4eca
11c0383
8cc0212
808e70d
624e557
c97825a
b908143
8db88bf
1de88dd
94ddb1c
5c778ec
8f8bc69
42cca39
2a36fb6
f5d33eb
5c1b994
05257d1
9317ea4
397d673
a9aee3a
52331ba
6d97572
d0d78bc
4e91227
e7a89b4
64e96d8
5114820
62c5269
0442071
bd0e1e0
f4602f6
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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/ |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,7 @@ | ||
MODULE_TOPDIR = ../.. | ||
|
||
PGM = i.saocom.geocode | ||
|
||
include $(MODULE_TOPDIR)/include/Make/Script.make | ||
|
||
default: script |
Original file line number | Diff line number | Diff line change | ||||||||
---|---|---|---|---|---|---|---|---|---|---|
@@ -0,0 +1,49 @@ | ||||||||||
<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. | ||||||||||
|
||||||||||
<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' ) | ||||||||||
Comment on lines
+35
to
+36
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||||||
</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>REFERENCES</h2> | ||||||||||
|
||||||||||
Reference or references to paper related to the tool or references which algorithm the tool is based on. | ||||||||||
sase1988 marked this conversation as resolved.
Show resolved
Hide resolved
|
||||||||||
|
||||||||||
<h2>SEE ALSO</h2> | ||||||||||
|
||||||||||
<a href="../i.saocom.import/i.saocom.import.html">i.saocom.import</a>, | ||||||||||
<a href="../../i.sar/i.sar.amplitude/i.sar.amplitude.html">i.sar.amplitude</a>, | ||||||||||
<a href="../../i.sar/i.sar.speckle/i.sar.speckle.html">i.sar.speckle</a>, | ||||||||||
<a href="../../i.sar/i.sar.pauliRGB/i.sar.pauliRGB.html">i.sar.pauliRGB</a> | ||||||||||
sase1988 marked this conversation as resolved.
Show resolved
Hide resolved
|
||||||||||
|
||||||||||
<h2>AUTHORS</h2> | ||||||||||
|
||||||||||
Santiago Seppi, CONAE, Argentina. |
Original file line number | Diff line number | Diff line change | ||||
---|---|---|---|---|---|---|
@@ -0,0 +1,303 @@ | ||||||
#!/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. | ||||||
# | ||||||
# | ||||||
############################################################################# | ||||||
Comment on lines
+3
to
+18
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. check indentation |
||||||
|
||||||
# %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 | ||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
|
||||||
# %end | ||||||
# %option | ||||||
# % key: pols | ||||||
# % 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']) | ||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Shouldn't the default answer be expressed as "plain grass", i.e., 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 grass | ||||||
from grass.pygrass.modules.shortcuts import general as g | ||||||
from grass.pygrass.modules.shortcuts import raster as r | ||||||
from zipfile import ZipFile | ||||||
import rasterio | ||||||
from rasterio.mask import mask | ||||||
from rasterio.vrt import WarpedVRT | ||||||
import geopandas as gpd | ||||||
import numpy as np | ||||||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
Suggested change
It is already imported above There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Perhaps rasterio, pandas, geopandas (maybe others too) could be added as lazy loading to allow for module installation even when they are not installed. See for example:
|
||||||
from xml.etree import ElementTree as ET | ||||||
from osgeo import gdal | ||||||
from affine import Affine | ||||||
import pandas as pd | ||||||
import shutil | ||||||
from grass.exceptions import ParameterError | ||||||
|
||||||
|
||||||
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 | ||||||
grass.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 = grass.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)}" | ||||||
) | ||||||
|
||||||
grass.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): | ||||||
grass.create_location( | ||||||
env["GISDBASE"], | ||||||
location, | ||||||
epsg=4326, | ||||||
desc="Location created by i.saocom.geocode", | ||||||
) | ||||||
|
||||||
grass.run_command("g.mapset", mapset="PERMANENT", location=location) | ||||||
grass.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["pols"] | ||||||
multilook = options["multilook"] | ||||||
basename = options["basename"] | ||||||
location = options["location"] | ||||||
outdir = options["dbase"] | ||||||
env = grass.gisenv() | ||||||
|
||||||
if not input_map and not data: | ||||||
grass.fatal(_("Either one of input map or data folder/zip must be specified")) | ||||||
|
||||||
if input_map and data: | ||||||
grass.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 | ||||||
grass.message(_("Running i.saocom.import")) | ||||||
grass.create_location( | ||||||
env["GISDBASE"], f"{basename}_XY_tempLocation", overwrite=1 | ||||||
) | ||||||
grass.run_command( | ||||||
"g.mapset", mapset="PERMANENT", location=f"{basename}_XY_tempLocation" | ||||||
) | ||||||
grass.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 = grass.list_grouped(type="raster", pattern=f"{basename}*")[ | ||||||
"PERMANENT" | ||||||
] | ||||||
for m in map_list: | ||||||
grass.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 | ||||||
grass.run_command( | ||||||
"g.mapset", mapset="PERMANENT", location=f"{basename}_XY_tempLocation" | ||||||
) | ||||||
# Go back to original location | ||||||
grass.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 = grass.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 | ||||||
grass.run_command( | ||||||
"g.mapset", mapset=env["MAPSET"], location=env["LOCATION_NAME"] | ||||||
) | ||||||
|
||||||
|
||||||
if __name__ == "__main__": | ||||||
options, flags = grass.parser() | ||||||
main() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.