Making Rasterio read CRS directly from PRJ file when CRS missing from metadata but present in PRJ file

Question

I'm trying to read a Raster file (GeoTiff) using Rasterio/Python. The file has fully populated metadata, but no CRS projection. This is present in a separate PRJ file (other files included are .aux and .rrd). The PRJ file mentions a User Defined Lambert Azimuthal Equal Area projection. Rasterio doesn't use the projections from .prj while reading the .tif file -- it just shows the CRS as None.

I've tried just entering the path to the PRJ file under the crs option in rasterio.open(file, crs) and that didn't work either.

Is there a way to make Rasterio read this information directly from the PRJ file?

I've seen a number of Q&As where the CRS may be clear from the PRJ file, but it's not obvious to me what the EPSG codes would be for a User Defined projection.

Extent of the raster:

>>>raster.bounds
BoundingBox(left=0.0, bottom=472.0, right=545.0, top=0.0)

Contents of the PRJ file: PROJCS["User_Defined_Lambert_Azimuthal_Equal_Area",GEOGCS["GCS_User_Defined",DATUM["D_User_Defined",SPHEROID["User_Defined_Spheroid",6378137.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",5071000.0],PARAMETER["False_Northing",3210000.0],PARAMETER["Central_Meridian",20.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]]

For what it's worth, this is the source - https://www.eea.europa.eu/data-and-maps/data/digital-elevation-model-of-europe

Answer 1

Your 9x9 (and the 3x3) tif are not georeferenced, the bounds/extent aren't in map coordinates, they are just the dimensions (columns, rows or pixel coordinates).

You can georeference using the 1x1 georeferencing using rasterio's rasterio.transform.from_gcps(gcps) method:

import rasterio as rio
from rasterio.transform import from_gcps
from rasterio.control import GroundControlPoint as GCP

# origin (upper left) map coordinate from elevation1x1.tif in x, y order 
xmin, ymax = (2399843.1035138000734150, 5500900.3833563001826406)  

cell_size = 9000

# Read the projection
prj = open("elevationeurope.prj").read()

with rio.open('elevation9x9.tif', "r+") as src: # open in r+ "append" mode
    # Get the georeferenced extent
    xmax = xmin + src.width * cell_size
    ymin = ymax - src.height * cell_size

    ul = xmin, ymax
    ll = xmin, ymin
    ur = xmax, ymax
    lr = xmax, ymin

    # calculate the geotransform using GCPs
    gcps = [
        GCP(0, 0, *ul),
        GCP(0, src.width, *ur),
        GCP(src.height, 0, *ll),
        GCP(src.height, src.width, *lr)
    ]

    transform = from_gcps(gcps)

    # Write the CRS and transform to the raster
    src.crs = prj
    src.transform = transform