CODEWITHSUNDEEP
pythongisgdalqgisgeotiff
Nathan Hui
Nathan Hui

Reputation: 26

How to ensure data is written to geotiff using GDAL?

I am working on a script that takes in a CSV file containing locations of measurements and creating a probability map based on those measurements. I am using GDAL to write the probabilities to a GeoTIFF, however, I am noticing that none of the data is properly written to the GeoTIFF. When generating the heatmap, I use numpy to store the data in an X by Y array, then write that array to the GeoTIFF. The current dataset I am testing on measure 260 x 262 meters, and the resulting GeoTIFF has the correct dimensions and is properly georeferenced. The heatmap data should have pixel values ranging from -51 to 0, however, when I load the GeoTIFF into QGIS, QGIS only has values ranging from -51 to -31.

The GeoTIFF should be stored as 32 bit float data and lossless compression, so it's not a data range or compression issue.

I've attached the relevant code below:

import sys
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plot
# USING utm 0.4.0 from https://pypi.python.org/pypi/utm
import utm
import os
import argparse
import fileinput
from osgeo import gdal
import osr
import math

# finalLon and finalLat are 1-D arrays containing the UTM coordinates for each measurement
tiffXSize = int(max(finalLon)) - int(min(finalLon)) + maxRange * 2
tiffYSize = int(max(finalLat)) - int(min(finalLat)) + maxRange * 2
pixelSize = 1
heatMapArea = np.zeros((tiffXSize, tiffYSize)) # y, x
refLat = min(finalLat) - maxRange
maxLat = max(finalLat) + maxRange
refLon = max(finalLon) + maxRange
minLon = min(finalLon) - maxRange
# some code to set values in the heatMapArea array

outputFileName = '%s/RUN_06d_COL_%06d.tiff' % (output_path, run_num, num_col)
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(
    outputFileName,
    tiffYSize,
    tiffXSize,
    1,
    gdal.GDT_Float32, ['COMPRESS=LZW'])

spatialReference = osr.SpatialReference()
spatialReference.SetUTM(zonenum, zone >= 'N')
spatialReference.SetWellKnownGeogCS('WGS84')
wkt = spatialReference.ExportToWkt()
retval = dataset.SetProjection(wkt)
dataset.SetGeoTransform((
    refLat,    # 0
    1,  # 1
    0,                      # 2
    refLon,    # 3
    0,                      # 4
    -1))
band = dataset.GetRasterBand(1)
band.SetNoDataValue(100)
print(tiffXSize)
print(tiffYSize)
print(np.amin(heatMapArea))
print(np.amax(heatMapArea))
print(np.mean(heatMapArea))
print(np.std(heatMapArea))
print((heatMapArea > -30).sum())
band.WriteArray(heatMapArea)
band.SetStatistics(np.amin(heatMapArea), np.amax(heatMapArea), np.mean(heatMapArea), np.std(heatMapArea))
dataset.FlushCache()
dataset = None

Upvotes: 0

Views: 1010

Answers (1)

Nathan Hui
Nathan Hui

Reputation: 26

Problem is in how QGIS computes the maximum and minimum of the GeoTIFF. If you go into the properties of the GeoTIFF layer, and open the Style tab, under Band Rendering, you can choose how to load the min/max values. By default, QGIS uses the 2% to 98% count as the minimum and maximum, however, for this, we want the min/max.

Upvotes: 1

Related Questions