#383 Support raster option on gdal2tiles_3.5.2/gdal2customtiles (#397)

* #383 gdal2tiles_3.5.2_v2 for raster support

* #383 Cleanup gdal2tiles scripts and improve documentation

* #383 gdal2customtiles extentworld working, overview tiles not working

* #383 Fix gdal2customtiles rasters with different pixel scales

* gdal2customtiles - fix width calc, still slight offset

* gdal2customtiles raster - fix rounding issues

auxiliary/gdal2customtiles/legacy/readme.md

@@ -0,0 +1,68 @@
+# gdal2customtiles
+
+This wraps together:
+
+- gdal2tiles4extent.py
+- gdal2tiles1bto4b_v3.py
+
+---
+
+## Raster Extents:
+
+Tile partial world rasters.
+
+**Requires:**
+
+- `-p raster`
+  - This is necessary for generating tiles with a custom extent.
+- `-x` OR `--extentworld` followed by values `ulx,uly,lrx,lry,pixel_resolution`
+  - The extentworld is the full bounding area of the projection for the planetary body. The extentworld is the full bounding area of the projection for the planetary body. Units are in meters using upper left (ul) and lower right (lr) order. These values are reported from gdalinfo. Units are in meters using upper left (ul) and lower right (lr) order. These values are reported from gdalinfo. Values are separated by commas with no spaces.
+
+**Example:**
+
+```
+python gdal2customtiles.py -p raster --extentworld -4022404.001,4022036.893,-4022036.893,4022404.001,367.108150109358121 input.tif output_dir
+```
+
+_Notes:_
+
+- Only works if set zoom (-z 0-10) encompasses the native zoom of the raster.
+- 'ERROR 5's are expected.
+
+---
+
+## Digital Elevation Model Tiles:
+
+Generate Digital Elevation Maps (DEMs) tiles.
+
+Any 32-bit image data can be encoded into the RGBA channels of a PNG. MMGIS uses this file type to create terrain meshes as well as for a data layer.
+
+- Certain resampling methods can corrupt `--dem` results.
+
+**Requires:**
+
+- `-m` or `--dem`
+
+**Example:**
+
+```
+python gdal2customtiles.py -p raster --extentworld -4022404.001,4022036.893,-4022036.893,4022404.001,367.108150109358121 --dem inputdem.tif output_dir
+```
+
+_Notes:_
+
+- Does not include the convenience of rasterstotiles.py yet.
+- Can only tile 32-bit images with --dem option.
+
+## gdal2tiles_3.5.2.py
+
+- `rasters2customtiles_3.5.2.py` and `gdal2tiles_3.5.2.py` support only the `--dem` option (and not `--raster` yet). `-m` no longer works and must be `--dem`. Tested with gdal 3.4.3. Upgraded to support multi-processes. See `python rasters2customtiles_3.5.2.py --help`. Unlike `gda2customtiles.py`, does not seam-match DEM tiles (better for Data Layers and Viewshed Tool, bad for 3D Globe).
+- Adds the resampling algorithm `near-composite` that uses nearest-neighbor and ovarlays the new tile onto the old tile (if any in output directory)
+- Certain resampling methods can corrupt `--dem` results.
+- To support the value 0, all 0 data values get mapped to to the value 2^31 (2147483648) (RGBA=79,0,0,0) and then decoded by the MMGIS reader back to 0. This avoids clashes with other nondata-like values writing to 0,0,0,0 in the outputted pngs.
+
+**Example:**
+
+```
+python gdal2tiles_3.5.2.py --dem input.tif output_dir --srcnodata=-9999 -r near-composite --tilesize=128
+```

auxiliary/gdal2customtiles/rasters2customtiles.py

@@ -0,0 +1,151 @@
+import sys
+import subprocess
+import optparse
+from osgeo import gdal, osr
+
+
+def optparse_init() -> optparse.OptionParser:
+    """Prepare the option parser for input (argv)"""
+
+    usage = "Usage: %prog [options] input_file [output]"
+    p = optparse.OptionParser(usage)
+    p.add_option(
+        "--dem",
+        action="store_true",
+        dest="isDEMtile",
+        help="Indicate if the input is a Digital Elevation Model"
+    )
+    p.add_option(
+        "--processes",
+        dest="processes",
+        type="int",
+        help="Number of processes to use for tiling",
+    )
+    p.add_option(
+        "--tilesize",
+        dest="tilesize",
+        metavar="PIXELS",
+        type="int",
+        help="Width and height in pixel of a tile. Defaults to 256 (or 32 for --dem)",
+    )
+    p.add_option(
+        "-z",
+        "--zoom",
+        dest="zoom",
+        help="Zoom levels to render (format:'2-5', '10-' or '10').",
+    )
+    p.add_option(
+        "-e",
+        "--resume",
+        dest="resume",
+        action="store_true",
+        help="Resume mode. Generate only missing files.",
+    )
+    p.add_option(
+        "-a",
+        "--srcnodata",
+        dest="srcnodata",
+        metavar="NODATA",
+        help="Value in the input dataset considered as transparent",
+    )
+    return p
+
+
+def GetExtent(gt, cols, rows):
+    ''' Return list of corner coordinates from a geotransform
+
+        @type gt:   C{tuple/list}
+        @param gt: geotransform
+        @type cols:   C{int}
+        @param cols: number of columns in the dataset
+        @type rows:   C{int}
+        @param rows: number of rows in the dataset
+        @rtype:    C{[float,...,float]}
+        @return:   coordinates of each corner
+    '''
+    ext = []
+    xarr = [0, cols]
+    yarr = [0, rows]
+
+    for px in xarr:
+        for py in yarr:
+            x = gt[0]+(px*gt[1])+(py*gt[2])
+            y = gt[3]+(px*gt[4])+(py*gt[5])
+            ext.append([x, y])
+        yarr.reverse()
+    return ext
+
+
+def ReprojectCoords(coords, src_srs, tgt_srs):
+    ''' Reproject a list of x,y coordinates.
+
+        @type geom:     C{tuple/list}
+        @param geom:    List of [[x,y],...[x,y]] coordinates
+        @type src_srs:  C{osr.SpatialReference}
+        @param src_srs: OSR SpatialReference object
+        @type tgt_srs:  C{osr.SpatialReference}
+        @param tgt_srs: OSR SpatialReference object
+        @rtype:         C{tuple/list}
+        @return:        List of transformed [[x,y],...[x,y]] coordinates
+    '''
+    trans_coords = []
+    transform = osr.CoordinateTransformation(src_srs, tgt_srs)
+    for x, y in coords:
+        x, y, z = transform.TransformPoint(x, y)
+        trans_coords.append([x, y])
+    return trans_coords
+
+
+def AutoGdalTranslate(geo_extent, cols, rows, raster):
+    gdal_translate = "gdal_translate -of VRT -a_srs EPSG:4326 -gcp 0 0 " + str(geo_extent[0][0]) + " " + str(geo_extent[0][1]) + " -gcp " + str(cols) + " 0 " + str(geo_extent[3][0]) + " " + str(
+        geo_extent[3][1]) + " -gcp " + str(cols) + " " + str(rows) + " " + str(geo_extent[2][0]) + " " + str(geo_extent[2][1]) + " " + raster + " " + raster[:-4] + ".vrt"
+    print(f"Running: {gdal_translate}\n")
+    subprocess.Popen(gdal_translate)
+
+
+def AutoGdal2Tiles(raster, options, outputdir):
+    dem = ""
+    if options.isDEMtile is True:
+        dem = " --dem"
+    processes = ""
+    if options.processes is not None:
+        processes = f" --processes={options.processes}"
+    tilesize = ""
+    if options.tilesize is not None:
+        tilesize = f" --tilesize={options.tilesize}"
+    zoom = ""
+    if options.zoom is not None:
+        zoom = f" --zoom={options.zoom}"
+    resume = ""
+    if options.resume is True:
+        resume = " --resume"
+    srcnodata = " --srcnodata=0,0,0"
+    if options.srcnodata is not None:
+        srcnodata = f" --srcnodata={options.srcnodata}"
+    output = ""
+    if outputdir is not None:
+        output = f" {outputdir}"
+    gdal2tiles = f"python gdal2customtiles.py -n{dem}{processes}{tilesize}{zoom}{resume}{srcnodata} {raster[:-4]}.vrt{output}"
+    print(f"Running: {gdal2tiles}\n")
+    subprocess.Popen(gdal2tiles)
+
+
+parser = optparse_init()
+options, args = parser.parse_args(args=sys.argv)
+
+raster = args[1]
+ds = gdal.Open(raster)
+
+gt = ds.GetGeoTransform()
+cols = ds.RasterXSize
+rows = ds.RasterYSize
+extent = GetExtent(gt, cols, rows)
+
+src_srs = osr.SpatialReference()
+src_srs.ImportFromWkt(ds.GetProjection())
+tgt_srs = src_srs.CloneGeogCS()
+
+geo_extent = ReprojectCoords(extent, src_srs, tgt_srs)
+
+AutoGdalTranslate(geo_extent, cols, rows, raster)
+AutoGdal2Tiles(raster, options, args[2])

auxiliary/gdal2customtiles/readme.md

@@ -1,43 +1,46 @@
-# gdal2customtiles
+# gdal2customtiles.py
 
-This wraps together:
+_Python 3.10.5_
 
-- gdal2tiles4extent.py
-- gdal2tiles1bto4b_v3.py
+Accepts all [gdal2tiles.py](https://gdal.org/programs/gdal2tiles.html) options. Built off of GDAL 3.5.2 and tested with GDAL 3.4.3, it adds the following new features and capabilities.
 
 ---
 
 ## Raster Extents:
 
-Tile partial world rasters.
+Tile partial world rasters. Useful for tiling non-mercator and non-geodetic projected data.
 
 **Requires:**
 
 - `-p raster`
   - This is necessary for generating tiles with a custom extent.
-- `-x` OR `--extentworld` followed by values `ulx,uly,lrx,lry,pixel_resolution`
+- `--extentworld` followed by values `ulx,uly,lrx,lry,pixel_resolution`
   - The extentworld is the full bounding area of the projection for the planetary body. The extentworld is the full bounding area of the projection for the planetary body. Units are in meters using upper left (ul) and lower right (lr) order. These values are reported from gdalinfo. Units are in meters using upper left (ul) and lower right (lr) order. These values are reported from gdalinfo. Values are separated by commas with no spaces.
 
-**Example:**
+#### Example:
 
 ```
-python gdal2customtiles.py -p raster --extentworld -4022404.001,4022036.893,-4022036.893,4022404.001,367.108150109358121 input.tif output_dir
+python gdal2customtiles.py -p raster --extentworld -931100.000,931100.000,931100.000,-931100.000,100 inputs/WAC_GLOBAL_P900S0000_100M.tif outputs/WAC_GLOBAL_P900S0000_100M
+python gdal2customtiles.py -p raster --extentworld -931100.000,931100.000,931100.000,-931100.000,100 inputs/ldem_87s_5mpp_hillshade.tif outputs/ldem_87s_5mpp_hillshade
 ```
 
+- `WAC_GLOBAL_P900S0000_100M.tif` is in Lunar South Polar projection (IAU2000:30120). Its data covers the full bounds of that projection's world-space (it's world extent/"extentworld") thus we use its bounds and pixel resolution directly from its metadata: `--extentworld -931100.000,931100.000,931100.000,-931100.000,100`
+
+  - _Note: If your basemap does not cover the full world-space, you would need to compute the world-space's bounds and its resolution relative to your datasets_
+
+- `ldem_87s_5mpp_hillshade.tif` is also in Lunar South Polar projection (IAU2000:30120). Its data only covers a small region of the projection's world-space. We still use the previous `--extentworld -931100.000,931100.000,931100.000,-931100.000,100`
+
 _Notes:_
 
 - Only works if set zoom (-z 0-10) encompasses the native zoom of the raster.
-- 'ERROR 5's are expected.
 
 ---
 
 ## Digital Elevation Model Tiles:
 
 Generate Digital Elevation Maps (DEMs) tiles.
 
-Any 32-bit image data can be encoded into the RGBA channels of a PNG. MMGIS uses this file type to create terrain meshes as well as for a data layer.
-
-- Certain resampling methods can corrupt `--dem` results.
+Any 32-bit image data can be encoded into the RGBA channels of a PNG. MMGIS uses this file type to create terrain meshes as well as for Data Layers.
 
 **Requires:**
 
@@ -51,18 +54,33 @@ python gdal2customtiles.py -p raster --extentworld -4022404.001,4022036.893,-402
 
 _Notes:_
 
-- Does not include the convenience of rasterstotiles.py yet.
 - Can only tile 32-bit images with --dem option.
-
-## gdal2tiles_3.5.2.py
-
-- `rasters2customtiles_3.5.2.py` and `gdal2tiles_3.5.2.py` support only the `--dem` option (and not `--raster` yet). `-m` no longer works and must be `--dem`. Tested with gdal 3.4.3. Upgraded to support multi-processes. See `python rasters2customtiles_3.5.2.py --help`. Unlike `gda2customtiles.py`, does not seam-match DEM tiles (better for Data Layers and Viewshed Tool, bad for 3D Globe).
-- Adds the resampling algorithm `near-composite` that uses nearest-neighbor and ovarlays the new tile onto the old tile (if any in output directory)
+- Current `--dem` tiles do not seam-match tile edges. This may or may not be desired (not seam-matching is better for Data Layers and the Viewshed Tool, but bad for MMGIS' 3D Globe/LithoSphere). If seam-matching is desired use `legacy/gdal2customtiles.py` or `legacy/gdal2customtiles_py27.py`
 - Certain resampling methods can corrupt `--dem` results.
 - To support the value 0, all 0 data values get mapped to to the value 2^31 (2147483648) (RGBA=79,0,0,0) and then decoded by the MMGIS reader back to 0. This avoids clashes with other nondata-like values writing to 0,0,0,0 in the outputted pngs.
 
+---
+
+## Compositing Tiles:
+
+Adds the resampling algorithm `near-composite` that uses nearest-neighbor resampling and overlays the new tile onto the old tile (if any in output directory). This makes it possible to accumulate or combine tilesets at the indivdual tile image level. Data in tiles can be overwritten by this process so be cognizant of run order and input extents.
+
 **Example:**
 
 ```
-python gdal2tiles_3.5.2.py --dem input.tif output_dir --srcnodata=-9999 -r near-composite --tilesize=128
+python gdal2customtiles.py -r near-composite --srcnodata=-9999 --processes=40 --tilesize=128 --dem input_A.tif output_dir
+python gdal2customtiles.py -r near-composite --srcnodata=-9999 --processes=40 --tilesize=128 --dem input_B.tif output_dir
 ```
+
+_Notes:_
+
+- Nodata values are treated as transparent and will not overwrite existing pixels in the output tile images.
+
+---
+
+# raster2customtiles.py
+
+A convience script that wraps gda2customtiles.py. Translates the input data into EPSG:4326 and sets proper ground control points. Might be outdated. Use gdal2customtiles directly for the most control.
+
+**Usage:**
+`rasters2customtiles.py [options] input_file [output]` or see `--help`