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

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_data`with 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?

• I think is better you reproject your raster in Qgis. Do you know what is your UTM Zone and EPSG? Mar 2, 2022 at 18:32
• it's a little complicated reproject a raster, look at this: github.com/rasterio/rasterio/blob/master/examples/reproject.py Mar 2, 2022 at 18:36
• I know how to write a raster, the data or at least. The problem is how to embed the geographic metadata. Mar 2, 2022 at 18:43
• Embedding georeferencing when writing a geotiff from a numpy array is exactly the same regardless of projection, you build a `Affine` transform. There's not much more we can say unless you edit your question and include more details about your arrays, in particular array shape and x,y coordinates of the bounds (upper left, lower right) which you can extract from your arrays of coordinates. Note if the array is actually in lamberts, then these coordinates won't/shouldn't be lon, lat but x, y Mar 2, 2022 at 20:53
• @SalimRodríguez there is no need to ask if the poster agrees. If you have an answer, add it. Mar 2, 2022 at 23:21

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)
height = np.shape(array)

# 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
``````
• Thanks, but assuming that I have another raster with the CRS is too much. Getting he CRS, building it from scratch, is actually my problem. Anyway I have edited my question trying to be more precise on where I'm and what I want. Mar 3, 2022 at 15:30

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)
``````
• Thanks, that looks promising. But I don't get it. What is `inputpath_raster`? My starting point are either an object of the `siphon.radarserver` library, or numpy arrays with distances. Mar 3, 2022 at 15:48
• inputpath_raster that is your raster band in geotiff Mar 3, 2022 at 15:50
• Sorry but we are not understanding each other. I do not have a raster. I want to create a raster. I have an object of `siphon.radarserver` I want to convert in GeoTIFF. Mar 3, 2022 at 15:51
• Reprojecting is a lossy process, this does not sound like a very good idea. Mar 3, 2022 at 16:49
• why do you think Reprojecting is a lossy process? I would like to know your opinion Mar 3, 2022 at 23:13

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