0

I'm using Python to process rasters that have a CRS. I want to copy the CRS to the raster that is produced.

The rasters are like the following. Sorry, I don't know how to create this from scratch, still learning how to do things in Python (these are mini versions):

tmin:
<xarray.DataArray (band: 2, y: 1, x: 1)>
[2 values with dtype=float32]
Coordinates:
  * band         (band) int32 1 2
  * x            (x) float64 1.426e+06
  * y            (y) float64 -3.82e+05
    spatial_ref  int32 0
Attributes:
    AREA_OR_POINT:       Area
    STATISTICS_MAXIMUM:  20.5
    STATISTICS_MEAN:     -9999
    STATISTICS_MINIMUM:  20.5
    STATISTICS_STDDEV:   -9999
    _FillValue:          nan
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           ('tmin_100', 'tmin_101')
tmin.values
array([[[3.5]],

       [[6. ]]], dtype=float32)
tmax:
<xarray.DataArray (band: 2, y: 1, x: 1)>
[2 values with dtype=float32]
Coordinates:
  * band         (band) int32 1 2
  * x            (x) float64 1.426e+06
  * y            (y) float64 -3.82e+05
    spatial_ref  int32 0
Attributes:
    AREA_OR_POINT:       Area
    STATISTICS_MAXIMUM:  3.5
    STATISTICS_MEAN:     -9999
    STATISTICS_MINIMUM:  3.5
    STATISTICS_STDDEV:   -9999
    _FillValue:          nan
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           ('tmax_100', 'tmax_101')
tmax.values
array([[[20.5]],

       [[19. ]]], dtype=float32)

A simplified version of the code is

import rasterio
import rioxarray as rxr
import numpy as np
import os

# I get an error about proj.db not found without this.
os.environ['PROJ_LIB'] = 'C:/OSGeo4W64/share/proj/'

# Simplified processing
def calculate_mean(tmin, tmax):
    # Read files
    tmax_src = rxr.open_rasterio(tmax_file)
    tmin_src = rxr.open_rasterio(tmin_file)
    # Change to np.array for calculate_mean
    tmin = tmin_src.values
    tmax = tmax_src.values
    rascrs = tmax_src.rio.crs
    
    result1 = (tmax + tmin) / 2
    return result1


    
tmax_file = "Q:/Shared drives/Data/tmax9.tif"
tmin_file = "Q:/Shared drives/Data/tmin9.tif"


rascrs = rxr.open_rasterio(tmax_file)   
rascrs = rascrs.rio.crs.wkt

tempraster = calculate_mean(tmax_file, tmin_file)

# Define the output file path
output_file = "Q:/Shared drives/Data/results/rhoutput.tif"

# Get the number of layers (365)
num_layers = tempraster.shape[0]

# Get the height and width of the raster
height, width = tempraster.shape[1], tempraster.shape[2]

# Create the output file using rasterio
with rasterio.open(output_file, 'w', driver='GTiff', height=height, width=width, count=num_layers, dtype=tempraster.dtype, crs=rascrs) as dst:
    # Loop through each layer and write it to the file
    for layer in range(num_layers):
        dst.write(tempraster[layer], layer + 1)  # Layers are 1-indexed in rasterio

The CRS:

rascrs 'PROJCS["unnamed",GEOGCS["unknown",DATUM["unnamed",SPHEROID["Spheroid",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",42.5],PARAMETER["central_meridian",-100],PARAMETER["standard_parallel_1",25],PARAMETER["standard_parallel_2",60],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'

When I try to write the raster to file I get this error:

C:\Users\japolo\Documents\code\pipenvs\pythonrast.venv\lib\site-packages\rasterio_init_.py:329: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned. dataset = writer(

I've seen other questions and resources that say to use GDAL library for Python. But I can't install GDAL in my Python environment because of some kind of problem due to Microsoft C++ versions. ( I have 2015-2019 and the error says "2014 or greater". Go figure. ) It seems like the crs=rascrs should work, but since it's not, my new raster files lack the geographic information. Is there a way to fix this?

1 Answer 1

0

A CRS alone is not enough to georeference a raster, you also need to specify the transform when opening the output file, by copying that from the input: https://rasterio.readthedocs.io/en/stable/api/rasterio.html#rasterio.open

1
  • The code would be something like ras_trans = rascrs.rio.transform() and the write statement like: with rasterio.open(output_file..., crs=crsras, transform=ras_trans).
    – John Polo
    Oct 3, 2023 at 18:14

Your Answer

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

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