Geographic Information Systems Asked on May 31, 2021
My specifications :
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
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
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP