0

I'm trying to resample an ALOS image of 12.5 m resolution to 1 m resolution with cubic interpolation using GDAL warp.

Here is the current code used in pyzo

import os
from os import path
import subprocess
import numpy as np
import numpy.ma as ma
import shutil

print (os.path.basename(__file__))
path0 = os.path.dirname(__file__)
print (path0)

# chemin vers les exécutables OTB, GDAL et OGR
path_otb = "C:/OTB-6.6.0-Win64/bin/"
path_gdal = "C:/OSGeo4W64/OSGeo4W.bat ; C:/OTB-6.6.0-Win64/bin/"
path_gdal_py = "C:/OSGeo4W64/OSGeo4W.bat ; python C:/OSGeo4W64/bin/"
path_ogr =  "C:\\OSGeo4W64\\OSGeo4W.bat ; C:/OSGeo4W64/bin/"

# chemin vers les données in et vers les résultats
path_data = path0.replace("prog","data")
path_out = path0.replace("prog","output")

image1 = path_data + "/AP_12775_FBD_F7170_RT1.dem.tif"
image2 = path_out + "/AP_12775_FBD_F7170_RT1_1M.dem.tif"
image3 = path_out + "/AP_12775_FBD_F7170_RT1_WM_2.dem.tif"

SRS33S = "EPSG:32733"
bbox= " 228444.7031250000000000 9901358.7500000000000000 308944.7031250000000000 9973996.2500000000000000"

if os.path.isfile(image1):
     os.remove(image1)
cmdline = "%sgdalwarp %s %s -te %s -tr %s -t_srs %s -r %s " \
             %(path_gdal, image1, image2, bbox, "1 1", SRS33S,
"cubic" )
return_code = subprocess.check_output(cmdline, shell=True)
print(return_code)

It returns me an error message saying :

"Traceback (most recent call last):
  File "C:\Users\Guillaume\Downloads\Cours_M2_TLD\TFE\prog\tp_gdal.py", line 553, in <module>
    return_code = subprocess.check_output(cmdline, shell=True)
  File "c:\python27\lib\subprocess.py", line 223, in check_output
    raise CalledProcessError(retcode, cmd, output=output)
CalledProcessError: Command 'C:/OSGeo4W64/OSGeo4W.bat ; C:/OTB-6.6.0-Win64/bin/gdalwarp C:\Users\Guillaume\Downloads\Cours_M2_TLD\TFE\data/AP_12775_FBD_F7170_RT1.dem.tif C:\Users\Guillaume\Downloads\Cours_M2_TLD\TFE\output/AP_12775_FBD_F7170_RT1_1M.dem.tif -te  228444.7031250000000000 9901358.7500000000000000 308944.7031250000000000 9973996.2500000000000000 -tr 1 1 -t_srs EPSG:32635 -r cubic ' returned non-zero exit status 2"

I cannot find the source of error, I have the same trouble using simple GDAL info command...

Do you have any idea how I can solve this problem?


The installation Python with Arcgis made my shell works, I still don't understand but it may help someone. Here is the final code that's working. Except forthe DEM slope results that are weird

#-*- coding: utf-8 -*-
#*************************************


import os
from os import path
import subprocess
import numpy as np
import numpy.ma as ma
import shutil

#pip install --upgrade pip
#install dbf

print (os.path.basename(__file__))
path0 = os.path.dirname(__file__)
print (path0)


# chemin vers les exécutables OTB, GDAL et OGR
path_otb = "C:/OTB-6.6.0-Win64/bin/"
path_gdal = "C:/OSGeo4W64/OSGeo4W.bat ; C:/OTB-6.6.0-Win64/bin/"
path_gdal_py = "C:/OSGeo4W64/OSGeo4W.bat ; python C:/OSGeo4W64/bin/"
path_ogr =  "C:\\OSGeo4W64\\OSGeo4W.bat ; C:/OSGeo4W64/bin/"

# chemin vers les données in et vers les résultats
path_data = path0.replace("prog","data")
path_out = path0.replace("prog","output")

    #******************************************************************************
# 17. GDAL_WARP
#  Rééchantillonnage d'image alos à des résolutions différentes
#  Choisir la bbox selon la bbox de l'image de base (aller voir dans image alos
#******************************************************************************

## L'emprise est trop grande, à l'aide de la calculatrice raster on peut la réduire
## en fonction de l'emprise de l'AAC

gr_alos_in = path_data + "/gr_alos_emprise_aac.tif"
gr_alos_out = path_out + "/gr_alos_1M.tif"

## GDALINFO pour la bbox
if 1==1:
print ("GDALINFO")

cmdline = path_gdal + "gdalinfo -stats " + gr_alos_in
return_code = subprocess.check_output(cmdline, shell=True)

print (return_code)

SRS33S = "EPSG:32733"

bbox= " 256788.0000000000000000 9923467.0000000000000000 271140.0000000000000000 9932682.0000000000000000"

# exemple 4.2 : masquer une image aux limites
# d'une couche de polygones ("clip raster" de ArcGIS)


## Rééchantillonnage

if os.path.isfile(gr_alos_out):
 os.remove(gr_alos_out)
cmdline = "%sgdalwarp %s %s -te %s -tr %s -t_srs %s -r %s " \
         %(path_gdal, gr_alos_in, gr_alos_out, bbox, "1 1", SRS33S,
"cubic" )
return_code = subprocess.check_output(cmdline, shell=True)
print(return_code)


#******************************************************************************
# 14. GDAL_DEM
#  Tools to process and analyze DEMs
#******************************************************************************

# input data
dem1 = path_out + "/gr_alos_1M.tif"

# 14.1 Slope
gr_slope = path_out + "/slope.tif"
cmdline = "%sgdaldem slope %s %s " \
       %(path_ogr,dem1,gr_slope)
return_code = subprocess.check_output(cmdline, shell=True)
print(return_code)
  • 1
    What is your goal, besides increasing the data volume by a factor of 150+? I wouldn't trust the values generated by cubic interpolation over such a distance. – Vince Feb 11 at 13:15
  • It's a M2 project in forestry. I want to assess the quality of DEM alos images compared to field value in Gabon. We did a first assessment using confusion matrices and linear regressions. But over the 96 values we gathered, we could only use half of them because the resolution of the raster (1 value for every 12.5 pixels) was lower than our study resolution (1 value each 10 m in a 150 m transect) wich led to errors in the results. – Guillaume Neve Feb 11 at 14:06
  • I do not have to use all the raster images, only the value contained in a certain region corresponding to the forest exploitation. I could change the bbox fitting this concession in order to reduce the amount of data – Guillaume Neve Feb 11 at 14:09
  • Resampling your ALOS to 1m will give you very little usable data. Also, be aware of what a radar 12.5m pixel is measuring over forestry, especially as radar resolution is different to optical resolution – Nathan Thomas Feb 18 at 4:30
  • Thanks for your answers @Vince The initial problem was the installation of Python, installing ArcGis and the linked python made my shell working... knowing that the Python versions were the same I don't know what has happen. – Guillaume Neve Feb 18 at 10:39

Your Answer

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

Browse other questions tagged or ask your own question.