TransWikia.com

Raster has a "none" CRS

Geographic Information Systems Asked on May 31, 2021

My specifications :

  • Python 3.9.1
  • gdal 3.1.4
  • rasterio 1.1.8
  • Windows 10 Business
  • conda version 4.9.2

I have a raster that has no crs :

os.chdir(path)
file_list= os.listdir()
file_list[0]
>>> RGEALTI_FXX_0605_6670_MNT_LAMB93_IGN69.asc

with rasterio.open(liste_fichiers[0]) as src:
    print(src.profile)
    print(src.crs)
    print(src.bounds)
    print(src.indexes)
>>> {'driver': 'AAIGrid', 'dtype': 'float32', 'nodata': -99999.0, 'width': 1000, 'height': 1000, 'count': 1, 'crs': None, 'transform': Affine(5.0, 0.0, 604997.5,
       0.0, -5.0, 6670002.5), 'tiled': False}
None
BoundingBox(left=604997.5, bottom=6665002.5, right=609997.5, top=6670002.5)
(1,)

For your information, I do have a GDAL_DATA environment variable. (I found this information here).

When I open the raster with QGIS, I can in fact notice that there is no CRS, but on the other hand I can easily select one.
I am looking forward to do the same with rasterio. Here is the code I tried to do so :

with rasterio.open(liste_fichiers[0], mode='r+') as src:
    print(src.profile)
    src.crs = rasterio.crs.CRS.from_epsg(2154)
>>> {'driver': 'AAIGrid', 'dtype': 'float32', 'nodata': -99999.0, 'width': 1000, 'height': 1000, 'count': 1, 'crs': None, 'transform': Affine(5.0, 0.0, 604997.5,
       0.0, -5.0, 6670002.5), 'tiled': False}

---------------------------------------------------------------------------
CPLE_AppDefinedError                      Traceback (most recent call last)
<ipython-input-191-a0546e13d796> in <module>
      1 with rasterio.open(liste_fichiers[0], mode='r+') as src:
      2     print(src.profile)
----> 3     src.crs = rasterio.crs.CRS.from_epsg(2154)

rasterio_base.pyx in rasterio._base.DatasetBase.__exit__()

rasterio_base.pyx in rasterio._base.DatasetBase.close()

rasterio_io.pyx in rasterio._io.BufferedDatasetWriterBase.stop()

rasterio_io.pyx in rasterio._io._delete_dataset_if_exists()

rasterio_err.pyx in rasterio._err.exc_wrap_int()

CPLE_AppDefinedError: Deleting RGEALTI_FXX_0605_6670_MNT_LAMB93_IGN69.asc failed: Permission denied

2 Answers

Taking LAMB93_IGN69 in the file name as a clue the CRS is most likely to be RGF93 / Lambert-93 + NGF-IGN69 height a compound CRS which in the EPSG registry is assigned code 5698 (https://epsg.org/crs_5698/RGF93-Lambert-93-NGF-IGN69-height.html).

SCOPE: Engineering survey, topographic mapping.

EXTENT: France - mainland onshore

HORIZONTAL CRS: RGF93 / Lambert-93 (https://epsg.org/crs_2154/RGF93-Lambert-93.html)

VERTICAL CRS: NGF-IGN69 height (https://epsg.org/crs_5720/NGF-IGN69-height.html)

Correct answer by nmtoken on May 31, 2021

The ASCII Grid format projection must be stored in a .prj file, which you don't have so the grid hasn't a projection.

Also, the driver isn't in-place writeable.

import rasterio

dst_crs = 'EPSG:2154'

with rasterio.open('without_crs.asc', 'r') as src:
    raster = src.read()
    kwargs = src.meta.copy()
    kwargs.update({
        'crs': dst_crs
    })

with rasterio.open('assigned_crs.asc', 'w', **kwargs) as dst:
    dst.write(raster)

You will get an assigned_crs.prj with the .asc file, i.e. the new grid file has a projection.


Don't need to do the same with all files, just copy and rename the same .prj file.

Answered by Gabriel De Luca on May 31, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP