Some Code Summaries Regarding GDAL Usage

1. Get TIF Extent

fun getExtent(path: String) {
            gdal.AllRegister()
            val driver = gdal.GetDriverByName("GTiff")

            val dataset = gdal.Open(path)
            val x = dataset.rasterXSize
            val y = dataset.rasterYSize
            val transforme = dataset.GetGeoTransform()


            val west = transforme[0]
            val dx = transforme[1]
            val north = transforme[3]
            val dy = transforme[5]

            val east = west + dx * x
            val sourth = north + dy * y

            println("${west},${east},${north},${sourth} ")

//            val ds = ogr.Open(path)
//            println(layer.GetExtent())


        }

2. Merge Multiple TIFs

      fun union(path0: String, name: String): String {
            gdal.AllRegister()

            val d0 = gdal.Open(path0)
            val d1 = gdal.Open("${base}img_cr/${name}")

            gdal.Warp("${base}union/${name}", listOf<Dataset>(d0, d1).toTypedArray(), null)

            // Resolve resource occupation issue
            gdal.GDALDestroyDriverManager()
            return "${base}union/${name}"

        }

3. Remove TIF Black Borders (TIF only)

    fun nearBlack(photo: Photo) {
          gdal.AllRegister()
          val data = gdal.Open("${base}tif\\${photo.File}.tif")
          val vector = Vector<String>()
          "-of GTiff -setalpha".split(" ").forEach { vector.add(it) }
          gdal.Nearblack("${base}tif\\${photo.File}.tif", data, NearblackOptions(vector))
      }

4. Convert to WKT in a Specific Coordinate System

     val proj4_ = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs"

     private fun toWKt(feature: Feature): String {
            val geom = feature.GetGeometryRef()
            val proj4 = geom.GetSpatialReference().ExportToProj4()
            if (!proj4.equals(proj4_)) {
                val sp = SpatialReference()
                sp.ImportFromProj4(proj4_)
                geom.TransformTo(sp)
            }
            geom.FlattenTo2D()
            return geom.ExportToWkt()
        }

5. Convert TIF to a 3-Band RGB TIF (Can resolve color distortion during merging)

fun handleListTif2Rgb(z: Int) {

    val vector = Vector<String>()
    vector.addElement("-expand")
    vector.addElement("RGB")

    val parent = "E:\\img2tif\\${z}_tif\\"
    val parent2 = "E:\\img2tif\\${z}_tif_rgb\\"
    for (it in File(parent).list()!!) {
        println("Start $it")
        val yyy = File(parent + it).list()!!
        val datas = mutableListOf<Dataset>()
        if (!File("${parent2}\\${it}").exists()) File("${parent2}\\${it}").mkdir()
        for (i in yyy.indices) {
            val it2 = yyy[i]
            val outPath = "${parent2}\\${it}\\${it2}"
            if (!File(outPath).exists()) {
                val dsrc = gdal.Open("${parent}${it}\\${it2}")
                gdal.Translate(outPath, dsrc, TranslateOptions(vector))
            }
            datas.add(gdal.Open(outPath))
        }
        gdal.Warp("E:\\img2tif\\14_union\\$it.tif", datas.toTypedArray(), null)
    }
}

6. Modify TIF Resolution via Resampling

• Using the warp method

  Reference: https://www.jianshu.com/p/821c741ff169

gdal.Warp(tif_out, tif_input, resampleAlg=gdalconst.GRA_NearestNeighbour, xRes=500, yRes=500)

• Using GeoTransform

  Citation: https://blog.csdn.net/gisuuser/article/details/106304155

from osgeo import gdal, gdalconst
import os
import numpy as np


def resampling(source_file, target_file, scale=1.0):
    """
    Image resampling
    :param source_file: Source file
    :param target_file: Output image
    :param scale: Pixel scaling ratio
    :return:
    """
    dataset = gdal.Open(source_file, gdalconst.GA_ReadOnly)
    band_count = dataset.RasterCount  # Number of bands

    if band_count == 0 or not scale > 0:
        print("Parameter exception")
        return

    cols = dataset.RasterXSize  # Number of columns
    rows = dataset.RasterYSize  # Number of rows
    cols = int(cols * scale)  # Calculate new number of columns
    rows = int(rows * scale)

    geotrans = list(dataset.GetGeoTransform())
    print(dataset.GetGeoTransform())
    print(geotrans)
    geotrans[1] = geotrans[1] / scale  # Pixel width becomes scale times the original
    geotrans[5] = geotrans[5] / scale  # Pixel height becomes scale times the original
    print(geotrans)

    if os.path.exists(target_file) and os.path.isfile(target_file):  # If an image with the same name exists
        os.remove(target_file)  # Then delete it

    band1 = dataset.GetRasterBand(1)
    data_type = band1.DataType
    target = dataset.GetDriver().Create(target_file, xsize=cols, ysize=rows, bands=band_count,
                                        eType=data_type)
    target.SetProjection(dataset.GetProjection())  # Set projection coordinates
    target.SetGeoTransform(geotrans)  # Set geographic transform parameters
    total = band_count + 1
    for index in range(1, total):
        # Read band data
        print("Writing band " + str(index))
        data = dataset.GetRasterBand(index).ReadAsArray(buf_xsize=cols, buf_ysize=rows)
        out_band = target.GetRasterBand(index)
        out_band.SetNoDataValue(dataset.GetRasterBand(index).GetNoDataValue())
        out_band.WriteArray(data)  # Write data to the new image
        out_band.FlushCache()
        out_band.ComputeBandStats(False)  # Compute statistics
    print("Writing complete")
    del dataset
    del target


if __name__ == "__main__":
    source_file = r"E:\商丘yx\相交4.tiff"
    target_file = r"E:\商丘yx\相交411.tiff"
    resampling(source_file, target_file, scale=1.1)
    target_file = r"E:\商丘yx\相交05.tiff"
    resampling(source_file, target_file, scale=0.5)

AI Translation | AI 翻译

This article was translated from Chinese to English by AI. If there are any inaccuracies, please refer to the original Chinese version.


本文由 AI 辅助从中文翻译为英文。如遇不准确之处，请以中文原版为准。