1

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

7
  • 1
    Try reading the contents of the PRJ file and passing it to: rasterio.readthedocs.io/en/latest/api/… Then after you open the file, see if you can set that as the crs attribute
    – mikewatt
    Commented Jul 25, 2022 at 20:15
  • And are the dimensions of the raster (rows & columns or height & width) 472 & 545? Looks to me like your raster has no georeferencing (which is different from projection info).
    – user2856
    Commented Jul 25, 2022 at 23:28
  • You're using the 9x9 tif which is not georeferenced correctly. Neither is the 3x3 tif. The 1x1 tif is georeferenced correctly though.
    – user2856
    Commented Jul 25, 2022 at 23:36
  • @user2856 thanks for this! Can I ask - how were you able to tell that the 9x9/3x3 are not georeferenced correctly?
    – sm2347
    Commented Jul 26, 2022 at 14:21
  • @user2856 sorry I hadn't seen your previous edit. For the 9x9, yes the row x column dims are 472 and 545. I wasn't aware that the georeferencing is separate from projection info.
    – sm2347
    Commented Jul 26, 2022 at 14:28

1 Answer 1

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

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.