0

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.

0
0

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.

enter image description here

0

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)

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.