gdalwarp produisant des rasters vides lorsqu'il est appelé à partir d'un script python, mais pas via la ligne de commande

8

J'essaie d'agréger certaines données raster (toutes .tif) jusqu'à une résolution plus grossière (de 0,05 degré à 0,25 degré) en utilisant gdalwarp en python, mais la commande ne fonctionne pas. Au lieu d'obtenir une sortie avec une large plage de valeurs, toutes les valeurs de la sortie sont 0. La résolution et la profondeur / le type de pixel sont corrects, mais les valeurs ne le sont pas.

Voici la documentation de la commande gdalwarp: http://www.gdal.org/gdalwarp.html

J'ai deux fichiers d'entrée que je souhaite agréger jusqu'à une résolution de 0,25 degré, produisant plusieurs sorties:

  • 'NDVI_raster': La première entrée est un raster signé 16 bits représentant NDVI, avec des valeurs allant de -10 000 à 10 000 et des valeurs nodata de -15 000.

  • 'nodata_mask': le second est un masque NoData, flottant 32 bits, où 1 = valeurs de "bonnes" données et 0 = NoData.

Avec «NDVI_raster» en entrée, je veux produire 7 sorties différentes, chacune représentant une statistique différente. Pour ce faire, j'appelle gdalwarp 7 fois, en définissant à chaque fois la méthode de rééchantillonnage (-r) sur l'une des valeurs suivantes: moyenne, mode, max, min, médiane, q1, q2. J'appellerai les sorties NDVI_ave.tif, NDVI_mode.tif, etc. En ce moment, j'utilise GDAL 1.10.1, qui n'autorise que la moyenne et le mode, donc je teste avec seulement ces deux statistiques maintenant.

En utilisant 'nodata_mask' comme entrée, je veux finalement produire un QAL (couche d'assurance qualité). Pour ce faire, j'utilise gdalwarp, avec le mode de rééchantillonnage réglé sur «moyenne» pour agréger jusqu'à 0,25 degré. Il en résulte que chaque pixel représente le rapport entre les bons pixels et le nombre total de pixels de l'entrée. Appelons la sortie QAL.

Voici ce qui est dans mon code (en utilisant le mode comme exemple pour la première entrée):

os.system('gdalwarp -tr .25 .25 -r mode -srcnodata -15000 %s %s' % (NDVI_raster, NDVI_mode))

Et pour la couche QA:

os.system('gdalwarp -tr .25 .25 -r average -srcnodata -15000 %s %s' % (nodata_mask, QAL))

Les résultats sont des rasters avec la résolution, la projection et la profondeur de pixels correctes, mais les valeurs des pixels sont toutes égales à 0.

Quiconque connaît python / gdal sait ce qui se passe?


Lorsque j'appelle la commande gdalwarp à partir de la ligne de commande (linux), j'obtiens le résultat souhaité. Lorsque j'utilise os.system pour appeler gdalwarp à partir de python, j'obtiens des rasters vides. Alors peut-être que quelque chose ne va pas avec mes liaisons gdal / python?


Au lieu d'appeler la commande via os.system, j'ai utilisé un sous-processus. L'outil via cette méthode semble également fonctionner correctement, mais le résultat est le même: un raster plein de 0.


J'ai essayé de mettre l'appel gdalwarp dans un script shell bash et d'appeler ce script shell à partir de python, mais le résultat est un tas de -1s au lieu de 0s. Curieusement, je l'avais déjà testé et je suis sûr que cela a fonctionné, mais le test a été supprimé de mon serveur et maintenant je ne peux pas le recréer pour une raison quelconque.


Mettre la commande gdalwarp dans un script shell bash puis appeler ce script shell à partir de la ligne de commande me donne le résultat souhaité. Cependant, appeler le même script shell à partir de python ne le fait pas. Il semble que quelque chose ne va pas avec python, mais comment et comment le corriger?

user20408
la source
2
Cela ne peut pas être un problème de liaison, vous ne les avez pas utilisés. Vous exécutez la commande gdalwarp. Vous devriez vérifier l'état d'erreur de la commande gdalwarp, voir: stackoverflow.com/questions/3791465/…
Zoltan
2
Pour développer ce que dit Zoltan, vous appelez essentiellement le même programme (gdalwarp). os.systemdémarre en fait un shell et exécute la commande donnée. En plus de la valeur de retour, vous devriez probablement vérifier les variables que vous essayez de transmettre (c.-à-d. Des problèmes de citation? Etc.)
Evil Genius
Vous devriez également envisager de passer du module le plus récent os.systemà celui subprocessqui contient un certain nombre d'améliorations (de sécurité).
Kersten
@Zoltan quel message d'erreur? Je n'ai pas reçu de message d'erreur car la commande fonctionnait techniquement, juste la sortie n'était pas correcte (valeurs de 0) .plz voir ma mise à jour. Merci!
user20408
2
Un problème est que je pense que vous devrez spécifier -ot Float32ou quelque chose de similaire pour le masque de qualité, car vous ne pouvez pas représenter des fractions avec un raster entier.
Nat Wilson

Réponses:

2

Vous ne dites pas comment démarrer votre script python / gdalwarp. J'ai découvert qu'un cronjob n'aura pas toujours le même environnement que mon environnement de ligne de commande. J'ai dû commencer à créer des environnements d'exécution pour ces types de scripts. Cela étant dit, si vous démarrez votre script à partir, disons, d'une icône sur votre bureau, il peut ne pas avoir le même environnement d'exécution que votre environnement de ligne de commande. "PYTHONPATH" peut être l'une des variables d'environnement que vous devez définir. De plus, vous devrez peut-être définir des variables pour gdalwarp. Enfin, vos fichiers de données peuvent ne pas être au bon endroit. Vous devrez peut-être définir un chemin absolu tel que / xxx / xxxx / NDVI_raster ou utiliser tilda ~ / NDVI_raster. Comme PYTHONPATH, vous devrez peut-être également configurer PATH et d'autres variables d'environnement."exporter ou source" les mêmes paramètres au début de votre script.

Greg
la source
1

J'ai eu ce problème également. J'ai finalement compris dans mon cas que c'était parce que je venais de créer l'image sur le disque à l'aide de gdalliaisons Python , mais je n'avais pas fermé l' gdal.Datasetobjet en mémoire, de sorte que l'écriture sur le disque n'était que partiellement terminée. Assez bizarrement, la seule façon que j'ai pu trouver pour fermer un gdal.Dataseten Python est: del variable_name_of_dataset- si moche!

Il Cexiste une GDALClose()méthode qui, pour le moment, n'est pas implémentée par l'API Python GDAL, mais le Rasteriofait: https://github.com/sgillies/rasterio/blob/876b9a1e2bf04e349b485e05ebc4a8674ace3cf0/rasterio/_io.pyx#L1463

Voir aussi: Pourquoi fermer un ensemble de données dans GDAL Python?

Hazzles
la source