Importing and exporting GeoTiff images with Python

• (2019/07) Updated some links *

Organize your notes. The code assumes an interactive shell.

Read as numpy array of pixel array data of image

read_tif.py

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

src = gdal.Open('inputfilepath', gdalconst.GA_ReadOnly) #loading tif(read only)
type(src) # "osgeo.gdal.Dataset"

src.RasterXSize #Number of horizontal pixels
src.RasterYSize #Number of pixels in the vertical direction
src.RasterCount #Number of bands

b1 = src.GetRasterBand(1).ReadAsArray() #1st band numpy array
b2 = src.GetRasterBand(2).ReadAsArray() #2nd band numpy array
b3 = src.GetRasterBand(3).ReadAsArray() #3rd band numpy array

dtid = src.GetRasterBand(1).DataType #Model number(ex: 6 -> numpy.float32)
gdal_array.GDALTypeCodeToNumericTypeCode(dtid) #Model number->Type name conversion

# others...

src.GetRasterBand(1).GetNoDataValue()
src.GetRasterBand(1).GetMinimum()
src.GetRasterBand(1).GetMaximum()
src.GetRasterBand(1).GetScale()
src.GetRasterBand(1).GetUnitType()

For the last 5 items, nothing was returned from the file I had. It may be referring to the value stored in the address in the form of metadata in advance. Of course, there are many cases where numpy processing after ReadAsArray () can be managed.

Acquisition of basic information

read_info.py

src.GetDescription() #file name

src.GetGeoTransform() #6 numbers on coordinates(See below)
src.GetProjection() #Coordinate system information

The meaning of the sequence output by GetGeoTransform () is [Start point longitude, West East resolution, Rotation (0 for north-south), Start point end latitude, Rotation (0 for north-south), North-south resolution (negative in the north-south direction)].

Get image metadata

See below, for example, to search for metadata concepts and items in Tiff. For the main metadata, get the metadata using the tag name written here. TIFF format (1) | JProgramer Overview of CG files, Chapter 5, Section 1 GTiff -- GeoTIFF File Format

read_meta.py

src.GetMetadata() #Dictionary-type metadata array
src.GetMetadataItem('itemname') #Item specification(ex: 'TIFFTAG_XRESOLUTION', 'TIFFTAG_DATETIME', ..)

GeoTiff image output

Follow GTiff --GeoTIFF File Format. In particular, please note that the correct operation may not be performed when using the image afterwards unless the settings are made properly, including those specified as options. For example, when processing the output image using some software, that software may refer to some of the metadata (tags) as basic information. If the information is missing, it may not work as expected, and the cause must be identified each time, except for a thorough understanding of the software in advance. (The reason why RGB is input to the Photometric Interpretation tag in the code below is that I expected it to be drawn as a color image in the software used, but there is a problem that it becomes a monochrome display using the band value of the first layer. Because there was.)

write_data.py

from osgeo import osr #Spatial reference module

dtype = gdal.GDT_Float32 #others: gdal.GDT_Byte, ...
band = 3 #Number of bands
output = gdal.GetDriverByName('GTiff').Create('outputfilepath', xsize, ysize, band, dtype, options = ['PHOTOMETRIC=RGB']) #Empty output file

output.SetGeoTransform((ul_x, h_res, 0, ul_y, 0, -v_res)) #Coordinate system specification
srs = osr.SpatialReference() #Spatial reference information
srs.ImportFromEPSG(32648) # WGS84 UTM_Specify the coordinate system for 48n
output.SetProjection(srs.ExportToWkt()) #Combine spatial information

output.GetRasterBand(1).WriteArray(b1)   #Export red band (b1 is numpy 2D array)
output.GetRasterBand(2).WriteArray(b2)   #Green band
output.GetRasterBand(3).WriteArray(b3)   #Blue band
output.FlushCache()                     #Write to disk
output = None