2

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

Basile
  • 617
  • 4
  • 15

2 Answers2

1

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)

nmtoken
  • 13,355
  • 5
  • 38
  • 87
1

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.

Gabriel De Luca
  • 14,289
  • 3
  • 20
  • 51