Problem with gdal.WarpOptions concerning reference

Question

I have a tiff file that I want to crop down to a shp cutline and reproject.

If the code is:

clip_spec=gdal.WarpOptions(format='GTiff',cutlineDSName=os.getcwd()+'/moksha/DEM/cutline2.shp',cropToCutline=True,dstNodata=-9999)
cutted=gdal.Warp(os.getcwd()+'/moksha/DEM/cutted22_m.tif',src,options=clip_spec,srcSRS='EPSG:4326',dstSRS='EPSG:32638')

I get a file with nodata.

If I include srs to options as:

clip_spec=gdal.WarpOptions(format='GTiff',cutlineDSName=os.getcwd()+'/moksha/DEM/cutline2.shp',cropToCutline=True,dstNodata=-9999,srcSRS='EPSG:4326',dstSRS='EPSG:32638')
cutted=gdal.Warp(os.getcwd()+'/moksha/DEM/cutted22_m.tif',src,options=clip_spec)

I get - TypeError: in method 'wrapper_GDALWarpDestName', argument 4 of type 'GDALWarpAppOptions *

Anywhere I looked on web - these options are treated correctly. gdalTranslate worked well for reprojecting, but I would like to use gdal.Warp, as I have to crop the file. Tried to type otherwise with **kwargs - same thing. What could be the problem here and how can I make it work finally?

Files are available through the link

UPD: The problem was only in Jupiter, commandline works nice.

Answer 1

Works for me with GDAL 3.1.4 and Python 3.7. I dropped the subdirectories from the path for simplicity.

from osgeo import gdal, osr
import os
src='merged22.tif'
clip_spec=gdal.WarpOptions(format='GTiff',cutlineDSName=os.getcwd()+'/cutline2.shp',cropToCutline=True,dstNodata=-9999,srcSRS='EPSG:4326',dstSRS='EPSG:32638')
cutted=gdal.Warp(os.getcwd()+'/cutted22_m.tif',src,options=clip_spec)

Test with gdalinfo

gdalinfo cutted22_m.tif -stats
Driver: GTiff/GeoTIFF
Files: cutted22_m.tif
Size is 4613, 7255
Coordinate System is:
PROJCRS["WGS 84 / UTM zone 38N",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["UTM zone 38N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",45,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - N hemisphere - 42┬░E to 48┬░E - by country"],
        BBOX[0,42,84,48]],
    ID["EPSG",32638]]
Data axis to CRS axis mapping: 1,2
Origin = (368922.764431581657846,6090209.774193479679525)
Pixel Size = (32.176682285386590,-32.173780105605999)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  368922.764, 6090209.774) ( 42d57'13.44"E, 54d56'29.78"N)
Lower Left  (  368922.764, 5856789.000) ( 43d 3'13.35"E, 52d50'41.55"N)
Upper Right (  517353.800, 6090209.774) ( 45d16'15.63"E, 54d57'30.69"N)
Lower Right (  517353.800, 5856789.000) ( 45d15'27.92"E, 52d51'37.97"N)
Center      (  443138.282, 5973499.387) ( 44d 8' 4.06"E, 53d54'24.48"N)
Band 1 Block=4613x1 Type=Int16, ColorInterp=Gray
  Minimum=94.000, Maximum=298.000, Mean=188.517, StdDev=38.004
  NoData Value=-9999
  Metadata:
    STATISTICS_MAXIMUM=298
    STATISTICS_MEAN=188.51715516792
    STATISTICS_MINIMUM=94
    STATISTICS_STDDEV=38.004354639842
    STATISTICS_VALID_PERCENT=95.22

The image has bit depth of 16 bits and maximum value is 298 so without histogram stretch it will appear as black. But with min/max stretch QGIS shows it fine. The blurred background is the source, merged22.tif.

Answer 2

A solution:

from osgeo import gdal
source = gdal.Open('merged22.tif')
# Projection of the source raster:
epsg = int(gdal.Info(source, format='json')['coordinateSystem']['wkt'].rsplit('"EPSG","', 1)[-1].split('"')[0])
print(epsg)
4326
#Crop
output= 'cutted22_m.tif'
cutted = gdal.Warp(output,
          source,
          format = 'GTiff',
          cutlineDSName = 'cutline2.shp',
          dstNodata = -9999,
          dstSRS= 'EPSG:32638')

And it works without error (polygon shapefile in red, cutted raster in grey, source raster in colors)

Projection of the saved cutted raster:

cutted_res= gdal.Open('cutted22_m.tif')
epsg = int(gdal.Info(cutted_res, format='json')['coordinateSystem']['wkt'].rsplit('"EPSG","', 1)[-1].split('"')[0])
print(epsg)
32638