1

Good Evening everyone. I have a question regarding the warping options in rasterio. I have two raster datasets from Albania. The first raster is the product of a clipping, and the second one is the country's boundaries. For the first, I have the following projection. enter image description here

For the second raster I have the following information: enter image description here

My goal is to reproject the raster#1 with the spatial information of raster 2. I made an attempt using rasterio, and the result shows the following spatial reference: Krassovsky_1942_Transverse_Mercator, and the Datum is unknown. This result was obtained after using the following code:

import numpy as np
import rasterio
from rasterio.warp import calculate_default_transform, reproject, Resampling


#-----------------------------------------------------------------------------  #
#                   OPENING IMAGE WITH THE CORRECT SPATIAL REFERENCE
#-----------------------------------------------------------------------------#
with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\raster\AL_020m_nat_buffer100m.tif") as proj:
dst_crs = proj.crs
print dst_crs
transform_proj, width_proj, height_proj = calculate

#-----------------------------------------------------------------------------#
#                   OPENING IMAGE WITH TO BE REPROJECTED
#-----------------------------------------------------------------------------#

with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\results\umprojected\ALB_adm0.tif") as src:
transform, width, height = calculate_default_transform(
    src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
    "crs": dst_crs,
    "transform": transform,
    "width": width,
    "height": height
})

#---------------------------------------------------------------------------
#                   CREATING THE OUTPUT IMAGE   
#---------------------------------------------------------------------------

with rasterio.open(r"C:\Users\Roger\Documents\git\ArcpyRecipes\Test_rasterio\test_albania\results\projected/Albita_projected_2.tif", "w", **kwargs) as dst:
    for i in range(1, src.count + 1):
        reproject(
            source = rasterio.band(src, i),
            destination=rasterio.band(dst, i),
            src_transform = src.transform,
            src_crs=src.crs,
            dst_transform = transform,
            dst_crs = dst_crs,
            resampling = Resampling.nearest)

I read that the EPGS code can be retrieved in ArcMap using arcpy.Describe, but I want to use the module rasterio to accomplish the task. I thought the option (dataset.crs) after opening an image with rasterio would give me the CRS to make the projection, but It was not the case. Could anyone help me with this?

Thank you very much

Roger Almengor
  • 442
  • 5
  • 16
  • I used the following post to retrieve the Well Known Text https://gis.stackexchange.com/questions/60371/gdal-python-how-do-i-get-coordinate-system-name-from-spatialreference. It retrieves the name of the Geographic Coordinate System and the Projected Coordinate System. The thing is, that in rasterio the warping documentation works with the EPGS codes and not with the names, as far as I know. – Roger Almengor May 05 '18 at 20:16
  • So your problem is that the CRS of the reprojected file is not the same as with the sample file you have, even though you got it from there? – j08lue Nov 06 '18 at 21:35
  • If so, I don't think it will help to get the CRS as WKT and then feed that back into rasterio. But since you asked -- yes you can create a `rasterio.crs.CRS.from_wkt`. – j08lue Nov 06 '18 at 21:40

0 Answers0