How to write a georeferenced GeoTIFF from known coordinates with Python rasterio

Question

I have three numpy arrays, one with data, and two more with explicit longitudes and latitudes. I'd like to write this data to make a GeoTIFF file with rasterio such that I can load it with QGIS. The data comes originally from a known lambert conformal projection, if it helps. How can I save this matrix in GeoTIFF format?

Trying to be more precise, I'm following this tutorial to get radar data. Following it, I have two matrices x and y of shape (720, 1832), and an additional matrix ref_datawith the same shape, which contains the actual data. The matrices x and y represent the distance in meters in the zonal and meridional distances to a given point of coordinates (-85.94388888888889, 37.97527777777778). Now, I'd like to diverge from the tutorial in that I do not want to get a plot using cartopy and matplotlib, but to save a GeoTIFF image with rasterio. For this I need the affine transformation (I believe). How can I get it and save this data as a raster properly georeferenced?

Is it possible to do this with Python and rasterio?

Answer 1

Assuming you have another raster layer with the same CRS and dimensions, you can use a function like this to convert your numpy array to a georeferenced raster using the extent and geotransform from the input layer:

import numpy as np
import gdal, osr 

def write_geotiff(array, gdal_obj, outputpath, dtype=gdal.GDT_UInt16, options=0, color_table=0, nbands=3, nodata=False):
    """
    Description:
    Writes a geotiff from a Numpy array with appended georeferencing from parent geotiff.
    
    Parameters:
    array: numpy array to write as geotiff
    gdal_obj: object created by gdal.Open() using a tiff that has the SAME CRS, transformation, and resolution as the array you're writing
    outputpath: path including filename.tiff
    dtype (OPTIONAL): datatype to save as
    nodata (default: False): set to any value you want to use for nodata; if False, nodata is not set
    """

    gt = gdal_obj.GetGeoTransform()

    width = np.shape(array)[1]
    height = np.shape(array)[0]

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(array)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    wkt = gdal_obj.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())

    # Close output raster dataset 
    dest = None

Answer 2

You can do that with rioxarray: First open the raster, second transform to geographic coordinates and finally transform to UTM coordiantes.

import rioxarray
import rasterio

band = rioxarray.open_rasterio(inputpath_raster)

band.plot()

print(band.rio.crs)

Reproject the raster to geographic coordinates

band_geographic = band.rio.reproject('EPSG:4326')

print(band_geographic.rio.crs)

band_geographic.plot()

Reproject the raster to UTM coordinates

band_utm = band_geographic.rio.reproject(band_geographic.rio.estimate_utm_crs())

print(band_utm.rio.crs)

band_utm.plot()

Save the raster

band_utm.rio.to_raster(outputpath_raster)

Answer 3

Here you go, you can try this.

I followed the tutorial you mentioned and I am using the same x y and ref variables. Your biggest challenge is making sure of the projection. From the tutorial it seems like it is an Lambert Conformal projection over Northern America. According to the tutorial's documentation (https://scitools.org.uk/cartopy/docs/latest/reference/projections.html#lambertconformal ) I assumed this one: https://epsg.io/102009 (your job to make 100% sure)

Here you can read about importing projections form the rasterio library https://rasterio.readthedocs.io/en/latest/api/rasterio.crs.html

import numpy as np
import rasterio
from rasterio.crs import CRS

#x ,y and ref variables from the provided tutorial

# Construct rasterio transform object using x and y coordiantes
xsize = abs(x[0,0] - x[0,1])
ysize = abs(y[0,0] - y[0,1])
transform = rasterio.transform.from_bounds(np.min(x), 
                                           np.min(y), 
                                           np.max(x), 
                                           np.max(y),
                                           width=xsize, 
                                           height=ysize)
# Set raster profile
profile = {
'driver': 'GTiff',
'dtype': np.float32,
'nodata': 9999,
'width': xsize,
'height': ysize,
'count': 1,
'crs': CRS.from_proj4('+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m no_defs'),
'transform': transform,
'tiled': True,
'compress': 'lzw'
}

# Write output raster
ref_filled = ref.filled(9999) # fill the masked array with the NoData value
with rasterio.open('path/to/my_radar_data.tif','w',**profile) as dst:
    dst.write(ref_filled,1)