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?
