Suite

Tuiles raster à partir de données vectorielles avec GDAL - comment éviter les artefacts d'aliasing ?

Tuiles raster à partir de données vectorielles avec GDAL - comment éviter les artefacts d'aliasing ?


Prérequis

J'ai des données vectorielles dans un ShapefileN47E007_CONTOURS.shpproduit avec

wget -O N47E007.zip http://e4ftl01.cr.usgs.gov/SRTM/SRTMGL1.003/2000.02.11/N47E007.SRTMGL1.hgt.zip unzip N47E007.zip # Convertir de WGS 84 en Web Mercator et de HGT to GeoTIFF gdalwarp -s_srs EPSG:4326 N47E007.hgt  -t_srs EPSG:3785 N47E007_WebMercator.tif  -co "COMPRESS=LZW" -of GTiff -r bilinéaire # Calculer les données vectorielles de contour gdal_contour -a elev -i 20 N47E007_WebMercator.tif .shp

(Lignes de contour basées sur un patch (N47E007) des données d'élévation NASA SRTM.)

Traiter

À partir de là, je génère des tuiles TMS/WMTS à utiliser par dépliant, en rastérisant d'abord les données vectorielles du Shapefile avecgdal_rasterizepuis le carreler avecgdal2tiles.py. Je répète les deux étapes pour chaque niveau de zoom, afin que je puisse pixelliser vers un plus grand raster pour des niveaux de zoom plus élevés, de sorte que les lignes aient la même épaisseur à l'écran sur chacune :

#!/bin/sh pour z dans $(seq 1 15) do gdal_rasterize  --config GDAL_CACHEMAX 1024  -ts $(calc "2851 << (${z} - 12)")  $(calc "4220 < < (${z} - 12)")  -ot byte -of GTiff -co COMPRESS=DEFLATE -co PREDICTOR=2 -co ZLEVEL=9  -ot Byte -burn 0 -a_nodata 255  -l N47E007_CONTOURS N47E007_CONTOURS.shp N47E007_CONTOURS_rasterized_z${z}.tif gdal2tiles.py  -r près de -z${z} -a 255 255 255  -e N47E007_CONTOURS_rasterized_z${z}.tif  N47E007_contours_tiled done

Adapter la résolution au niveau de zoom

(Voir aussi ma réponse à Problème d'épaisseur des lignes de contour dans des niveaux de zoom plus importants.)

calc "2851 << (${z} - 12)"etcalc "4220 << (${z} - 12)"calculer la taille de l'image à pixelliser (N47E007_CONTOURS_rasterized_z${z}.tif) en fonction du niveau de zoom$z.

<<est l'opérateur binaire de décalage à gauche, donc cela se traduit par 2851⋅2 z-12 et 4220⋅2 z-12. Cela signifie que

  • pour le niveau de zoom 12, l'image raster intermédiaireN47E007_CONTOURS_rasterized_z12.tifsera de 2851×4220 pixels,
  • pour le niveau de zoom 13 (N47E007_CONTOURS_rasterized_z13.tif) ce sera le double dans chaque sens (5702×8440) et
  • pour le niveau de zoom 11 (N47E007_CONTOURS_rasterized_z11.tif) la moitié de celle du niveau de zoom 12, soit 1425×2110.

Paramétrage imprécis

Les nombres2851et4220(la taille du zoomelevel 12) ont été devinés ou mesurés à l'écran ou autrement dérivés de manière imprécise.

Problème

Les tuiles résultantes montrent des artefacts qui ne sont pas visibles dans les images rastérisées intermédiaires (non tuilées) :

  • segments de lignes à double largeur,
    par exemple.13/4270/5331.png"> vs.

  • lacunes dans les lignes,
    par exemple. sur les bords des carreaux entre

    14/8536/10665.png">

    extrait correspondant deN47E007_CONTOURS_rasterized_z14.tifen comparaison:

Si je saute-r près(donc la rastérisation utilise le mode de rééchantillonnage par défaut 'moyenne') À la place, j'obtiens des pixels délavés qui rendent les lignes floues.

Ma théorie

Les systèmes de coordonnées de référence correspondent déjà (en raison de lagdalwarpétape avant de calculer les contours), de sorte qu'aucune transformation/reprojection générale ne devrait être requise. Mais comme la taille à rastériser n'est qu'approximative de la taille à l'écran de la zone géographique couverte par leN47E007correctif de données, un redimensionnement doit encore se produire dansgdal2tiles.py, pour produire des carreaux sans soudure correctement placés et dimensionnés.

Avec-r prèscela provoque une sorte d'effet de moiré macroscopique, car les grilles de pixels des tuiles et de la trame intermédiaire ne correspondent presque qu'à peu. (Aditionellement,gdal2tiles.pypeut ne pas regarder à travers les bords de tuiles lors de la recherche du pixel le plus proche dans l'image source, provoquant des lacunes.) Avecmoyenne(la valeur par défaut pour-r), les mêmes causes affectent l'aspect délavé car les valeurs de plus d'un pixel source sont « mélangées » pour obtenir la valeur d'un pixel de tuile donné.

Ainsi, si je connaissais la taille correcte à l'écran (en fonction du niveau de zoom) de la zone géographique du patch de données et que je l'utilisais pour la rastérisation, les artefacts devraient disparaître.

Des questions

  • Ma théorie est-elle correcte ?
  • Si oui, comment puis-je déterminer la précision-tsparamètres pourgdal_rasterize(Ou bien,-trparamètres) de sorte qu'aucun rééchelonnement (ou un rééchelonnement sans effet) ne sera effectué pargdal2tiles.py?

Ma théorie est-elle correcte ?

Il semble que oui, car j'obtiens des résultats impeccables avec la méthode décrite ci-dessous.

Si oui, comment puis-je déterminer le-tsparamètres pourgdal_rasterize(Ou bien,-trparamètres) de sorte qu'aucun rééchelonnement (ou un rééchelonnement sans effet) ne sera effectué pargdal2tiles.py?

Tandis quegdal2tiles.pyest destiné à être utilisé comme une commande, nous pouvons le charger en tant que module python et utiliser les mêmes fonctions qu'il utilise lui-même pour les calculs de taille. Ainsi, je me suis fait ce petit script python :

de __future__ importer print_function importer sys de gdal2tiles importer GlobalMercator def main() : lat, lon = [ float(arg) for arg in sys.argv[1:3] ] zoom = int(sys.argv[3]) px, py = lat_lon_to_pixels(lat, lon, zoom) print(px, py) def lat_lon_to_pixels(lat, lon, zoom): gm = GlobalMercator() mx, my = gm.LatLonToMeters(lat, lon) return tuple( int(round(p )) pour p dans gm.MetersToPixels(mx, my, zoom) ) if __name__ == '__main__': main()

J'utilise ensuite les coordonnées WGS 84 degdalinfo N47E007.hgt(en les extrayant avec les expressions régulières/^s*En haut à droites*(s*(?d+.?d*),s*(?d+.?d*))et/^s*En bas à gauches*(s*(?d+.?d*),s*(?d+.?d*))/) appelerlat_lon_to_pixelsà partir du script python ci-dessus pour le coin nord-est et le coin sud-ouest et calculez les différences, que j'utilise ensuite comme argument pour le-tspossibilité degdal_rasterize. (Cela ne fonctionne que parce que ce fichier HGT a WGS 84 comme système de référence de coordonnées et WGS 84 et WebMercator sont alignés les uns sur les autres.)


Voir la vidéo: Stockage, manipulation et analyse de données matricielles avec PostGIS Raster