diff --git a/README.md b/README.md
new file mode 100644
index 0000000..b76d2e5
--- /dev/null
+++ b/README.md
@@ -0,0 +1,42 @@
+# Overview
+## Summary
+This repository contains GRASS GIS Addons, bash scripts, and a detailed step-by-step documentation of a GRASS GIS based georectification process of satellite images from the KEYHOLE reconnaissance missions KH-4A/B (CORONA), KH-7 (GAMBIT), and KH-9 (HEXAGON). Further background can be found [here](https://www.mundialis.de/en/georeferenzierung-von-corona-spionagesatellitendaten/).
+## Repo Content
+ - **documentation**: contains a detailed step-by-step guide to rectify KEYHOLE reconnaissance missions based on GRASS GIS. See [documentation/README.md](documentation/README.md) to get started.
+ - **grass_addons**: contains two addons for GRASS GIS that are required for the rectification process
+ - **scripts**: contains shell scripts that are required for the rectification process
+
+
+ ## KEYHOLE missions overview:
+
+ ##### KH-4A, KH-4B - CORONA
+
+ - Panorama camera
+     - -> orthorectification with panorama correction
+ - Film size: 70 mm x 756.9 mm
+ - flight direction seems to be south
+ - in the scanned film, north is approximately up
+
+ ##### KH-7 - GAMBIT
+
+ - strip camera, acquired imagery in continuous lengthwise sweeps of the terrain.
+     - -> orthorectification does not really work, bad geolocation
+     - -> panorama correction not applicable
+     - -> standard rectification
+ - Film size: 9 inch (228.6 mm) x variable length
+ - flight direction seems to be south
+ - in the scanned film, north is approximately right
+
+ ##### KH-9 - HEXAGON
+
+ - Panorama camera
+     -> orthorectification with panorama correction
+ - flight direction seems to be south
+ - in the scanned film, north is approximately up
+ - focal_length: 152.4 cm
+
+ ##### Further KEYHOLE missions
+
+ - KH-6 Lanyard: Lanyard was a failed mission - no data with acceptable quality available
+ - KH-5 Argon: ground resolution too low - irrelevant
+ - KH-8 Gambit 3: no data available
diff --git a/documentation/0_preparation.md b/documentation/0_preparation.md
new file mode 100644
index 0000000..6e93356
--- /dev/null
+++ b/documentation/0_preparation.md
@@ -0,0 +1,70 @@
+# Preparation
+## Operating System
+The described process was developed and tested using Linux Mint 20. It is in general optimized for use with Linux. No tests on other operating systems were performed.
+
+## Required Software
+- [GRASS GIS >= v7.9](https://grasswiki.osgeo.org/wiki/Installation_Guide)
+- [QGIS](https://www.qgis.org/en/site/forusers/alldownloads.html#)
+- (optional): [Hugin](https://hugin.sourceforge.io/) for the automatic stitching process:
+    - for Ubuntu: use `apt-get install`
+            --> `apt-get install hugin`
+    - for Fedora: use `dnf install`
+            --> `dnf install hugin`
+
+- GRASS GIS Addons:
+    - Start GRASS GIS and install the GRASS modules [i.ortho.position](../grass_addons/i.ortho.position) and [i.ortho.corona](../grass_addons/i.ortho.corona) from this repository:
+    ```
+    g.extension extension=i.ortho.position url=/path/to/this/repo/grass_addons/i.ortho.position
+    g.extension extension=i.ortho.corona url=/path/to/this/repo/grass_addons/i.ortho.corona
+    ```
+    - install the [r.in.nasadem](https://grass.osgeo.org/grass82/manuals/addons/r.in.nasadem.html) addon from the official GRASS Addons repository:
+    ```
+    g.extension extension=r.in.nasadem
+    ```
+    To use this addon, you will need a registered user account at [NASA EARTHDATA](https://urs.earthdata.nasa.gov/)
+
+
+## Data preparation
+**0)** Download the raw Keyhole Scene you would like to process and the vector file containing its rough extent and metadata from the USGS Earth Explorer.
+- if the vector dataset is not available yet:
+  - Browse https://earthexplorer.usgs.gov/ and navigate to the scene of interest using the ENTITY ID (in Step 2: "Data Sets" select Declass I for CORONA, Declass II for GAMBIT or Declass III for HEXAGON; enter the ENTITY ID in Step 3: "Additional Criteria"). Click at "Click here to export your results" (above the search results) to export a .zip-file with the vector dataset of the scene
+  - Note: For some HEXAGON scenes there is no specific vector dataset that can be downloaded. Instead the .zip-file will contain footprints for all HEXAGON scenes around the world. To cut out only the desired footprint you need to open the attribute table of this dataset in e.g. QGIS, search for the desired footprint using the ENTITY ID of the scene and create a new layer that only contains this footprint. Save this layer as a vector dataset and continue at the step above.  
+
+**1)** Launch GRASS GIS and create a GRASS xy location, e.g. corona_xy (you can name it e.g. gambit_xy for GAMBIT scenes - Note: in the walkthrough below the location is always called corona_xy)
+- xy means no CRS, no georeferencing information as the original USGS scans do not have any CRS information
+
+**2)** create a GRASS location in the projection of the target UTM zone
+- example for zone 32N:
+    `grass -c epsg:32632 ~/grassdata/corona_utm32n`
+
+**3)** download/import the DEM (Note: only applies to KH-4A/KH-4B CORONA and KH-9 HEXAGON scenes)
+- if the DEM is already available, import it into the UTM location using `r.import input=/path/to/the/DEM.tif output=dem_<country>`
+- if no DEM is available yet:
+   - import the vector dataset of the rough scene extent using `v.import`: `v.import input=/path/to/scene_extent.shp output=<scene>_extent`
+   - set the region to this vector and add some km of buffer around it (due to inaccuracies of the footprints):
+    `g.region vector=<scene>_extent n=n+60000 s=s-60000 w=w-60000 e=e+60000 res=30 -pa`
+   - import the nasadem:
+    `r.in.nasadem user=<USER> password=<PASSWORD> output=dem_<COUNTRY> memory=6000 resolution=30 method=bilinear`
+Note: The NASA only provides DEMs that cover an area of the Earth that consists at least proportionally of land mass. All parts of the earth that lie in water are not covered by the DEMs. Due to the enlargement of the region by 60km, it may happen that `r.in.nasadem` wants to download DEMs that are completely in water, especially when the area of interest is a coastal region - these DEMs are not available! Consequently, `r.in.nasadem` will raise an error. In order to get the DEMs anyway, either the region in the corresponding cardinal direction must be set smaller than 60km, or the DEMs must be downloaded manually for the corresponding area under https://search.earthdata.nasa.gov/.
+   - set DEM NULLs to 0 to avoid nodata in oceans:
+    `r.null dem_<COUNTRY> null=0`
+   - round the DEM to integer values to save disk space:
+    ` r.mapcalc expression="dem_<country> = round(dem_<country>)" --o`
+
+Note: In the walkthrough below, it is assumed that the target UTM Zone is 32N.
+For other zones, this has to be adapted in steps 2_gcp_collection, 3_georectification and 3_orthorectification
+
+**4)** collect scene specifications
+
+A scene specification document is not required as input, but it is helpful to keep track of the processing in it.
+1. Get the four approximate corner coordinates.
+    - Use the scene extent vector dataset
+    - OR, check the USGS earth explorer (https://earthexplorer.usgs.gov/):
+        - select Dataset -> Declassifed -> declassI (for CORONA), declassII (for GAMBIT) or declassIII (for HEXAGON)
+        - additional criteria: scene id
+        - get metadata + browse overlay
+2. Check if the image needs to be flipped --> use the USGS earth explorer thumbnails to find out
+3. For CORONA and HEXAGON note the look direction of the camera:
+    - It's indicated in the scene name (F for forward, A for aft)
+
+**Continue** with the stitching of the individual scene scans in [1_stitching.md](1_stitching.md) for manual stitching or [1_automatic_stitching.md](1_automatic_stitching.md) for automatic stitching using the Hugin panorama software
diff --git a/documentation/1_automatic_stitching.md b/documentation/1_automatic_stitching.md
new file mode 100644
index 0000000..583cb89
--- /dev/null
+++ b/documentation/1_automatic_stitching.md
@@ -0,0 +1,33 @@
+# Automatic Stitching
+The Keyhole scenes are delivered in several pieces: CORONA scenes usually in 2-4 pieces, GAMBIT in 2 pieces and HEXAGON in 4-14 pieces. These pieces are not georeferenced and must be stitched together.
+
+Automatic stitching is successfully tested for CORONA (KH-4A/KH-4B) scenes, GAMBIT 1 (KH-7) scenes and HEXAGON (KH-9) scenes with a maximum of 4 pieces. Scenes with a very high proportion of undifferentiable land cover (>80%), e.g. deserts or ocean, or scenes with poor image quality may not be successfully stitched.
+Automatic stitching is particularly suitable for scenes that can also be easily stitched manually and scenes that do not have clearly identifiable infrastructure but do have clearly recognisable land structures such as rock formations.
+
+The automatic stitching is processed with [Hugin](https://hugin.sourceforge.io/).
+
+
+## Requirements
+
+- 30++ GB RAM for CORONA scenes
+- 55++ GB RAM for HEXAGON scenes
+- 65++ GB RAM for GAMBIT scenes
+
+Note: some HEXAGON scenes are delivered in more than 4 pieces. The automatic stitching process for these scenes would possibly need more than 60 GB RAM, which has not been successfully tested yet. Splitting scenes with e.g. 7 pieces in two stitching groups has not been successful too due to the same RAM problem. For those scenes it is possibly practical to enlarge the physical RAM with a swap file to a total of > 80 GB.
+
+## Stitching process
+
+**1)** Create a folder in your local file system that only contains the [hugin_automatic_stitching.sh](../scripts/hugin_automatic_stitching.sh) script and the pieces of the scene that you would like to stitch together. Note: The script has to be adapted for HEXAGON scenes with more or fewer pieces (currently for scene pieces `_a` to `_f`) --> lines 35, 45 in the script
+
+**2)** Navigate to the created folder in your terminal and run the `hugin_automatic_stitching.sh` script passing the name of the satellite mission (`CORONA`, `GAMBIT` or `HEXAGON`) as an argument:
+
+e.g. `bash hugin_automatic_stitching.sh CORONA`
+
+It outputs `stitched_D*.tif` which is the result of the stitching process.
+
+**3)** The stitching process produces a `.tif` that contains 2 bands. Only the first band is relevant for further processing steps. Therefore only import the first band to the GRASS xy location using the `band=1`parameter.
+
+`r.in.gdal input=stitched_D*.tif output=stitched_D* band=1`
+Note: Make sure the output name does not include `-`, otherwise `r.mapcalc` in the next step  will not work. Use the `g.rename` command if required.
+
+**Continue** with the cropping of the remaining physical film border in [2_cropping.md](2_cropping.md).
diff --git a/documentation/1_stitching.md b/documentation/1_stitching.md
new file mode 100644
index 0000000..c34f950
--- /dev/null
+++ b/documentation/1_stitching.md
@@ -0,0 +1,48 @@
+# Stitching
+The Keyhole scenes are delivered in several pieces: CORONA scenes usually in 2-4 pieces, GAMBIT in 2 pieces and HEXAGON in 4-14 pieces. These pieces are not georeferenced and must be stitched together:
+
+**1)** Start GRASS GIS in the xy location that was created in the first step, select the mapset `PERMANENT`
+
+**2)** import the westernmost part using `r.in.gdal` (usually *_d for CORONA, but check with the USGS earth explorer thumbnail):
+   `r.in.gdal in=DS1043-2201DF006_d.tif out=DS1043-2201DF006_d`
+
+**3)** Move the remaining parts to the correct place, while the westermost part stays the reference. For each part after the westernmost, do the following:
+
+ * create a new mapset named after the part, e.g. *_c:
+    `g.mapset -c DS1043-2201DF006_c`
+ * import the part to stitch, e.g. *_c:
+     `r.in.gdal in=DS1043-2201DF006_c.tif out=DS1043-2201DF006_c`
+ * create an image group named after and consisting of the raster, e.g. *_c:
+     `i.group group=DS1043-2201DF006_c in=DS1043-2201DF006_c`
+ * In the GUI: Launch the "Georectify" dialog. source location: corona_xy,
+     source mapset: *_c, e.g. DS1043-2201DF006_c. click next, keep the settings
+     in the next dialog, click next, source-map: the map to be stitched, e.g. *_c in
+     the corresponding mapset. target map: The previous reference map, e.g. *_d in the
+     PERMANENT mapset.
+ * Collect 4-8 GCPs that you find in the overlap area of the two scene pieces, e.g. four at the edge, and two within the images to compensate for rotation. Use distinct pixels.
+ * Under "Georectifier Settings" --> "Georectification", make sure 1.st order is ticked
+ * Tick all GCPs, and click on "Recalculate RMS error". Since we don't have a CRS,
+     the unit of the errors is pixels. It should be < 1.0.
+ * click on "Georectify"
+ * Sometimes the writing to the PERMANENT mapset does not work, if so: copy the "rectified"
+     raster to the PERMANENT mapset and if necessary rename it, e.g.:
+     `g.mapset PERMANENT`
+     `g.copy DS1043-2201DF006_c_georect17884@DS1043-2201DF006_c,DS1043-2201DF006_c`
+ * Load the result in the GUI check if it fits. Therefore, right-click on the reference map, e.g. *_d and choose "Display layer". Then navigate in the GUI to "Layers" (at the bottom) and click "Add Raster Map" and choose the georectified raster.
+ * Delete the temporary mapset in your local file system inside the grassdb
+ * Repeat step **3)** for the next pieces of the scene (e.g. *_b), until all pieces have been moved to the correct position.
+
+**4)** Run the script [raster_patch_xy.sh](../scripts/raster_patch_xy.sh). It patches the pieces together and outputs a raster in the GRASS location/mapset, as well as a GeoTiff and a .vrt with a rough georeferencing.
+* modify in the script the parameters `ORDEREDLIST` (line 21) and `FLIP` (line 28) according to the respective scene. `ORDEREDLIST` has to begin with the part that was stitched last (e.g. *_a, see the comments in the script)
+* run the script by passing the parameters of the scene ID and the corner coordinates:
+
+   `bash raster_patch_xy.sh SCENEID NWLON NWLAT NELON NELAT SWLON SWLAT SELON SELAT`
+
+   e.g.:
+
+   `bash raster_patch_xy.sh DS1043-2201DF006 7.639 36.330 10.720 36.792 7.642 36.116 10.816 36.591`
+
+   (latitude and longitude coordinates can be found in the attribute table of the vector map)
+   Note: KH-7 Gambit scenes may be flipped by rougly 90 degrees so the coordinates passed to the script may not be correct. This however only affects the output .vrt, which is just a visual aid for GCP collection and not crucial for the process.
+
+**Continue** with the cropping of the remaining physical film border in [2_cropping.md](2_cropping.md).
diff --git a/documentation/2_cropping.md b/documentation/2_cropping.md
new file mode 100644
index 0000000..eff4795
--- /dev/null
+++ b/documentation/2_cropping.md
@@ -0,0 +1,32 @@
+# Cropping
+The scenes are delivered with black borders that should be removed for asthetic and functional reasons.
+
+**0)** If you manually stitched the raster: import the stitched raster using `r.import`, e.g.:
+    `r.import in=DS1036-2155DF050_stitched.tif out=DS1036_stitched_raster` Note: Make sure the output name does not include `-`, otherwise `r.mapcalc` in the next step  will not work. Use the `g.rename` command if required.
+
+**1)** Manually set the region in the GUI to the extent of the actual image (only the relevant pixels, no black border etc.):
+
+Therefore right-click on the stitched layer --> "Display Layer" --> "Different Zoom Settings" --> "Set computational region extent interactively" --> Draw a rectangle only around the image data (leave out the black film border)
+
+OR
+
+Use "Display Layer" --> "Analyze Map" --> "Measure Distance" to measure the distance of the four edges of the raster to the start of actual image data, and then set the region accordingly, e.g.: `g.region n=n-4593 s=s+710 e=e-3075 w=w+12865`
+
+**2)** Use `r.mapcalc` to recalculate the stitched raster map. Black areas (with a pixel value of 0) are recalculated to 1 to use the 0 as NoData later on: `r.mapcalc expression="DS1036_stitched_raster_calc = if(DS1036_stitched_raster == 0,1,DS1036_stitched_raster)"`
+
+**3)** `r.mapcalc` changed the colors of the raster map - set it back to grey scale using `r.colors`.
+    `r.colors map=DS1036_stitched_raster_calc color=grey255`
+
+**4)** Run `r.info` to get the basic information about the raster that are used in the next step (rows & columns of the image):
+    `r.info DS1036_stitched_raster_calc`
+
+**5)** Use the values of the rows and columns from the previous step to reset the image coordinates of the raster map so that the image coordinate with the value 0,0 is in the lower left corner of the image. Use the rows value for `north` and the columns value for `east`:
+    `r.region map=DS1036_stitched_raster_calc s=0 n=<ROWS_VALUE> w=0 e=<COLUMNS_VALUE>`
+
+**6)** Set the region to the raster:
+    `g.region raster=DS1036_stitched_raster_calc`
+
+**7)** Export the cropped raster map as a GeoTIFF:
+    `r.out.gdal -cm overviews=5 createopt="COMPRESS=LZW,PREDICTOR=2,TILED=YES" in=DS1036_stitched_raster_calc out=DS1036_stitched_raster_calc.tif`
+
+**Continue** with the selection of GCPs in QGIS in [3_gcp_collection.md](3_gcp_collection.md)
diff --git a/documentation/3.1_gcps_overlapping_scenes.md b/documentation/3.1_gcps_overlapping_scenes.md
new file mode 100644
index 0000000..5fcd00c
--- /dev/null
+++ b/documentation/3.1_gcps_overlapping_scenes.md
@@ -0,0 +1,16 @@
+## GCPS for overlapping scenes
+If you want to rectify scenes from the same overflight that have overlap areas to each other, there is a way to quicker rectifying them than doing it for each scene individually.
+
+Example: You want to rectify the scenes `DS1036-2155DF050`, `DS1036-2155DF051`, `DS1036-2155DF052` and `DS1036-2155DF053`. Instead of collecting GCPs for each scene individually, you can use the overlap area between two scenes to process them quicker. In the following workflow it is assumed that you want to rectify the scenes in the direction from top to bottom (`DS1036-2155DF050` to `DS1036-2155DF053`) - of course it is also possible the other way round.
+
+**1)** Collect GCPs for the first scene (e.g. `DS1036-2155DF050`). Set 10-15 GCPs at the top of the scene and 10-15 at the bottom. The GCPs at the bottom needs to be in an area of the image that overlaps with the next scene (in this example: `DS1036-2155DF051`). Due to little differences in the width of the scenes, which is a result of the flight direction of the satellite, and due to the possibility of poor image quality in the corners of the image you should make sure that you leave a free space at the left and the right border of the scene where you do not set any GCPs (~2,000 pixels at each side of the image).
+
+**2)** Use the collected GCPs to continue the process at the GCP conversion in `3_gcp_collection.md` for the first scene (`DS1036-2155DF050`).
+
+**3)** After finishing the rectification process of the first scene (`DS1036-2155DF050`), you can start off the GCP collection process for the next scene (in this example: `DS1036-2155DF051`) with already half a set of GCPs: Open the file with the collected GCPs from the first scene (the `.points` file) and delete all points that were set in the top area of the image. Save this file and open it in the QGIS-georectify window of the new scene (here: `DS1036-2155DF051`). You now already have the GCPs that belong to the top area of this scene. Those GCPs already have the correct reference coordinates. The x-image coordinates are also almost correct and give a good orientation about the actual position of the points in the image. Only the y-image coordinates are completely wrong. Using the "move" function, the GCPs in the image can now be moved up to their actual position.
+
+**4)** Afterwards, collect further 10-15 GCPs in the lower image area, which are located in the overlap area to the next scene (in this example: `DS1036-2155DF052`). Use the collected GCPs to process the geo- and orthorectification of the second scene (`DS1036-2155DF051`).
+For the third scene the same is done as described in steps 1 to 4.
+Repeat steps 1 to 4 until you have rectified all the required scenes of the same flyover.
+
+By reusing the points from the previous scene, the processing time of the following scene is significantly reduced. If several scenes of the same flyover are to be rectified, this can lead to great time savings with the same quality of results compared to the individual GCP collection for each scene.
diff --git a/documentation/3_gcp_collection.md b/documentation/3_gcp_collection.md
new file mode 100644
index 0000000..12660a7
--- /dev/null
+++ b/documentation/3_gcp_collection.md
@@ -0,0 +1,48 @@
+# GCP Collection
+The collection of Ground Control Points (GCPs) is the most time-consuming step. QGIS is used here as it provides a simpler handling of the GCP collection process.
+
+## Preparation
+**1)** Launch QGIS, activate the GDAL-georeferencing plug-in (one time only). Make sure you are in EPSG:4326 (WGS84)
+
+**2)** Load reference data:
+   - OSM: e.g. use the Terrestris OSM WMS (http://ows.terrestris.de/osm/service)
+   - Google Imagery: In the QGIS Browser, right click XYZ-Tiles, New Connection:
+     [http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}](http://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z})
+   - BING Aerial: In the QGIS Browser, right click XYZ-Tiles, New Connection:
+   [http://ecn.t3.tiles.virtualearth.net/tiles/a{q}.jpeg?g=1](http://ecn.t3.tiles.virtualearth.net/tiles/a{q}.jpeg?g=1)
+   - ESRI World Imagery: In the QGIS Browser, add a connection to ArcGisMapServer:
+     [http://server.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer](http://server.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer)
+
+**3)** Open the Georeferencing tool
+- Open the Georeferencing window ("Raster" --> "Georeferencing"). Load the newly created .tif raster (The result of [2_cropping.md](2_cropping.md))
+- Inside the Georeferencing Window: Open Settings --> Configure Georeferencer --> Click "Show Georeferencer Window docked" to directly have the window of the target Raster (the .tif) under the window of your reference data which makes it easier to quickly collect the GCPs.
+
+## GCP collection
+Use the footprint that was delivered together with the scene pieces by the USGS or the .vrt-file to get a rough idea of where the scene is located.
+Start collecting GCPs:
+Identify a point in the target raster that you would like to have as GCP. Click at this point --> in the new window "Enter Map Coordinates" select "From Map Canvas" --> click at the same point in the reference data to set the reference GCP. Set ~20 GCPs at points where
+the reference data agree well. Try to distribute the GCPs as evenly as possible across the image. Good points are street junctions, because there you can
+compare OSM and imagery data. If no OSM data is available, pick points that agree in the different imagery sources, e.g. river junctions. Avoid mountain tops or ares with
+steep terrain, because the imagery might be a bit shifted or distorted there.
+Setting 20 GCPs for one image takes roughly 1.5 hours. Note that 20 GCPs are usually enough as the improvement in accuracy when using more GCPs is only marginal.
+IMPORTANT: Manually SAVE your progress regularly by clicking "Save GCPs as". Save the
+final GCPs accordingly.  
+
+## Multiple overlapping scenes
+If you have multiple scenes from the same overflight with overlapping areas, there is a possibility to rectify them more quickly than doing this for each scene individually. Therefore look at [3.1_gcps_overlapping_scenes.md](3.1_gcps_overlapping_scenes.md)
+
+## GCP conversion
+The GCPs collected in QGIS have to be adapted to be processable by GRASS, which is done by the script [gcp_qgis2grass.sh](../scripts/gcp_qgis2grass.sh)
+
+**1)** Start GRASS in the location with the projection of the target UTM zone (e.g. corona_utm32).
+
+**2)** Run the script passing the qgis-gcp file and the DEM as arguments, e.g.:
+
+`bash gcp_qgis2grass.sh xyz.points dem_<country> 32`
+
+The last argument is the UTM zone. The script outputs a GRASS-optimized GCP-file with the file-ending ".grass" next to the original file.
+
+## For KH-7 (GAMBIT-1)
+The .grass file has to be manually edited. Open the file in a Text Editor and remove each entry of z-coordinates (presumably the 3rd entry in each row). Save this file as `POINTS`.
+
+**Continue** with the rectification in [4_orthorectification.md](4_orthorectification.md) (for CORONA or HEXAGON) or [4_georectification.md](4_georectification.md) (for GAMBIT)
diff --git a/documentation/4_georectification.md b/documentation/4_georectification.md
new file mode 100644
index 0000000..71dd1bf
--- /dev/null
+++ b/documentation/4_georectification.md
@@ -0,0 +1,42 @@
+# Georectification for KH-7 (GAMBIT-1)
+The final step is the rectification of the scene using GRASS commands.
+Note: This readme applies only to KH-7 (GAMBIT-1). For KH-4A/KH-4B CORONA and KH-9 HEXAGON go to [4_orthorectification.md](4_orthorectification.md).
+
+## Rectification
+**1)** Launch GRASS GIS in the gambit_xy location.
+
+**2)** Create a group in the gambit_xy location. The input raster is the stitched raster:
+    `i.group group=<NAME> in=<NAME_OF_STITCHED_RASTER_MAP>`
+
+**3)** Run i.target to target the group into the gambit UTM location:
+    `i.target group=<NAME> location=gambit_utm33n mapset=PERMANENT`
+
+**4)** Go to your local file system and manually insert the `POINTS`-file into the group-folder in the GRASS database (GRASSDATA/<LOCATIONS>/<MAPSETS>/group/<GROUP>)
+
+**5)** Run `m.transform` to calculate the RMS of the GCPs:
+    `m.transform -s group=<NAME> order=2`
+(As a reference: in a first test the RMS was ~27m)
+
+
+m.transform with `order=2` gives better results but requires evenly distributed GCPs; use `order=1` if evenly distributed GCPS can not be obtained
+
+**6)** Run i.rectify to process the actual georectification:
+    `i.rectify group=<NAME> extension=_rectified order=2 resolution=0.5 method=linear`
+
+
+#### File Saving
+
+**7)** Start GRASS in the location with the projection of the target UTM zone (e.g. gambit_utm33n)
+
+**8)** Set the region to the new georectified raster, e.g.:
+    `g.region raster=DZB00401700059H005001_rectified`
+
+**9)** Round the raster back to integer:
+    `r.mapcalc expression="DZB00401700059H005001_rectified_rounded = round(DZB00401700059H005001_rectified)"`
+
+**10)** Export the result as GeoTiff.
+    Since the raster is very large, activate compression, tiles, and overviews.
+    `export COMPRESS_OVERVIEW="LZW"`
+    `r.out.gdal -cm overviews=5 createopt="COMPRESS=LZW,PREDICTOR=2,TILED=YES,BIGTIFF=YES" in=DZB00401700059H005001_rectified_rounded out=DZB00401700059H005001.tif nodata=0`
+
+**11)** Analyze the GeoTiff in QGIS, e.g. compare it to OSM data.
diff --git a/documentation/4_orthorectification.md b/documentation/4_orthorectification.md
new file mode 100644
index 0000000..1d344f4
--- /dev/null
+++ b/documentation/4_orthorectification.md
@@ -0,0 +1,76 @@
+## Orthorectification for KH-4A/KH-4B CORONA and KH-9 HEXAGON
+The final step is the rectification of the scene using GRASS commands. Note: This readme applies only to KH-4A, KH-4B (CORONA) and KH-9 (HEXAGON). For KH-7 GAMBIT go to [4_georectification.md](4_georectification.md).
+
+For an exhaustive explanation of the underlying processing steps, see the [i.ortho.photo help](https://grass.osgeo.org/grass82/manuals/i.ortho.photo.html).
+
+
+## Preparation & Rectification
+**1)** Launch GRASS GIS in the corona_xy location
+
+**2)** create a group of the raster to be orthorectified, e.g.
+
+   `i.group group=DS1043-2201DF006 in=DS1043-2201DF006`
+Note: Make sure the group name is not too long (not longer than ~20 digits), otherwise an error can occur: `buffer overflow detected`
+
+**3)** use `i.ortho.corona` to check the GCPs and rectify the image
+
+
+  Use `i.ortho.corona` as follows:
+
+  Use the -t flag first to examine the RMS of the GCPs. If the RMS are within limits (total Reverse RMS < 25m), run i.ortho.corona without the -t flag, but use the logfile parameter to save the pointwise error statistics.
+
+**3.1)** Estimate the RMS
+   - use the `i.ortho.corona --help` command to get a better understanding of the individual parameters
+
+   - for CORONA (KH-4A/KH-4B), the following parameters have to be adapted:
+  ```
+   i.ortho.corona -t group=<GROUP> grasspoints=path/to/adapted/grasspoints nw=<NW_COORDINATES> ne=<NE_COORDINATES> sw=<SW_COORDINATES> se=<SE_COORDINATES> targetdem=<DEM> targetlocmapset=<MAPSET_WITH_PROJECTION_OF_TARGET_UTM_ZONE> proj=<PROJ_STRING_WITH_UTM_ZONE> camheading=<CAMHEADING> flightdir=<FLIGHT_DIRECTION> target_res=2 memory=5000
+
+    e.g.:
+   i.ortho.corona -t group=DS1043-2201DF006 grasspoints=path/to/adapted/gcps nw=7.639,36.330 ne=10.720,36.792 sw=7.642,36.116 se=10.816,36.591 targetdem=nasadem_tun_alg_complete targetlocmapset=corona_utm32/PERMANENT proj="+proj=utm +zone=32 +datum=WGS84 +units=m" camheading=forward flightdir=south target_res=2 memory=5000
+   ```
+   - for HEXAGON (KH-9), the following parameters have to be adapted:
+   ```
+   i.ortho.corona -t group=<GROUP> grasspoints=path/to/adapted/grasspoints nw=<NW_COORDINATES> ne=<NE_COORDINATES> sw=<SW_COORDINATES> se=<SE_COORDINATES> targetdem=<DEM> targetlocmapset=<MAPSET_WITH_PROJECTION_OF_TARGET_UTM_ZONE>
+   proj=<PROJ_STRING_WITH_UTM_ZONE> camheading=<CAMHEADING> flightdir=<FLIGHT_DIRECTION> target_res=1 memory=5000 focal_length=1052.4
+
+   e.g.
+   i.ortho.corona -t group=D3C1207-400442F004 grasspoints=path/to/adapted/gcps nw=7.639,36.330 ne=10.720,36.792 sw=7.642,36.116 se=10.816,36.591 targetdem=nasadem_tun_alg_complete targetlocmapset=hexagon_utm32/PERMANENT
+   proj=+proj=utm +zone=32 +datum=WGS84 +units=m camheading=forward flightdir=south target_res=1 memory=5000 focal_length=1052.4
+   ```
+   - If the error attributed to individual GCPs is too high, deactivate them in the points file in `grassdata/corona_xy/PERMANENT/group/>GROUPNAME</CONTROL POINTS`. Change the 1 in the last column of the respective point to a 0 and rerun `i.ortho.corona -t ...`
+
+
+**3.2)** Start the actual rectification and indicate the logfile
+
+   - Uses the same command as in **3.1**, but remove -t flag and add `logfile=~/logfile` at the end of the command and add `target_res=1` for HEXAGON:
+
+   - for CORONA (KH-4A/KH-4B):
+   ```
+   i.ortho.corona group=<GROUP> grasspoints=path/to/adapted/grasspoints nw=<NW_COORDINATES> ne=<NE_COORDINATES> sw=<SW_COORDINATES> se=<SE_COORDINATES> targetdem=<DEM> targetlocmapset=<MAPSET_WITH_PROJECTION_OF_TARGET_UTM_ZONE> proj=<PROJ_STRING_WITH_UTM_ZONE> camheading=<CAMHEADING> flightdir=<FLIGHT_DIRECTION> target_res=2 memory=5000 logfile=~/logfile
+   ```
+   - for HEXAGON (KH-9):
+   ```
+   i.ortho.corona group=<GROUP> grasspoints=path/to/adapted/grasspoints nw=<NW_COORDINATES> ne=<NE_COORDINATES> sw=<SW_COORDINATES> se=<SE_COORDINATES> targetdem=<DEM> targetlocmapset=<MAPSET_WITH_PROJECTION_OF_TARGET_UTM_ZONE>
+   proj=<PROJ_STRING_WITH_UTM_ZONE> camheading=<CAMHEADING> flightdir=<FLIGHT_DIRECTION> target_res=1 memory=5000 focal_length=1052.4 logfile=~/logfile
+   ```
+
+Note: If the used imagery has a different scan resolution than 1432 pixels/cm, this can be adapted using the `scan_res` parameter.
+
+`i.ortho.corona` without the `-t` flag will run the rectification and create a result raster in the indicated target location/mapset.
+
+## File Saving
+
+**4)** Start GRASS in the location with the projection of the target UTM zone (e.g. corona_utm32)
+
+**5)** Set the region to the new orthorectified raster, e.g.:
+
+   `g.region raster=DS1043-2201DF006_utm_v1`
+
+**6)** Export the result as GeoTiff. Since the raster is very large, activate compression, tiles, and overviews.
+
+   `export COMPRESS_OVERVIEW="LZW"`
+
+   `r.out.gdal -cm overviews=5 createopt="COMPRESS=LZW,PREDICTOR=2,TILED=YES,BIGTIFF=YES" in=DS1043-2201DF006_utm_v1 out=DS1043_2201DF006_utm.tif nodata=0`
+
+**7)** Analyze the GeoTiff in QGIS, e.g. compare it to OSM data.
diff --git a/documentation/README.md b/documentation/README.md
new file mode 100644
index 0000000..1463e1c
--- /dev/null
+++ b/documentation/README.md
@@ -0,0 +1,12 @@
+# CORONA / HEXAGON / GAMBIT processing overview
+This directory contains detailed READMEs for each step of the rectification process. See below for a short description:
+
+- [0_preparation.md](0_preparation.md): setup of required software and data
+- [1_stitching.md](1_stitching.md): stitching of the 2 (or more) partial GeoTIFF scenes into one single GeoTIFF (applicable for CORONA; HEXAGON and GAMBIT scenes)
+    OR
+- [1_automatic_stitching.md](1_automatic_stitching.md): automatic stitching for CORONA and HEXAGON scenes with a maximum of 4-5 scene pieces (scenes with more pieces will need >60GB RAM)
+- [2_cropping.md](2_cropping.md): crop the raster to an extent where only relevant image pixels are left and residual pixels from the physical film scan are removed
+- [3_gcp_collection.md](3_gcp_collection.md): identification of GCPs between unrectified scene and reference data
+- [3.1_gcps_overlapping_scenes.md](3.1_gcps_overlapping_scenes.md): additional hints for GCP collection of an entire swath of scenes
+- [4_orthorectification.md](4_orthorectification.md): orthorectification of CORONA and HEXAGON scenes (incl. correction of panoramic distortion)
+- [4_georectification.md](4_georectification.md): georectification of GAMBIT scenes (no correction of panoramic distortion required)
diff --git a/grass_addons/README.md b/grass_addons/README.md
new file mode 100644
index 0000000..b07b982
--- /dev/null
+++ b/grass_addons/README.md
@@ -0,0 +1,13 @@
+## Overview
+
+- i.ortho.position estimates the camera postion from the corner coordinates of the scene in longitudes and latitudes, omega and flight height.
+- i.ortho.corona runs the required steps from i.ortho.photo for the orthorectification of KH-4A/B Corona and KH-9 HEXAGON scenes.
+
+## Installation
+
+Install these addons in GRASS GIS (https://grass.osgeo.org/) with
+
+```
+g.extension extension=i.ortho.position url=/path/to/i.ortho.position
+g.extension extension=i.ortho.corona url=/path/to/i.ortho.corona
+```
diff --git a/grass_addons/i.ortho.corona/Makefile b/grass_addons/i.ortho.corona/Makefile
new file mode 100644
index 0000000..7cd0508
--- /dev/null
+++ b/grass_addons/i.ortho.corona/Makefile
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.ortho.corona
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
diff --git a/grass_addons/i.ortho.corona/i.ortho.corona.html b/grass_addons/i.ortho.corona/i.ortho.corona.html
new file mode 100644
index 0000000..dc8dd5a
--- /dev/null
+++ b/grass_addons/i.ortho.corona/i.ortho.corona.html
@@ -0,0 +1,41 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.ortho.corona</em>  runs the required steps from
+<a href="i.ortho.photo.html">i.ortho.photo</a> for the
+orthorectification of KH-4A/B Corona and KH-9 Hexagon scenes.
+
+<p>
+<em>i.ortho.corona</em> requires the image group to be already created. Based on
+the input parameters it then runs <a href="i.ortho.target.html">i.ortho.target</a>
+, <a href="i.ortho.elev.html">i.ortho.elev</a>, <a href="i.ortho.camera.html">i.ortho.camera</a>
+, <a href="i.ortho.position.html">i.ortho.position</a>, <a href="i.ortho.init.html">i.ortho.init</a>,
+<a href="i.ortho.transform.html">i.ortho.transform</a>, and <a href="i.ortho.rectify.html">i.ortho.rectify</a>.
+If the <b> -t </b> flag is applied, it only outputs the result of
+<a href="i.ortho.transform.html">i.ortho.transform</a> and does not continue to
+run <a href="i.ortho.rectify.html">i.ortho.rectify</a>. If the <b>logfile</b> parameter is indicated,
+the output of <a href="i.ortho.tranform.html">i.ortho.transform -t</a>
+will be written to a logfile.
+<p>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.ortho.photo.html">i.ortho.photo</a>,
+<a href="i.ortho.target.html">i.ortho.target</a>,
+<a href="i.ortho.elev.html">i.ortho.elev</a>,
+<a href="i.ortho.camera.html">i.ortho.camera</a>,
+<a href="i.ortho.position.html">i.ortho.position</a>,
+<a href="i.ortho.init.html">i.ortho.init</a>,
+<a href="i.ortho.transform.html">i.ortho.transform</a>,
+<a href="i.ortho.rectify.html">i.ortho.rectify</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Guido Riembauer, <a href="https://www.mundialis.de/">mundialis</a>
+
+<!--
+<p>
+<i>Last changed: $Date$</i>
+-->
diff --git a/grass_addons/i.ortho.corona/i.ortho.corona.py b/grass_addons/i.ortho.corona/i.ortho.corona.py
new file mode 100644
index 0000000..c0f68f1
--- /dev/null
+++ b/grass_addons/i.ortho.corona/i.ortho.corona.py
@@ -0,0 +1,446 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+############################################################################
+#
+# MODULE:     i.ortho.corona
+# AUTHOR(S):  Guido Riembauer
+# PURPOSE:    Runs steps for i.ortho.photo for KH-4A/B Corona and KH-9 Hexagon scenes
+#
+# COPYRIGHT:  (C) 2020 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: Runs steps for i.ortho.photo for KH-4A/B Corona and KH-9 Hexagon scenes
+# % keyword: imagery
+# % keyword: orthorectification
+# %end
+
+# %option G_OPT_I_GROUP
+# % key: group
+# % guisection: Input
+# % required: yes
+# %end
+
+# %option G_OPT_F_INPUT
+# % key: grasspoints
+# % label: path to GCP file in grass format
+# % description: path to GCP file in grass format
+# % guisection: Input
+# % required: yes
+# %end
+
+# %option G_OPT_M_COORDS
+# % key: nw
+# % label: North-West corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+
+# %option G_OPT_M_COORDS
+# % key: ne
+# % label: North-East corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+
+# %option G_OPT_M_COORDS
+# % key: sw
+# % label: South-West corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+
+# %option G_OPT_M_COORDS
+# % key: se
+# % label: South-East corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+
+# %option
+# % key: proj
+# % type: string
+# % label: PROJ definition of the target CRS
+# % description: e.g. +proj=utm +zone=32 +datum=WGS84 +units=m
+# % required: no
+# % answer: +proj=utm +zone=32 +datum=WGS84 +units=m
+# %end
+
+# %option
+# % key: camheading
+# % type: string
+# % label: camera heading (forward or aft)
+# % required: yes
+# % options: forward,aft
+# % answer: forward
+# %end
+
+# %option
+# % key: flightdir
+# % type: string
+# % label: flight direction (north or south)
+# % required: yes
+# % options: north,south
+# % answer: south
+# %end
+
+# %option
+# % key: targetdem
+# % type: string
+# % label: Name of DEM in target location/mapset
+# % required: yes
+# %end
+
+# %option
+# % key: targetlocmapset
+# % type: string
+# % label: target location and mapset in the format location/mapset
+# % required: yes
+# %end
+
+# %option G_OPT_MEMORYMB
+# %end
+
+# %option
+# % key: map_extension
+# % type: string
+# % label: extension to be added to the rectified map in the target mapset
+# % required: no
+# % answer: _utm_v1
+# %end
+
+# %option
+# % key: target_res
+# % type: string
+# % label: resolution of rectified map (in meters)
+# % required: no
+# % answer: 2
+# %end
+
+# %option
+# % key: focal_length
+# % type: string
+# % label: focal length of the camera
+# % required: no
+# % answer: 609.602
+# %end
+
+# %option
+# % key: scan_res
+# % type: string
+# % label: scan resolution of films in pixels/cm
+# % required: no
+# % answer: 1432.3507
+# %end
+
+# %option
+# % key: default_height
+# % type: string
+# % label: altitude of the satellite in meters
+# % required: no
+# % answer: 160000
+# %end
+
+# %option
+# % key: logfile
+# % type: string
+# % label: path to store logfile of i.ortho.transform
+# % required: no
+# %end
+
+# %flag
+# % key: t
+# % description: Only run until i.ortho.transform to analyse GCPs - don't run rectification
+# %end
+
+# import library
+import os
+import sys
+import shutil
+import psutil
+import grass.script as grass
+import subprocess
+
+
+def freeRAM(unit, percent=100):
+    """ The function gives the amount of the percentages of the installed RAM.
+    Args:
+        unit(string): 'GB' or 'MB'
+        percent(int): number of percent which shoud be used of the free RAM
+                      default 100%
+    Returns:
+        memory_MB_percent/memory_GB_percent(int): percent of the free RAM in
+                                                  MB or GB
+
+    """
+    # use psutil cause of alpine busybox free version for RAM/SWAP usage
+    mem_available = psutil.virtual_memory().available
+    swap_free = psutil.swap_memory().free
+    memory_GB = (mem_available + swap_free) / 1024.0 ** 3
+    memory_MB = (mem_available + swap_free) / 1024.0 ** 2
+
+    if unit == "MB":
+        memory_MB_percent = memory_MB * percent / 100.0
+        return int(round(memory_MB_percent))
+    elif unit == "GB":
+        memory_GB_percent = memory_GB * percent / 100.0
+        return int(round(memory_GB_percent))
+    else:
+        grass.fatal("Memory unit <%s> not supported" % unit)
+
+
+def getfiducials(rows, cols, scan_resolution):
+    height_mm = rows * 10 / scan_resolution
+    width_mm = cols * 10 / scan_resolution
+    fid_nw = -width_mm / 2, height_mm / 2
+    fid_ne = width_mm / 2, height_mm / 2
+    fid_se = width_mm / 2, -height_mm / 2
+    fid_sw = -width_mm / 2, -height_mm / 2
+    return fid_nw, fid_ne, fid_se, fid_sw
+
+
+def string_refpoints(fid_nw, fid_ne, fid_se, fid_sw, rows, cols):
+    string = (
+        "# Ground Control Points File\n#\n# target location: XY\n# target mapset: "
+        + "source_mapset\n# source  target  status\n# east north east north "
+        + "(1=ok, 0=ignore)\n#-------------------------------------------------------------\n"
+        + "0 %s     %s %s     1\n" % (rows, fid_nw[0], fid_nw[1])
+        + "%s %s     %s %s     1\n" % (cols, rows, fid_ne[0], fid_ne[1])
+        + "%s 0     %s %s     1\n" % (cols, fid_se[0], fid_se[1])
+        + "0 0     %s %s     1\n" % (fid_sw[0], fid_sw[1])
+    )
+    return string
+
+
+def main():
+
+    # input parameters
+    group = options["group"]
+    grasspoints = options["grasspoints"]
+    nw = options["nw"]
+    ne = options["ne"]
+    sw = options["sw"]
+    se = options["se"]
+    proj = options["proj"]
+    cameraheading = options["camheading"]
+    flightdirection = options["flightdir"]
+    targetloc = options["targetlocmapset"].split("/")[0]
+    targetmapset = options["targetlocmapset"].split("/")[1]
+    targetdem = options["targetdem"]
+    memory = int(options["memory"])
+    map_extension = options["map_extension"]
+    target_resolution = options["target_res"]
+    focal_length = options["focal_length"]
+    scan_resolution = float(options["scan_res"])
+    default_height = options["default_height"]
+    transform_only = flags["t"]
+
+    # get omega
+    if cameraheading == "forward":
+        omega = 15
+    elif cameraheading == "aft":
+        omega = -15
+    if flightdirection == "south":
+        omega = -omega
+    omega = str(omega)
+
+    # get current GRASS environment
+    env = grass.gisenv()
+    start_gisdbase = env["GISDBASE"]
+    start_location = env["LOCATION_NAME"]
+    start_mapset = env["MAPSET"]
+
+    # check for i.ortho.position
+    if not grass.find_program("i.ortho.position", "--help"):
+        grass.fatal(
+            _("The 'i.ortho.position' module was not found, install it first:")
+            + "\n"
+            + "g.extension i.ortho.position url=/path/to/i.ortho.position"
+        )
+
+    # check for free ram
+    free_ram = freeRAM("MB", 100)
+    if free_ram < memory:
+        grass.warning(
+            "Using %d MB but only %d MB RAM available. Setting memory to free RAM"
+            % (memory, free_ram)
+        )
+        memory = free_ram
+
+    # get name of the raster to be rectified
+    rastername_keys = list(
+        grass.parse_command("i.group", group=group, flags="lg").keys()
+    )
+    if len(rastername_keys) > 1:
+        grass.fatal(_("Group contains more than one raster, Exiting..."))
+    rastername = rastername_keys[0]
+
+    # get raster size
+    raster_rows = grass.parse_command("r.info", map=rastername, flags="g")["rows"]
+    raster_cols = grass.parse_command("r.info", map=rastername, flags="g")["cols"]
+
+    # run i.ortho.target
+    grass.run_command(
+        "i.ortho.target",
+        group=group,
+        target_location=targetloc,
+        mapset_loc=targetmapset,
+    )
+
+    # run i.ortho.elev
+    grass.run_command(
+        "i.ortho.elev",
+        group=group,
+        elev=targetdem,
+        location=targetloc,
+        mapset=targetmapset,
+        units="meters",
+    )
+
+    # get fiducials
+    fid_nw, fid_ne, fid_se, fid_sw = getfiducials(
+        int(raster_rows), int(raster_cols), scan_resolution
+    )
+
+    fidstring = "%s,%s,%s,%s,%s,%s,%s,%s" % (
+        fid_nw[0],
+        fid_nw[1],
+        fid_ne[0],
+        fid_ne[1],
+        fid_se[0],
+        fid_se[1],
+        fid_sw[0],
+        fid_sw[1],
+    )
+    # run i.ortho.camera
+    grass.run_command(
+        "i.ortho.camera",
+        group=group,
+        camera=group,
+        name=group,
+        id=group,
+        clf=focal_length,
+        pp="0,0",
+        fid=fidstring,
+    )
+
+    # write REF_POINTS file
+    ref_points_string = string_refpoints(
+        fid_nw, fid_ne, fid_se, fid_sw, raster_rows, raster_cols
+    )
+    ref_points_path = "%s/%s/%s/group/%s/REF_POINTS" % (
+        start_gisdbase,
+        start_location,
+        start_mapset,
+        group,
+    )
+
+    if os.path.isfile(ref_points_path):
+        grass.message(_("REF_POINTS file already in group folder."))
+    else:
+        try:
+            with open(ref_points_path, "w") as file:
+                file.write(ref_points_string)
+        except Exception as e:
+            grass.fatal(_("Unable to write REF_POINTS file."))
+
+    # copy the GCP file into the group folder
+    control_points_path = "%s/%s/%s/group/%s/CONTROL_POINTS" % (
+        start_gisdbase,
+        start_location,
+        start_mapset,
+        group,
+    )
+    if os.path.isfile(control_points_path):
+        grass.message(_("CONTROL_POINTS file already in group folder."))
+    else:
+        try:
+            shutil.copyfile(grasspoints, control_points_path)
+        except Exception as e:
+            grass.fatal(_("Unable to copy CONTROL_POINTS file."))
+
+    # run i.ortho.position
+    keylist = list(
+        grass.parse_command(
+            "i.ortho.position",
+            group=group,
+            omega=omega,
+            height=default_height,
+            nw=nw,
+            ne=ne,
+            sw=sw,
+            se=se,
+            proj=proj,
+        ).keys()
+    )
+    omega_est = keylist[1].split(": ")[1]
+    phi_est = keylist[2].split(": ")[1]
+    kappa_est = keylist[3].split(": ")[1]
+    xc_est = keylist[4].split(": ")[1]
+    yc_est = keylist[5].split(": ")[1]
+    zc_est = str(float(default_height))
+
+    # run i.ortho.init
+    grass.run_command(
+        "i.ortho.init",
+        group=group,
+        xc=xc_est,
+        yc=yc_est,
+        zc=zc_est,
+        omega=omega_est,
+        phi=phi_est,
+        kappa=kappa_est,
+        xc_sd="1000",
+        yc_sd="1000",
+        zc_sd="1000",
+        omega_sd="0.1",
+        phi_sd="0.1",
+        kappa_sd="0.1",
+        flags="r",
+    )
+
+    # run i.ortho.transform
+    grass.run_command("i.ortho.transform", group=group, flags="sp", verbose=True)
+
+    if not transform_only:
+        if options["logfile"]:
+            # write logfile
+            # i.ortho.transform writes to stderr and stdout, we want to catch both
+            stdout_stderr = grass.Popen(
+                "i.ortho.transform -sp group=%s --v" % group,
+                shell=True,
+                stderr=subprocess.PIPE,
+                stdout=subprocess.PIPE,
+            ).communicate()
+            camerastats = stdout_stderr[1].decode("utf-8")
+            pointstats = stdout_stderr[0].decode("utf-8")
+            logstring = camerastats + pointstats
+            logfile = options["logfile"]
+            if os.path.isfile(logfile):
+                grass.warning(
+                    _("File %s already exists, will be overwritten..." % logfile)
+                )
+            with open(logfile, "w") as file:
+                file.write(logstring)
+        # run i.ortho.rectify
+        grass.run_command(
+            "i.ortho.rectify",
+            group=group,
+            extension=map_extension,
+            resolution=target_resolution,
+            mem=memory,
+            flags="ap",
+            verbose=True,
+        )
+
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    sys.exit(main())
diff --git a/grass_addons/i.ortho.position/Makefile b/grass_addons/i.ortho.position/Makefile
new file mode 100644
index 0000000..c4a9607
--- /dev/null
+++ b/grass_addons/i.ortho.position/Makefile
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = i.ortho.position
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
diff --git a/grass_addons/i.ortho.position/i.ortho.position.html b/grass_addons/i.ortho.position/i.ortho.position.html
new file mode 100644
index 0000000..f4619b1
--- /dev/null
+++ b/grass_addons/i.ortho.position/i.ortho.position.html
@@ -0,0 +1,26 @@
+<h2>DESCRIPTION</h2>
+
+<em>i.ortho.position</em> estimates the camera postion from the corner
+coordinates of the scene in longitudes and latitudes, omega and flight height.
+
+<p>
+The estimated parameters are kappa and the coordinates of the camera
+at the time of exposure. These values can be used for <em>i.ortho.init</em>
+<p>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="i.ortho.photo.html">i.ortho.photo</a>,
+<a href="i.ortho.init.html">i.ortho.init</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Markus Metz, <a href="https://www.mundialis.de/">mundialis</a>
+
+<!--
+<p>
+<i>Last changed: $Date$</i>
+-->
diff --git a/grass_addons/i.ortho.position/i.ortho.position.py b/grass_addons/i.ortho.position/i.ortho.position.py
new file mode 100644
index 0000000..da1be97
--- /dev/null
+++ b/grass_addons/i.ortho.position/i.ortho.position.py
@@ -0,0 +1,241 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+############################################################################
+#
+# MODULE:     i.ortho.position
+# AUTHOR(S):  Markus Metz
+# PURPOSE:    Estimate sensor position from scene center and angles
+#
+# COPYRIGHT:  (C) 2020 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: Estimate sensor position from scene center and angles
+# % keyword: imagery
+# % keyword: orthorectification
+# %end
+# %option G_OPT_I_GROUP
+# % guisection: Input
+# % required: yes
+# %end
+# %option G_OPT_M_COORDS
+# % key: nw
+# % label: North-West corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+# %option G_OPT_M_COORDS
+# % key: ne
+# % label: North-East corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+# %option G_OPT_M_COORDS
+# % key: sw
+# % label: South-West corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+# %option G_OPT_M_COORDS
+# % key: se
+# % label: South-East corner coordinates in ll (lon,lat)
+# % description: must be decimal degrees
+# % required: yes
+# %end
+# %option
+# % key: proj
+# % type: string
+# % label: PROJ definition of the target CRS
+# % description: e.g. +proj=utm +zone=32 +datum=WGS84 +units=m
+# % required: yes
+# %end
+# %option
+# % key: omega
+# % type: double
+# % description: Omega (pitch): Raising or lowering of the aircraft's front (turning around the wings' axis)
+# % required: no
+# %end
+# %option
+# % key: height
+# % type: double
+# % description: Flight height
+# % required: yes
+# %end
+
+# import library
+import os
+import sys
+import threading
+import math
+import re
+import grass.script as gscript
+from grass.script.utils import separator, parse_key_val, encode, decode
+from grass.script import core as gcore
+
+
+class TrThread(threading.Thread):
+    def __init__(self, ifs, inf, outf):
+        threading.Thread.__init__(self)
+        self.ifs = ifs
+        self.inf = inf
+        self.outf = outf
+
+    def run(self):
+        while True:
+            line = self.inf.readline()
+            if not line:
+                break
+            line = line.replace(self.ifs, " ")
+            line = encode(line)
+            self.outf.write(line)
+            self.outf.flush()
+
+        self.outf.close()
+
+
+def main():
+    # check if GISBASE is set
+    if "GISBASE" not in os.environ:
+        # return an error advice
+        print("You must be in GRASS GIS to run this program.")
+        sys.exit(1)
+
+    # input group
+    group = options["group"]
+    # input nw
+    nw = options["nw"]
+    # input ne
+    ne = options["ne"]
+    # input sw
+    sw = options["sw"]
+    # input se
+    se = options["se"]
+    # input proj
+    projstring = options["proj"]
+    # input flight height
+    height = float(options["height"])
+    # input omega
+    omega = 0
+    if options["omega"]:
+        omega = float(options["omega"])
+
+    # check for cct
+    have_cct = True
+    coordsep = " "
+    if not gcore.find_program("cct"):
+        have_cct = False
+        coordsep = "\t"
+        if not gcore.find_program("proj"):
+            gcore.fatal(
+                _(
+                    "Neither cct nor proj program found, install PROJ first: \
+                https://proj.org"
+                )
+            )
+
+    nw_east = float(nw.split(",")[0])
+    nw_north = float(nw.split(",")[1])
+    ne_east = float(ne.split(",")[0])
+    ne_north = float(ne.split(",")[1])
+    sw_east = float(sw.split(",")[0])
+    sw_north = float(sw.split(",")[1])
+    se_east = float(se.split(",")[0])
+    se_north = float(se.split(",")[1])
+
+    # check south
+    if sw_east > se_east:
+        gcore.fatal("se must be east of sw")
+    # check north
+    if nw_east > ne_east:
+        gcore.fatal("ne must be east of nw")
+    # check west
+    if sw_north > nw_north:
+        gcore.fatal("nw must be north of sw")
+    # check east
+    if se_north > ne_north:
+        gcore.fatal("ne must be north of se")
+
+    scene_center_east = (nw_east + ne_east + sw_east + se_east) / 4.0
+    scene_center_north = (nw_north + ne_north + sw_north + se_north) / 4.0
+
+    dxn = ne_east - nw_east
+    dyn = ne_north - nw_north
+    dxs = se_east - sw_east
+    dys = se_north - sw_north
+
+    kappan = math.degrees(math.atan2(dyn, dxn))
+    kappas = math.degrees(math.atan2(dys, dxs))
+    kappa = (kappan + kappas) / 2.0
+
+    # convert ll coordinates to target CRS
+    # shell: echo 9.19625 36.47025 | cct -z0 -t0 +proj=utm +zone=32 +datum=WGS84 +units=m
+    tmpfile = gcore.tempfile()
+    fd = open(tmpfile, "w")
+    fd.write("%s %s\n" % (scene_center_east, scene_center_north))
+    fd.close()
+    inf = open(tmpfile)
+
+    outf = sys.stdout
+
+    if have_cct:
+        cmd = ["cct"] + ["-z0"] + ["-t0"] + projstring.split()
+    else:
+        cmd = ["proj"] + projstring.split()
+
+    p = gcore.Popen(cmd, stdin=gcore.PIPE, stdout=gcore.PIPE)
+
+    ifs = " "
+    tr = TrThread(ifs, inf, p.stdin)
+    tr.start()
+
+    x = y = 0
+
+    for line in p.stdout:
+        try:
+            # incredibly unfriendly output format of cct
+            line = re.sub(" +", " ", decode(line).strip())
+            outcoords = line.split(coordsep)
+            x = outcoords[0]
+            y = outcoords[1]
+        except ValueError:
+            gcore.fatal(line)
+        # only one line
+        break
+
+    scene_center_east = float(x)
+    scene_center_north = float(y)
+
+    p.wait()
+
+    if p.returncode != 0:
+        gcore.warning(_("Projection transform probably failed, please investigate"))
+
+    tan_omega = math.tan(math.radians(omega))
+    ground_offset = height * tan_omega
+    dx = ground_offset * math.sin(math.radians(kappa))
+    dy = -ground_offset * math.cos(math.radians(kappa))
+    camera_x = scene_center_east + dx
+    camera_y = scene_center_north + dy
+
+    outf.write("summary:\n")
+    outf.write("omega: %.3f\n" % omega)
+    outf.write("phi: 0\n")
+    outf.write("kappa: %.3f\n" % kappa)
+    outf.write("camera east: %.2f\n" % camera_x)
+    outf.write("camera north: %.2f\n\n" % camera_y)
+
+    outf.write("Parameters for i.ortho.init:\n")
+    outf.write(
+        "i.ortho.init -r group=%s xc=%.2f yc=%.2f zc=%.2f omega=%.3f phi=0 kappa=%.3f xc_sd=1000 yc_sd=1000 zc_sd=1000 omega_sd=0.1 phi_sd=0.1 kappa_sd=0.1\n"
+        % (group, camera_x, camera_y, height, omega, kappa)
+    )
+
+
+if __name__ == "__main__":
+    options, flags = gscript.parser()
+    sys.exit(main())
diff --git a/scripts/README.md b/scripts/README.md
new file mode 100644
index 0000000..12a8e01
--- /dev/null
+++ b/scripts/README.md
@@ -0,0 +1,6 @@
+## Bash scripts
+- `raster_patch_xy.sh`: patches the pieces together and outputs a raster in the GRASS location/mapset, as well as a GeoTiff and a .vrt with a rough georeferencing. Required in step [../documentation/1_stitching.md](../documentation/1_stitching.md)
+
+- `hugin_automatic_stitching.sh`: automatically stitches the pieces together and outputs a 2-band-raster (1st band: image, 2nd band: alpha channel). Required in step [../documentation/1_automatic_stitching.md](../documentation/1_automatic_stitching.md)
+
+- `gcp_qgis2grass.sh`: the GCPs collected in QGIS have to be adapted to be processable by GRASS, which is done by this script. Required in step [../documentation/3_gcp_collection.md](../documentation/3_gcp_collection.md)
diff --git a/scripts/gcp_qgis2grass.sh b/scripts/gcp_qgis2grass.sh
new file mode 100755
index 0000000..861ebb6
--- /dev/null
+++ b/scripts/gcp_qgis2grass.sh
@@ -0,0 +1,51 @@
+#!/bin/sh
+
+# run inside GRASS
+# run in the target location
+
+qgisgcpfile="$1"
+dem=$2
+utmzone=$3
+
+grassgcpfile="${qgisgcpfile}.grass"
+
+rm -f "$grassgcpfile"
+
+g.region rast=$dem
+
+counter=0
+
+# cct available ?
+HAVE_CCT=`which cct`
+
+# clean comments from input file
+cat ${qgisgcpfile} | grep -v '^#' > ${qgisgcpfile}.clean
+
+
+while read LINE ; do
+  if [ $counter -gt 0 ] ; then
+    mapX=`echo $LINE | cut -d ',' -f1`
+    mapY=`echo $LINE | cut -d ',' -f2`
+    pixelX=`echo $LINE | cut -d ',' -f3`
+    pixelY=`echo $LINE | cut -d ',' -f4`
+
+    # convert to utm32/31n
+    if [ "$HAVE_CCT" ] ; then
+      UTMCOORDS=`echo $mapX $mapY | cct -z0 -t0 +proj=utm +zone=$utmzone +datum=WGS84 +units=m`
+    else
+      UTMCOORDS=`echo $mapX $mapY | proj +proj=utm +zone=$utmzone +datum=WGS84 +units=m`
+    fi
+    mapX=`echo $UTMCOORDS | tr -s ' ' | cut -d ' ' -f1`
+    mapY=`echo $UTMCOORDS | tr -s ' ' | cut -d ' ' -f2`
+  
+    height=`r.what map=$dem coordinates=${mapX},${mapY} | cut -d '|' -f4`
+
+    # for i.ortho.rectify
+    echo "$pixelX $pixelY 0.0    $mapX $mapY $height    1" >>"$grassgcpfile"
+    # for i.rectify
+    #echo "$pixelX $pixelY    $mapX $mapY    1" >>"$grassgcpfile"
+  
+  fi
+  counter=`expr $counter + 1`
+done < "$qgisgcpfile"
+
diff --git a/scripts/hugin_automatic_stitching.sh b/scripts/hugin_automatic_stitching.sh
new file mode 100644
index 0000000..72101bc
--- /dev/null
+++ b/scripts/hugin_automatic_stitching.sh
@@ -0,0 +1,76 @@
+#!/bin/sh
+
+## Hugin based stitching
+
+# fail immediately on error
+set -e
+
+# Requirements: 30++ GB RAM for CORONA scenes, 65++ GB RAM for GAMBIT scenes, 55++ GB RAM for HEXAGON scenes
+# Ubuntu: apt-get install hugin
+# Fedora: dnf install hugin
+
+# NOTE: only a single scene per directory!
+
+###
+# variables to adjust
+MISSION=$1
+
+if [ "$MISSION" == "GAMBIT" ]; then
+    BASENAME=$(basename DZ*_a.tif _a.tif)
+elif [ "$MISSION" == "CORONA" ]; then
+    BASENAME=$(basename DS*_a.tif _a.tif)
+elif [ "$MISSION" == "HEXAGON" ]; then
+    BASENAME=$(basename D3*_a.tif _a.tif)
+else
+    echo "Automatic stitching does not support this satellite mission."
+fi
+
+# generate Hugin project file
+# Note: if there are more less scene pieces the following part has to be adjusted
+if [ "$MISSION" == "GAMBIT" ]; then
+    pto_gen -f 1 ${BASENAME}_b.tif ${BASENAME}_a.tif
+elif [ "$MISSION" == "CORONA" ]; then
+    pto_gen ${BASENAME}_d.tif ${BASENAME}_c.tif ${BASENAME}_b.tif ${BASENAME}_a.tif
+elif [ "$MISSION" == "HEXAGON" ]; then
+    pto_gen ${BASENAME}_f.tif ${BASENAME}_e.tif ${BASENAME}_d.tif ${BASENAME}_c.tif ${BASENAME}_b.tif ${BASENAME}_a.tif
+fi
+
+# search for GCPs
+METHOD=linearmatch
+if [ "$MISSION" == "GAMBIT" ]; then
+    cpfind -o output_${BASENAME}.pto ${BASENAME}_b-${BASENAME}_a.pto
+elif [ "$MISSION" == "CORONA" ]; then
+    cpfind --$METHOD -o output_${BASENAME}.pto ${BASENAME}_d-${BASENAME}_a.pto
+elif [ "$MISSION" == "HEXAGON" ]; then
+    cpfind --$METHOD -o output_${BASENAME}.pto ${BASENAME}_f-${BASENAME}_a.pto
+fi
+
+# pruning GCPs
+cpclean -o project_${BASENAME}.pto output_${BASENAME}.pto
+
+# Optimising positions and geometry
+autooptimiser -a -l -s -m -o project_${BASENAME}.pto project_${BASENAME}.pto
+
+# modify panorama output parameters
+if [ "$MISSION" == "GAMBIT" ]; then
+pano_modify --output modify_${BASENAME}.pto --canvas=AUTO project_${BASENAME}.pto
+elif [ "$MISSION" == "CORONA" ]; then
+pano_modify --output modify_${BASENAME}.pto --canvas=100% --crop=AUTO project_${BASENAME}.pto
+elif [ "$MISSION" == "HEXAGON" ]; then
+pano_modify --output modify_${BASENAME}.pto --canvas=100% --crop=AUTO project_${BASENAME}.pto
+fi
+
+# stitching to a single TIFF
+hugin_executor --stitching --prefix=stitched_${BASENAME} modify_${BASENAME}.pto
+echo "Written <stitched_${BASENAME}.tif>"
+
+echo "If needed, rotate scene by 180 degree, using:
+convert stitched_${BASENAME}.tif  -rotate 180 -compress lzw stitched_${BASENAME}_flipped.tif
+"
+
+####### Alternative in GRASS GIS:
+
+# i.points.auto
+# https://grass.osgeo.org/grass7/manuals/addons/i.points.auto.html
+
+# see Abbildung 4 in https://www.grassbook.org/wp-content/uploads/neteler/papers/neteler2005_IJG_051-061_draft.pdf
diff --git a/scripts/raster_patch_xy.sh b/scripts/raster_patch_xy.sh
new file mode 100755
index 0000000..8325133
--- /dev/null
+++ b/scripts/raster_patch_xy.sh
@@ -0,0 +1,99 @@
+#!/bin/bash
+
+# patch different parts of a CORONA scene together
+# the order of patching must match the order of alignment
+# i.e. if map a is in the display on top of map b
+# the order must be map a, map b
+
+# variables to adjust
+
+SCENEID=$1
+NWLON=$2
+NWLAT=$3
+NELON=$4
+NELAT=$5
+SWLON=$6
+SWLAT=$7
+SELON=$8
+SELAT=$9
+# odered list assuming a on top of b (on top of c on top of d)
+# 4 image parts
+ORDEREDLIST="${SCENEID}_a,${SCENEID}_b,${SCENEID}_c,${SCENEID}_d"
+# 2 image parts
+#ORDEREDLIST="${SCENEID}_a,${SCENEID}_b"
+
+# set FLIP to non-zero if the image needs to be
+# flipped top-bottom and left-right
+# applies to KH-4A aft
+FLIP=0
+
+# nothing to change below
+export GRASS_OVERWRITE=1
+
+g.region -pa res=1 rast=$ORDEREDLIST
+r.patch in=$ORDEREDLIST out=${SCENEID}_patch
+
+south=0
+west=0
+eval `r.info -g ${SCENEID}_patch`
+
+if [ "$south" != "0" ] ; then
+  # invert sign
+  if [ "${south:0:1}" = "-" ] ; then
+    south="+${south:1}"
+  else
+    south="-${south}"
+  fi
+  r.region map=${SCENEID}_patch n=n$south s=s$south
+fi
+
+if [ "$west" != "0" ] ; then
+  # invert sign
+  if [ "${west:0:1}" = "-" ] ; then
+    west="+${west:1}"
+  else
+    west="-${west}"
+  fi
+  r.region map=${SCENEID}_patch w=w$west e=e$west
+fi
+
+g.region rast=${SCENEID}_patch
+r.mapcalc "\"${SCENEID}_zero\" = if(isnull(\"${SCENEID}_patch\"), 0, \"${SCENEID}_patch\")"
+g.remove rast name=${SCENEID}_patch -f
+
+if [ $FLIP -ne 0 ] ; then
+  r.flip -b in=${SCENEID}_zero out=${SCENEID}
+  g.remove rast name=${SCENEID}_zero -f
+else
+  g.rename rast=${SCENEID}_zero,${SCENEID}
+fi
+
+r.out.gdal -cm overviews=4 createopt="COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES,BIGTIFF=YES" in=${SCENEID} out=${SCENEID}.tif
+
+# get rows and cols from r.info
+eval $(r.info -g map=$SCENEID)
+
+# GDAL counts rows from top to bottom:
+# upper left (ul): 0 0
+# lower left (ll): 0 nrows
+# upper right (ur): ncols 0
+# lower right (lr): ncols nrows
+
+# KH-4A:
+# GDAL ul = NW
+# GDAL ll = SW
+# GDAL ur = NE
+# GDAL lr = SE
+
+# KH-7:
+# GDAL ul = SW
+# GDAL ll = SE
+# GDAL ur = NW
+# GDAL lr = NE
+
+# run gdal_translate
+gdal_translate -of VRT -a_srs EPSG:4326 \
+               -gcp 0 0 $NWLON $NWLAT \
+               -gcp 0 $rows $SWLON $SWLAT \
+               -gcp $cols 0 $NELON $NELAT \
+               -gcp $cols $rows $SELON $SELAT ${SCENEID}.tif ${SCENEID}.vrt