On suppose dans cette procédure que nos données constituent un maillage régulier et qu’il n’est pas nécessaire de faire d’interpolation. Le fichier ASCII XYZ est composé de trois colonnes X, Y et Z, ex.
| 880000 |
6340000 |
828 |
| 880025 |
6340000 |
831 |
| 880050 |
6340000 |
833 |
| 880075 |
6340000 |
835 |
Procédure avec ArcGIS
1-Changer l’extension du fichier en .txt ou .csv
2- Modifier le fichier schema.ini
Plus d’infos sur ce fichier schema.ini à http://www.tectonics.caltech.edu/gi…
Ex. de fichier schema.ini correspondant au fichier input.csv ci-dessus : [input.csv]
ColNameHeader=False
Format=Delimited( )
Col1=X Double
Col2=Y Double
Col3=Z Double
- si le fichier contient une ligne d’entête, on aura ColNameHeader=True
- si le séparateur est la virgule on aura Format=CSVDelimited
- si le séparateur est la tabulation on aura Format=TabDelimited
- si le séparateur est l’espace on aura Format=Delimited( )
- si le séparateur est | on aura Format=Delimited(|)
- … etc avec d’autres séparateurs
Si le fichier n’est pas bien reconnu par ArcCatalog, il est possible que ce soit un problème de codage des caractères. On pourra le modifier en ajoutant la ligne suivante dans schema.ini :
CharacterSet=UNICODE
Changer UNICODE par le codage correct (ANSI, autre …)
3- Vérifier la reconnaissance du fichier par ArcCatalog
- Vérifier dans ArcCatalog que votre fichier est bien reconnu. Pour cela votre fichier doit s’afficher dans l’onglet Aperçu.
- Vérifier aussi que le type des champs est correct. Pour cela faire un clic droit sur le fichier et vérifier l’onglet Champs
4- Créer le shapefile de points
- Pour créer les points, utiliser la fonction d’Arctoolbox « Outils 3D Analyst Depuis un fichier – ASCII 3D vers classe d’entités »
5-Convertir le shapefile de points en raster
- Utiliser la fonction d’ArcToolbox « Outils de conversion – Vers raster – Entité vers raster »
- Prendre le champ Shape.Z ou le champ Z comme valeur de raster
- Définir la taille de cellule selon l’espacement des points en entrée.
Procédure avec Linux et GDAL
ex. conversion d’un fichier XYZ ASCII en TIF et assignation de la référence spatiale RGF93-Lambert93 (Code EPSG 2154)
gdal_translate -a_srs EPSG:2154 input.XYZ output.tif
en batch on pourra faire un script du type :
#!/bin/bash
#Ce script assigne le code EPSG saisi en paramètre aux images en sortie
# Jules Fleury, SIGeo/CEREGE, mars 2011# saisir le code epsg en parametre
# Assign_SRS_to_TIF.sh
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326) : ”
read c_EPSG
for FILE in *.XYZ # modifier ici en cas d'autre format
do
BASE=`basename $FILE .XYZ` # modifier ici en cas d'autre format
NEWFILE=${BASE}_c.tif
gdal_translate -a_srs EPSG:$c_EPSG $FILE $NEWFILE
done
Avec Linux et GDAL : cas ou les données en entrée ne sont pas bien triées en X ou Y
solution prise ici : http://osgeo-org.1560.x6.nabble.com/gdal-dev-Opening-gridded-xyz-data-that-is-out-of-order-td5334341.html
if we obtain this kink of error:
> ERROR 1: Ungridded dataset: At line 4125001, change of Y direction
> gdalinfo failed – unable to open ‘DTM_swissALTI3D_XYZ.txt’.
You must sort the data beforehand
– the y position must be sorted properly (ASC or DESC)
Standard linux sort program:
sort -k2 -n -k1 input.xyz -o output_sort.xyz
(2nd column (y) as numeric, then first column (x) -o = output file
puis on relance gdal_translate
gdal_translate input.xyz output.tif -a_nodata -9999
Remplacer par la bonne valeur pour nodata
Si cette méthode ne fonctionne pas, on fera une interpolation comme ci-dessous
Avec Linux et GDAL : cas ou les données ne couvrent pas une dalle entière
On utilise la fonction d’interpolation de GDAL gdal_grid pour ne prendre pour chaque pixel qu’un seul point en entrée. La méthode consiste à premièrement créer un CSV, puis un VRT et utiliser ce fichier comme source de gdal_grid.
ex. gdal_grid -a nearest:radius1=1:radius2=1:nodata=-9999.0 -ot Float32 -a_srs EPSG:2154 -txe 1074000 1079000 -tye 6325000 6330000 -outsize 1000 1000 -l grid grid.vrt grid_out.tif
traitement en batch
1ère étape : conversion des XYZ en CSV avec création de VRT associé
ex. de script
#!/bin/bash
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) : ”
read c_EPSG
for IFile in *.xyz # modifier ici en cas d’autre format
do
BASE=`basename $IFile .xyz`
OFile=${BASE}.csv
sed ‘1i\X Y Z’ $IFile > $OFile
TFile=”Temp”
tr ‘ ‘ ‘\t’ < $OFile > $TFile
mv $TFile $OFile
#création du fichier VRT
OFileVRT=${BASE}.vrt
echo "<OGRVRTDataSource>" > $OFileVRT
echo " <OGRVRTLayer name=\"$BASE\">" >> $OFileVRT
echo " <SrcDataSource>$OFile</SrcDataSource>" >> $OFileVRT
echo " <GeometryType>wkbPoint</GeometryType>" >> $OFileVRT
echo " <LayerSRS>EPSG:$c_EPSG</LayerSRS>" >> $OFileVRT
echo " <GeometryField encoding=\"PointFromColumns\" x=\"X\" y=\"Y\" z=\"Z\"/>" >> $OFileVRT
echo " </OGRVRTLayer>" >> $OFileVRT
echo "</OGRVRTDataSource>" >> $OFileVRT
done
2ème étape : conversion des VRT (avec leur CSV) en TIF
on utilise gdal_grid,ex. de script
#!/bin/bash
#selon différents algorithmes et paramètres.
#le fichier grid en sortie est stocké en TIF
#puis rééchantilloné selon la résolution voulue.
#avril 2013
echo “Input spatial reference EPSG (exemple: for WGS84 enter 4326, for Lambert93 2154, for L2E 27572) :”
read c_EPSG
for IFile in *.vrt # modifier ici en cas d’autre format
do
echo “Processing $IFile”
BASE=`basename $IFile .vrt`
OFile=${BASE}.tif
val_nodata=-9999 #valeur pour nodata
gdal_grid -a nearest:radius1=0.5:radius2=0.5:nodata=$val_nodata -ot Float32 -a_srs EPSG:$c_EPSG -l $BASE $IFile $OFile
#ré-échantillonage du fichier en sortie
ressize=10 #valeur de ré-échantillonage
OFileRes=${BASE}_$ressize.tif #fichier ré-échantilloné en sortie
gdalwarp -tr $ressize $ressize -r cubic $OFile $OFileRes
echo "Done"
done
n -32768 : pour fixer les valeurs nodata