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?

  • I think is better you reproject your raster in Qgis. Do you know what is your UTM Zone and EPSG?
    – Helios
    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
    – Helios
    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.
    – Pythonist
    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
    – user2856
    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.
    – user2856
    Mar 2, 2022 at 23:21

3 Answers 3


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):
    Writes a geotiff from a Numpy array with appended georeferencing from parent geotiff.
    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)
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:


    if nodata is not False:

    # Set transform and projection
    wkt = gdal_obj.GetProjection()
    srs = osr.SpatialReference()

    # Close output raster dataset 
    dest = None
  • 1
    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.
    – Pythonist
    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)



Reproject the raster to geographic coordinates

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



Reproject the raster to UTM coordinates

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



Save the 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.
    – Pythonist
    Mar 3, 2022 at 15:48
  • inputpath_raster that is your raster band in geotiff
    – Helios
    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.
    – Pythonist
    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
    – Helios
    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), 
# 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:

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.