Suite

Pourquoi gdalwarp renvoie-t-il des coordonnées de coin incorrectes pour mon géotiff ?


J'ai une image source png qui est dans la zone égale azimutale Lambert et j'ai l'étendue long/lat. J'ai donc exécuté gdal_translate pour créer des points de liaison et les convertir en gtif. Ensuite, j'ai exécuté gdalwarp pour reprojeter sur Web Mercator. gdalinfo montre que les coordonnées de délimitation sont correctes pour la sortie gdal_translate, mais pas pour la sortie gdalwarp. Que pasa ?

Voici mon opération gdal_translate sur le fichier source :

gdal_translate -of GTiff -a_srs "+proj=laea +lat_0=40. +lon_0=-105 +x_0=0 +y_0=0 +a=6370997 +rf=298.257" ^ -gcp 0 0 -119.344674 47.148697 ^ -gcp 525 0 -90.655326 47.148697 ^ -gcp 525 431 -93.609897 31,438244 ^ -gcp 0 431 -116.390103 31,438244 ^ -a_ullr -119.344674 47.148697 -93.609897 31,438244 ^ src.png">

La raison pour laquelle vous semblez obtenir de très petites valeurs est que votre image d'origine est projetée, elle utilise donc des unités linéaires - qui, dans ce cas, PROJ4, je pense, sont par défaut des mètres - plutôt que des degrés. Vous spécifiez donc une partie d'environ 680 mètres carrés de l'ouest des États-Unis.

Maintenant, comment déterminer les bonnes coordonnées à utiliser est une question délicate car votre image couvre une zone si vaste que vous obtiendrez beaucoup de distorsion sur les bords.

Désolé, ce n'est pas vraiment une solution - quelqu'un d'autre veut prendre le relais ? Si j'ai un éclair d'inspiration, je le posterai ici.


Figure 1. Exemple de zone de couverture et d'étendue spatiale du domaine nord-américain pour les prévisions de fumée RAP. Bientôt sur RealEarth.

De nombreux chercheurs en sciences de la Terre ont été frustrés par la difficulté de cartographier la sortie du modèle d'actualisation rapide (RAP) des centres nationaux de prévision environnementale (NCEP) de la NOAA dans les SIG (systèmes d'information géographique). Voici la bonne nouvelle : nous avons maintenant un moyen de le faire pour le domaine nord-américain de résolution de 13,5 km (Figure 1) en utilisant les bibliothèques standard récemment mises à jour : PROJ v6.x, GDAL v3.x. Pourquoi est-ce important ? Parce que, alors que la recherche en sciences de l'atmosphère commence à converger avec la recherche en sciences de la terre, des obstacles techniques au partage des données existent toujours et entravent l'innovation et la collaboration. Cette solution supprime l'un de ces obstacles. Il donne aux praticiens du SIG un moyen d'accéder à la sortie du modèle du domaine nord-américain du RAP en utilisant leurs outils dans leur environnement informatique familier. — Sam Batzli, CIMSS/SSEC, Université du Wisconsin-Madison

Reconnaissance

Je travaille sur ce problème par intermittence depuis quelques mois. j'ai tendu la main à Kevin Sampson à UCAR parce que je suis tombé sur un code qu'il avait écrit pour essayer de résoudre ce problème. Nous avons échangé des exemples, des tests et des réflexions dans une série d'e-mails, nous rapprochant chaque fois d'une solution. Nous n'aurions pas pu résoudre ce problème sans travailler ensemble. Il en est de même de Eric James, une filiale de la NOAA au Earth System Research Lab (ESRL) impliquée dans l'intégration des prévisions de fumée dans la sortie du RAP (les données brutes sont encore expérimentales et ne sont pas encore diffusées publiquement). Il m'a présenté les progiciels wgrib2 et NCL (NCAR Command Language) et m'a fourni des paramètres de domaine non publiés qui ont confirmé notre solution ultime.

Arrière-plan

Mon entrée dans ce casse-tête est venue d'une demande du JPSS-PGRR (Joint Polar Satellite System Proving Ground and Risk Reduction), Fire and Smoke Initiative pour mettre des prévisions horaires de fumée du domaine complet RAP North American dans notre plate-forme de visualisation “RealEarth. ” “Pas de problème,” je pense. « Ce sera aussi facile qu'il l'était de mettre dans les prévisions de fumée HRRR CONUS que nous faisons maintenant depuis plus d'un an. » Je me suis trompé. Le problème découlait d'une incompatibilité fondamentale dans la description des projections géographiques entre les disciplines universitaires.

La façon dont l'étendue géographique et la forme du domaine nord-américain complet du RAP (également connu sous le nom de « grille E décalée d'Arakawa sur une grille de latitude/longitude pivotée ») est définie dans l'environnement de modélisation est comme une projection cylindrique équidistante (ou Plate Carrée), dans laquelle la Terre (représentée comme une sphère parfaite) a été tournée vers le sud de telle sorte que le pôle Nord et l'ensemble du cercle arctique sont visibles ainsi que des parties de la Sibérie, de l'Europe du Nord, d'Hawaï et même des parties de l'Amérique du Sud sous l'équateur (voir la figure 1). Pour la modélisation météorologique, il y a des avantages distincts à utiliser une grille Arakawa. Comme indiqué dans son article sur Wikipédia, il offre un moyen de représenter et de calculer des quantités physiques orthogonales (en particulier des quantités liées à la vitesse et à la masse) sur des grilles rectangulaires utilisées pour les modèles du système terrestre pour la météorologie et l'océanographie. fondamentalement erroné de représenter la terre comme une grille décalée sur une sphère en rotation, il n'est pas conventionnel dans le monde de la géodésie, de la géographie et des SIG. Les géographes en sciences de la terre construisaient généralement des cartes comme celle-ci avec une projection azimutale équidistante. L'approche du pôle tourné (appelée par les géographes une transformation oblique générale), fonctionne bien dans un environnement de modélisation, mais est rarement utilisée par les géographes car elle dépend d'un modèle sphérique de la Terre (pas l'ellipsoïde préféré et plus précis, WGS84) et implique deux transformations (une projection cylindrique et une rotation). Il n'était tout simplement pas, et n'est toujours pas, entièrement pris en charge dans la plupart des logiciels SIG. De plus, l'ensemble complet des paramètres nécessaires au SIG pour décrire la "grille E décalée d'Arakawa sur une grille de latitude/longitude tournée" avec une résolution de 13,5 km, sur l'Amérique du Nord, n'a pas été publié en un seul endroit (jusqu'à présent). Passez aux étapes 4, 5 et 6 si vous souhaitez accéder directement aux réponses.

Progiciels et outils plus couramment utilisés dans les sciences de l'atmosphère tels que MATLAB, le module GEOGRID du système de prétraitement des prévisions météorologiques (WRF) (WPS), les opérateurs de données climatiques (CDO), les outils de cartographie génériques (GMT), wgrib2 et NCL offrent un support partiel ou interne pour les projections de poteaux en rotation. À l'exception peut-être de MATLAB, ceux-ci ne sont pas couramment utilisés dans les environnements SIG. L'utilisation de formats de regroupement de données hiérarchiques et multidimensionnelles super efficaces tels que GRIB2, HDF et netCDF comme formats de sortie de données du RAP et de modèles similaires complique encore le problème. Alors que le logiciel SIG s'améliore dans l'analyse des métadonnées pour ces formats, il fait un vide lorsqu'il rencontre des projections et des descriptions non standard de la terre, comme un pôle en rotation. Donc, si vous souhaitez superposer les données de la grille nord-américaine de RAP sur les données du même emplacement géographique dans un SIG sans installer et apprendre un tout nouvel ensemble d'outils, ce n'est pas si facile.

Heureusement, la combinaison logicielle de transformation de PROJ et GDAL peut désormais s'adapter à la notion de "pôle tourné". Malheureusement, des étapes supplémentaires sont nécessaires pour traduire les paramètres du pôle tourné RAP North American Domain en termes compréhensibles et attendus par PROJ ou GDAL. Le reste de cet article décrit les étapes que nous avons suivies pour le faire fonctionner.

Décrire la projection en termes conviviaux pour les SIG

Notre objectif ici est de produire un GeoTIFF avec une projection bien définie que nous pouvons ensuite reprojeter dans quelque chose que les logiciels de cartographie SIG comme QGIS et ArcGIS (qui sont actuellement construits sur les anciennes versions de GDAL et PROJ) peuvent gérer. Ainsi, même lorsque nous définissons correctement le GeoTIFF de sortie (ce que nous faisons à l'étape 5 ci-dessous), il ne fonctionnera pas dans le SIG tant que les futures versions du logiciel n'intégreront pas les dernières fonctionnalités de GDAL et PROJ. Nous appliquons donc une projection supplémentaire à l'étape 6 pour rendre le GeoTIFF rétrocompatible.

Étape 1 : Obtenez des exemples de données

Nous avons d'abord besoin d'un ensemble de données d'exemple avec lequel travailler. Les données du modèle RAP accessibles au public sont disponibles ici : ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/rap/prod/rap.<yyyymmdd>/ (remplacez la date d'aujourd'hui’s par <yyyymmdd> ou accédez à un répertoire de votre choix). Je recommande de vous connecter avec “user:anonymous” et “password:<your email>” en utilisant un client FTP comme FileZilla ou un utilitaire de ligne de commande comme wget ou curl pour récupérer un exemple de jeu de données. Étant donné que les exemples de données que nous voulons sont au format binaire GRIB2, assurez-vous de télécharger avec le type de transfert défini sur “binary.”. Si vous naviguez dans FileZilla, vous remarquerez de nombreux fichiers dans le répertoire “prod”. La plupart sont associés à des grilles AWIPS spécifiques (le système logiciel fermé utilisé par le National Weather Service). Ces fichiers sont des sous-ensembles du domaine nord-américain complet, généralement dans une projection conique conforme de Lambert standard qui fonctionnera déjà dans un SIG car leurs paramètres de projection sont reconnaissables par un SIG, mais ce que nous voulons, c'est un fichier d'exemple du domaine RAP complet avec le modèle “rap.t06z.wrfprsf02.grib2”. Traduire ce nom de fichier… “rap” est le nom de l'implémentation du modèle “t06z” représente le pas de temps horaire en sortie du modèle exécuté en UTC/Zulu (il y a 21 pas de temps d'une heure dans chaque l'exécution du modèle “wrfprs” correspond à la sortie du modèle WRF avec les niveaux de pression, comprenant essentiellement un tableau vertical de niveaux d'altitude, du domaine nord-américain pour la plupart des variables “f02” correspond à l'heure en UTC de l'exécution de la prévision (le modèle est exécuté toutes les heures) et “grib2” est le format de sortie.

Je n'ai testé cette méthodologie que sur les fichiers publics ci-dessus et les fichiers expérimentaux GRIB2 non encore publics qui contiennent des variables supplémentaires telles que la fumée de surface, la fumée intégrée à la colonne. En théorie, cette approche fonctionnera pour tous les fichiers de sortie GRIB2 pour la "grille électronique décalée d'Arakawa sur une grille de latitude/longitude tournée" d'une résolution de 13,5 km, domaine nord-américain complet avec fichiers de sortie GRIB2.

Étape 2 : E Pluribus Unum – Lecture des métadonnées et extraction d'un fichier TIFF

Tout d'abord, voyons ce que nous avons. L'exécution de gdalinfo sur le fichier nouvellement téléchargé nous donnera un résumé. Assurez-vous que vous utilisez la dernière version (ou 3.0 ou supérieure) des outils GDAL. Les versions antérieures ne pourront pas lire ces fichiers GRIB2.

Maintenant, lancez gdalinfo avec le drapeau -proj4 (proj est resté bloqué en v4 pendant si longtemps qu'il est devenu une partie du nom pendant un certain temps). Le haut de la sortie ressemblera à ceci :

Pensez donc à cela comme un cube de données avec plusieurs couches de données représentant un pas de temps de sortie de prévision : quatre dimensions (latitude, longitude, altitude et temps). C'est la beauté de GRIB2. C'est un moyen très efficace d'organiser, de stocker et de déplacer beaucoup d'informations. Comme d'autres formats de données scientifiques hiérarchiques tels que HDF et NetCDF, il s'agit essentiellement d'une base de données de fichiers autonome et auto-descriptive. La difficulté avec GRIB2 dans le cas de “wrfprs” et de tout le domaine nord-américain, c'est que jusqu'à présent, vous aviez besoin d'un logiciel spécialisé commun uniquement dans les sciences de l'atmosphère pour l'ouvrir. Même ainsi, de mon point de vue de novice, il faut une véritable ceinture noire dans le logiciel wgrib2 ou NCL pour obtenir quoi que ce soit de significatif. Mon objectif est de permettre aux spécialistes des sciences de la Terre d'utiliser des outils qui leur sont familiers et qui font partie de leur flux de travail pour accéder à ces importants ensembles de données.

À partir de la commande “gdalinfo” ci-dessus, nous avons maintenant une définition de la sphère utilisée dans cette version du modèle : “+proj=longlat +a=6371229 +b=6371229 +no_defs” (“ +a” est le rayon équatorial et “+b” est le rayon polaire — puisqu'il s'agit d'une sphère, c'est la même chose +R=6371229, “+no_defs” indique à PROJ de ne pas définir valeurs par défaut pour les variables non définies). Nous pouvons également voir que le “Size est 953, 834”. Il s'agit des dimensions en pixels des couches raster en largeur et en hauteur. A part ça, gdalinfo ne nous donne vraiment rien qui nous permette de situer cela dans l'espace géographique. Par exemple, gdalinfo ne reconnaît pas l'existence de coordonnées de coin et ne peut donc pas calculer la taille des pixels (résolution spatiale) dans ce fichier GRIB2 de pôles tournés. Ainsi, ces données ne peuvent tout simplement pas être cartographiées par GDAL/PROJ et les logiciels SIG telles quelles. Nous devrons nous tourner vers les outils wgrib2 ou NCL pour nous aider à collecter plus de paramètres à partir du fichier GRIB2.

Fait amusant : l'UNITE[“degree”,0.0174532925199433]] est la définition standard d'un degré. Il est calculé en divisant Π (pi) par 180 radians.

Si vous avez exécuté la commande gdalinfo ci-dessus, vous aurez sans doute remarqué qu'il y a 779 “Bands” ou ensembles de données dans le fichier GRIB2, et vous avez dû faire défiler les 11 730 lignes pour atteindre le haut de la sortie juste pour voir les informations générales ci-dessus. C'est beaucoup d'informations. Simplifions. Sur beaucoup bandes, extrayons une en TIFF pour voir à quoi ça ressemble. En utilisant une suggestion de Kevin Sampson, nous allons extraire la variable “TSOIL” (température du sol) car contrairement à d'autres variables comme les nuages ​​ou le vent, TSOIL montre suffisamment de caractéristiques géographiques fixes et d'éléments côtiers pour que nous puissions voir avec quoi nous travaillons. alors que nous essayons de capturer ces données dans la réalité géographique. À partir de gdalinfo ci-dessus, nous constatons que TSOIL est situé en tant que bande 628. Maintenant, extrayons un TIFF (je lui donne le nom “unnav” car il est “un-navigated” qui est un terme atmosphérique les scientifiques et les météorologues ont l'habitude de dire qu'il n'a pas de système de coordonnées géographiques ou d'informations de définition, qu'il n'est pas navigué ou qu'il est perdu. Les géographes diraient qu'il n'est pas géolocalisé.

Nous obtenons un TIFF à virgule flottante 64 bits à bande unique. Après avoir joué avec le contraste et l'échelle dans Photoshop et converti en un PNG compatible avec le Web, cela ressemble à ceci (Figure 2) :

Figure 2 : Sortie brute du RAP North American Domain. Notez que le nord et le sud sont inversés.

La première chose que nous remarquons est peut-être qu'il s'agit d'une image miroir de ce à quoi nous nous attendions pour le domaine nord-américain, comme le montre la figure 1. L'emplacement du Groenland et les faibles contours de la côte nord-américaine montrent que les données ont un -haut” plutôt que l'orientation “nord-haut”. Ce n'est pas conventionnel, mais pas « faux » en soi. Pas de jugement. On vient de le noter et de revenir sur le blog, avec une douce gentillesse. Nous pouvons trouver de la joie et de la gratitude dans les dimensions et la forme générales. Ceux-ci ont l'air bien. Appuyons sur.

Étape 3 : extraire l'eau d'une pierre – Obtenez des paramètres supplémentaires

Nous avons besoin de plus d'informations. Pour définir cette projection en termes PROJ, nous utiliserons une transformation oblique générale d'une projection cylindrique équidistante basée sur une sphère. La commande de transformation que nous voulons utiliser ressemble à ceci : >proj +proj=ob_tran +o_proj=eqc +o_lon_p=? +o_lat_p=? +lon_0=? +R=6371229. Par conséquent, pour cette partie, nous avons besoin de connaître l'emplacement du nouveau pôle (+o_lon_p=? +o_lat_p=?) et la longitude centrale de la carte résultante (+lon_0=?). Pour créer notre raster en sortie, nous aurons également besoin de la latitude du centre de projection et des points d'angle. Ainsi, pour une définition complète de cette grille au format SIG-friendly, nous avons besoin d'un minimum de dix paramètres (dont une seconde version de la longitude du centre de projection dont nous parlerons plus loin). A partir de ceux-ci, nous pourrons faire un GeoTIFF compatible SIG :

  • Type de transformation : +proj=ob_tran
  • Type de projection initiale. +o_proj=eqc
  • Rayon de la sphère de projection +R=6371229
  • Dimensions en pixels du raster : -outsize 953 834
  • Nouvelle longitude du pôle : +o_lon_p=?
  • Nouvelle latitude du pôle : +o_lat_p=?
  • Longitude du centre de projection : +lon_0=? (la source)
  • Longitude du centre de projection : +lon_0=? (cibler)
  • Latitude du centre de projection : +lat_0=?
  • Points d'angle (en haut à gauche, en bas à droite) en unités projetées (mètres) : -a_ullr ulx? uly? lrx ? lry?

Après de nombreuses recherches sur Internet et de nombreuses consultations avec Eric James, j'ai obtenu des points de vue qui semblaient plausibles. Mais je voulais les confirmer et avoir un moyen plus fiable d'obtenir ces informations. Sur la suggestion d'Eric, je me suis tourné vers NCL et un outil appelé “ncl_filedump” (avec l'option -c pour les informations de coordonnées), pour vider plus de métadonnées sur ce domaine afin d'obtenir les quatre paramètres restants qui décrivent le modèle de grille de pôles tourné . En faisant cela, j'obtiens ce qui suit :

Nous avons maintenant des points d'angle officiels en longitude-latitude. Comme paires longitude/latitude, nous avons quatre points : 1) -139.0858 -10.5906, 2) -72.91415 -10.5906, 3) 22.66102 46.59194, 4) 125.339 46.59194.

Étape 4: Miracles et pensée magique – Logique et chance – Essais et erreurs

Avec une projection oblique, il est difficile de savoir quel point va où à moins de les tracer comme sur la figure 3. Ici nous utilisons un azimutal équidistant (défini par +proj=aeqd +lat_0=54 +lon_0=-106 +x_0=0 +y_0=0 +a=6371229 +b=6371229 +units=m +no_defs) comme une approximation du domaine nord-américain du pôle tourné, juste pour tester les points d'angle et de centre . Voici ce que nous voyons dans QGIS : LL) -139.0858 -10.5906, LR) -72.91415 -10.5906, UR) 22.66102 46.59194, UL) 125.339 46.59194 Center) -106 54

Figure 3. La cartographie des points d'angle et du centre dans une projection approximative (azimutale équidistante) permet d'identifier quel point est lequel, en termes de supérieur, inférieur, gauche et droit.

Pour préserver notre santé mentale alors que nous plongeons dans la réalité abstraite, il est utile d'utiliser la pensée magique et de considérer ces paramètres en termes de « condition source », ou de la façon dont les données sont mappées nativement (Figure 2) et « condition cible », ” ou comment il apparaît finalement sur la carte projetée (Figure 1). Les points ont fière allure sur la carte de la figure 3, mais nous savons que nos données sont inversées, nous devons donc ici décrire la condition source en termes cibles et inverser les latitudes ou les valeurs y (uniquement). Le coin supérieur gauche devient le coin inférieur gauche et le coin supérieur droit devient le coin inférieur droit comme ceci : UL) -139.0858 -10.5906, LR) 22.66102 46.59194

La longitude centrale est donnée comme 254 et utilise un modèle à 360 degrés du monde plutôt que le plus conventionnel -180 à 180 autour du premier méridien ou “Greenwich.” La soustraction de 360 ​​à 254 nous donne le plus familier -106 (ouest) , qui, avec 54 (nord) a du sens en tant que point central de notre condition cible lorsque nous regardons la figure 3. Pour que le centre de notre projection soit à 54 degrés nord plutôt qu'à l'équateur (ce serait la latitude centrale de la projection cylindrique non tournée), nous devons faire tourner la Terre à l'intérieur de ce cylindre vers le sud de 36 degrés (90-36=54). Poursuivant notre réflexion magique, nous soupçonnons que le modèle de latitude utilise 180 plutôt que le plus conventionnel -90 à 90 et dans notre condition de source inversée, nous reconnaissons que la latitude du pôle nord dans ce modèle est 0, l'équateur est 90 et le pôle sud est de 180. Puisque la latitude du pôle est tournée de 36 degrés, soustrayons 36 de 180, en tournant de 36 degrés au nord du pôle sud (180-36=144). Puisque nous ne voulons pas faire pivoter le pôle est ou ouest et ne savons pas quoi en faire d'autre, nous mettons notre foi en notre puissance supérieure et laissons la source "Nouvelle longitude du pôle" à un 180 générique, avec une volonté de le changer pour un autre générique comme 0 ou 360 si nécessaire. Ce n'était pas nécessaire. C'est le seul paramètre mystérieux que je ne peux pas expliquer. 180 est la seule chose qui semble faire l'affaire dans les terres à l'envers. Peut-être que le pôle sud d'origine était défini comme 180, 180 et puisque nous ne faisons qu'une rotation oblique, cette longitude arbitraire (puisque toute longitude est valable pour le pôle d'une sphère) reste la même. J'accepte les explications suggérées. De ma correspondance avec Eric James, je peux confirmer que c'est bien la valeur attribuée à ce paramètre dans le fichier de démarrage d'entrée du modèle qui lance RAP. [Alerte spoiler pour la longitude de la source du centre de projection : c'est 74 ! J'ai essayé d'utiliser Center Longitude comme -106, mais cela n'a pas fonctionné. J'ai essayé d'utiliser 254 et cela n'a pas fonctionné non plus. J'ai remarqué que dans l'un de nos échanges d'e-mails, mon collègue Keven Sampson essayait d'utiliser 74 et a supposé que peut-être le modèle de longitude comptait jusqu'à 360 dans l'autre sens (…puisque les données sont inversées, cela a du sens : 254-180=74). Mais si je n'avais pas remarqué ce que Kevin a essayé, je doute que j'aurais pensé à cela et cet article de blog n'existerait pas.]

Nous avons donc maintenant tous nos dix paramètres, ensemble au même endroit :

  • Type de transformation : +proj=ob_tran (source)
  • Type de projection initiale. +o_proj=eqc (source)
  • Rayon de la sphère de projection +R=6371229 (source)
  • Dimensions en pixels du raster : -outsize 953 834 (cible)
  • Nouvelle longitude du pôle : +o_lon_p=180 (source)
  • Nouvelle latitude du pôle : +o_lat_p=144 (source)
  • Longitude du centre de projection : +lon_0=74 (source)
  • Longitude du centre de projection : +lon_0=-106 (cible)
  • Latitude du centre de projection : +lat_0=54 (cible)
  • Points d'angle (en haut à gauche, en bas à droite) : -a_ullr -6448701.88 -5642620.27 6448707.45 5642616.94 (source)

“Mais attendez, ” vous demandez, “comment avez-vous obtenu les points d'angle de lon/lat en mètres x/y ?” Bonne question ! Nous avons transformé les points en utilisant PROJ avec les autres paramètres de projection… et un peu d'amorçage …et quelques conjectures …et un peu plus de pensée magique. Dans ses propres mots, “PROJ est un logiciel de transformation de coordonnées générique qui transforme les coordonnées géospatiales d'un système de référence de coordonnées (CRS) à un autre. Cela inclut les projections cartographiques ainsi que les transformations géodésiques. PROJ comprend des applications en ligne de commande pour une conversion facile des coordonnées à partir de fichiers texte ou directement à partir d'une entrée utilisateur.” C'est en partie ce que GDAL, et le fonctionnement interne de nombreux systèmes SIG, utilisent pour manipuler les rasters, pixel par pixel. Nous pouvons utiliser l'application en ligne de commande pour nos points d'angle individuels.

J'ai également découvert que je pouvais l'utiliser comme moyen de tester et de régler plusieurs de nos paramètres de transformation décrits ci-dessus. Dans le processus de recherche d'une solution pour le point central, j'ai essentiellement procédé à une ingénierie inverse de la logique qui a produit les paramètres de latitude et de longitude du pôle et du centre ci-dessus. J'ai testé diverses combinaisons de paramètres en utilisant la transformation oblique générale. Je crois comprendre qu'une projection cylindrique équidistante mappe généralement un globe centré sur lon/lat 0,0 sur un carré bidimensionnel en mètres avec une origine cartésienne de 0,0 et s'étendant vers l'extérieur avec +/-nords sur le l'axe des y et +/- abscisse sur l'axe des x. J'ai donc supposé qu'un point central correct pour cette projection pivotée serait également 0,0. Cela s'est avéré être une supposition chance-logique importante. Après de nombreux essais et erreurs, voici les commandes et les résultats :

C'est vraiment génial que GDAL 3.x reconnaisse et traite la transformation oblique générale décrite par PROJ. Je pense qu'au fil du temps, cela aura un impact significatif sur l'interopérabilité des données et la collaboration entre les sciences de la Terre et de l'atmosphère.

Étape 5 : extraire un GeoTIFF du GRIB2 et appliquer la définition de projection

Donc, tout cela se résume en une seule commande “gdal_translate”. Je devrais peut-être mettre ça en haut.

Étape 6 : Reprojeter pour une compatibilité descendante avec le SIG

Cette étape n'est nécessaire que jusqu'à ce que le logiciel SIG intègre GDAL 3.x et PROJ 6.x. J'ai choisi une projection azimutale équidistante centrée sur -106, 54 comme sur la figure 3 afin que nous puissions admirer une projection similaire dans QGIS.

Le résultat final ressemble à ceci (Figures 4 et 5) et montre un excellent ajustement à la Terre Naturelle 1:10m “Admin 0 – Pays.” Il est temps de lâcher les ballons !


Aide géoTIFF

Geotiff_Information :
Version 1
Révision_clé : 1.0
Tagged_Information :
ModèleTiepointTag (2,3) :
0 0 0
-87.1435563 1.47207446 0
ModelPixelScaleTag (1,3) :
8.97580174e-006 8.97580174e-006 0
End_Of_Tags.
Informations_clés :
GTModelTypeGeoKey (Court,1) : ModelTypeGeographic
GTRasterTypeGeoKey (Court,1) : RasterPixelIsArea
GeographicTypeGeoKey (Court,1) : GCS_WGS_84
GeogCitationGeoKey (Ascii,7) : "WGS 84"
GeogAngularUnitsGeoKey (Court,1) : Angular_Degree
End_Of_Keys.
End_Of_Geotiff.

GCS : 4326/WGS 84
Référence : 6326/Système géodésique mondial 1984
Ellipsoïde : 7030/WGS 84 (6378137.00,6356752.31)
Premier méridien : 8901/Greenwich (0.000000/0d 0'0.00"E)

Coordonnées du coin :
En haut à gauche ( 87d 8'36.80"W, 1d28'19.47"N)
En bas à gauche ( 87d 8'36.80"W, 1d22'2.05"N)
En haut à droite ( 87d 3'25.34"W, 1d28'19.47"N)
En bas à droite ( 87d 3'25.34"W, 1d22'2.05"N)
Centre ( 87d 6' 1,07"W, 1d25'10,76"N)

ce n'est pas du tout là où il doit être converti en wgs84 et rééchantillonné. Il devrait être autour de 34N et 93W. Qu'ai-je fait de mal ?


FSXA Exportation GeoTIFF lat/long

Vous avez peut-être raison. J'ai fait comme suggéré, et resample compile actuellement les cellules (6437 d'entre elles). C'était certainement beaucoup plus simple que je ne le pensais.

Néanmoins, je suis toujours un peu perplexe quant à savoir pourquoi la reprojection dans Global Mapper semble déformée, pourquoi SBX ne reconnaît pas les données dans le fichier .prj et pourquoi SBX pense que je n'ai sélectionné aucun élément à exporter.

Je vais vérifier ces réponses en supposant que votre méthode beaucoup plus simple fait le travail. Merci de m'avoir aidé.

Cela vaut pour tous ceux qui ont eu la gentillesse de répondre ici aussi : Merci !

Kack911

Il a bien compilé et il est activé dans FSX, mais. C'est tout rouge. On dirait Mars. Une idée de quoi il s'agit ?

Scott967

Dans GM, vous devez l'enregistrer en 24 bits, je pense qu'il est par défaut en palette 8 bits, avez-vous vérifié cela?

Hcornée

À partir de la mémoire, la compilation de photos de paysages dans SBX n'est implémentée que pour les fichiers BMP.

J'ai tendance à l'utiliser comme arrière-plan pour les données vectorielles uniquement et à compiler les données manuellement dans un fichier inf.

Les métadonnées sont supprimées par Photoshop et autres de toute façon, c'est donc important pour mon flux de travail.

Je n'ai pas eu à convertir Datum, car toutes mes images ont été en WGS84 . et en tant que tel, n'a jamais eu de difficulté avec SBX à lire les métadonnées des fichiers exportés par Globalmapper.

Je pense qu'il est normal qu'il ne compile pas de photos pour vous.

Saisissez les informations et coupez/collez-les dans un fichier inf.

Kack911

Je me souviens avoir vu cette option, et je pense que je l'ai laissée à son réglage par défaut. Je vais revenir en arrière et essayer l'option 24 bits. Merci pour la suggestion.

La documentation SBX fait plusieurs références au format geoTIFF, mais je suppose que cela pourrait être un héritage. Pour être tout à fait honnête, la méthode .inf utilisant Resample semble de toute façon plus simple. J'utiliserai probablement encore SBX pour faire le trafic et quelques autres fonctionnalités vectorielles, cependant.

Kack911

Je voulais juste vous faire savoir que l'option d'exportation 24 bits dans Global Mapper a résolu le problème du canal RVB. Tout est dans la sim et est joli.

Scott967

Bon à entendre. GM peut également vous fournir un DEM qui fonctionne tout aussi facilement. Exportez simplement en tant qu'élévation geotiff 16 bits.

Kack911

J'ai un nouveau problème, mais je ne pensais pas que je devrais commencer un tout nouveau post.

Resample digère mon GeoTiff, et j'ai une belle résolution de 1 pied pour mon aéroport dans FSX. Mais.

Mes pistes et voies de circulation créées par 3ds Max ne s'alignent pas correctement. Je sais que la réponse évidente serait que j'ai simplement fait quelque chose de mal dans Max, mais j'ai vérifié plusieurs fois et je ne trouve aucune erreur.

J'ai créé les polys au sol en utilisant les mêmes images que j'ai finalement « rééchantillonnées » dans FSX comme arrière-plan. Dans Max 8, cela se superpose parfaitement, mais dans FSX, les choses ne s'alignent pas très bien du tout. Si je répare une extrémité de piste, l'autre extrémité est décalée de près de 10 pieds. Aucune rotation de mon modèle ne va le réparer, semble-t-il.

C'est comme si l'image (dans FSX) était légèrement étirée, donc les lignes parallèles ne sont plus parallèles. C'est très subtil, mais les extrémités opposées d'une piste de 11 000 pieds sont à près de 10 pieds.

J'ai la nette impression que cela est lié à la reprojection que Global Mapper a effectuée sur mes images. J'ai parcouru le SDK et essayé de jouer avec les paramètres xDim et yDim, mais cela n'a pas semblé faire de différence.

Comment puis-je battre ça? Changer les longueurs de mes pistes Max 8 à des valeurs autres que les longueurs du monde réel ne me semble pas être une solution, ce qui me laisse le soin de corriger ou de modifier les images.

Dois-je retirer les informations de géoréférencement du TIF et spécifier manuellement mes propres coordonnées dans le .inf ? A quoi ressemblerait un tel .inf ?

J'espérais utiliser la fonction de quadrillage de Global Mapper pour essayer d'éliminer certaines des images supplémentaires sur les bords de ma zone d'exportation. Sinon, il me restera beaucoup d'espace blanc après avoir fait mon masque alpha. Sans oublier, je veux réduire l'exclusion autogène.

Mais si je vais jouer avec les données de coin du "main" Tiff, comment pourrais-je gérer toutes les tuiles que je veux créer ?

Désolé pour tant de questions, et j'apprécie toutes les idées que vous pourriez avoir.

Scott967

Je doute un peu que le fait de retirer les données georef du Geotiff va changer quoi que ce soit, à moins qu'une erreur d'arrondi ne soit introduite. Il existe un programme gratuit, l'examinateur d'en-tête geotiff geotife.exe qui peut extraire les données d'un geotiff, mais il est tout aussi facile de créer un fichier tfw dans GM et d'obtenir les mêmes données.

AFAIK GM utilise des algorithmes de reprojection "standards de l'industrie". Une chose que je voudrais vérifier que vos paramètres de projection / référence sont corrects (par exemple, vérifiez l'unité "US survey foot") bien que je ne pense pas que cela provoquerait une erreur de 10 pieds sur la longueur d'une piste.

Une chose dans GM est qu'il existe une option "créer des pixels carrés" à l'exportation qui est activée par défaut. Je décoche toujours ça.

J'ai obtenu de bons résultats en peignant mon image dans un programme de peinture avec du violet dans la zone que je ne veux pas, puis en définissant le violet comme valeur nulle dans le fichier .inf. Mon programme de peinture détruit les balises geotiff, mais je viens de renvoyer l'image via GM pour la remettre en place.

Je n'ai jamais essayé de superposer des polys au sol générés par GMAX, donc je ne sais pas quels problèmes de positionnement il pourrait y avoir. Je ne sais pas comment vous déterminez la "vérité sur le terrain" lorsque vous travaillez avec une précision aussi élevée.

Rumbaflappy

Administrateur

Je ne sais pas comment vous placez les images de piste, mais sachez que les caps dans FS sont magnétiques. pas de vrais titres. Vous devez utiliser un programme comme TCalcX pour obtenir un véritable cap dans la simulation.

De plus, l'image rééchantillonnée, maintenant en projection géographique, n'est plus une représentation de la longueur ou de la largeur réelles.

Les projections géographiques sont en degrés de latitude et de longitude. Max ou GMax utilise des mesures linéaires de mètres ou de pieds. Donc, vous ne pouvez vraiment pas utiliser une image de ce type comme arrière-plan max.

Vous pourrez peut-être utiliser une image d'arrière-plan de cette projection UTM, qui est une mesure linéaire (je pense). Peut-être que quelqu'un d'autre peut intervenir ici.

Kack911

@Dick, j'ai oublié de mentionner que tout ce que j'ai conçu dans Max correspond exactement à la documentation du plan directeur de l'aéroport. Les relèvements vrais, les dimensions de la piste, les largeurs des voies de circulation, les distances d'axe à axe et les rayons de congé sont tous spécifiés dans la documentation, et j'ai créé mes polys au sol fidèlement à ces dimensions.

Pour sauvegarder ces données, j'ai importé l'imagerie aérienne que j'ai mentionnée au début de ce fil (toujours en projection State Plane et à 1 pied/px), en la plaçant sur un plan et en l'utilisant comme toile de fond lors de la conception de la chaussée. Le résultat correspondait parfaitement aux dimensions publiées.

In a nutshell that's why I think the reprojected imagery in FSX is the culprit.

@Scott, I've re-exported numerous times with various settings in order to try and isolate any setting that I might have set incorrectly. The "square pixels" option was one of them. I also check the Feet (Survey) setting, and double checked the arc degrees per pixel, all to no avail. Nothing I changed seemed to make any difference in FSX.

I also went through the SDK and fiddled with the INF file. I tried the xDim and yDim entries, but that didn't have any effect either.

Thats what led me to believe that if I manually tell resample to slightly move the Southeast corner of my image, perhaps I can get it to line up properly.

When comparing the imagery to the surrounding vector features in FSX, certain road line up exactly, while others seem to be off slightly. on the order of about 10 feet in places. It's not really about rotation, more like a slight distortion. I just can't tell if it's in the x or the y dimension.

I'd love to post some pics to make what I'm seeing more clear, but I can't really do that, unfortunately.

If you guys have any more ideas, I'm all ears. I'll keep fooling around with it in the mean time. Thanks again.

GaryGB

Somebody please correct me if I'm wrong (as I am not a GMAX user).

Imagery to be used as source data in any file format for FS ground textures / background images / maps / vector content etc. in QUELCONQUE FS utility, whether FS SDK Resample, SHP2VEC, AFCAD, ADE, AFX /Airport Studio, SBuilder, SBuilderX, Slartibartfast. and/or GMAX, 3dsMAX, Sketchup etc. should premier be projected in the required Geographic Lat-Lon WGS84 datum.

The same is true of imagery intended to be used for texturing large areas ex: long RWYs or other terrain regions such as custom ground poly in any of the above FS utilities. just because you can "warp" or "stretch"imagery in an app doesn't mean that will produce the best result compared to an image already projected closer to the final "shape" that FS will require.


FYI: This does not apply to 3D scenery objects such as buildings. vegetation, SimObjects, Aircraft etc. to be placed "above" ground.


Yes imagery will look "squished", and squares will look like rectangles with longer dimensions along the horizontal (East-West) direction. that is what WGS84 devrait "look like" but FS will 'warp' it back into "normal" looking 'square' quad matrix cell tiles at run time, so no worries !


Fais ne pas do any significant manual "stretching" of imagery intended for use in apps with such a 'rubber sheet' feature option (ex: AFX /Airport Studio).

Likewise, do ne pas do any processing of imagery with 'arbitrary' parameters via GDALWarp or other such tools to "warp" an image into what you pense "looks right". il doit be projected correctly as WGS84. or IMHO you are better off NOT attempting to use it.

Administrator

For use in GMax/FSDS/SketchUp the imagery should not be in lat/long projection, but in a local flat earth projection. Since that is the grid those programs work in.

Most other programs will indeed require WGS84 lat/long.

GaryGB

But I thought Christopher Columbus discovered that the Earth n'était pas flat ? <. just kidding, of course ! >


Apparently Sketchup performs the required re-projection when one geo-locates a (max 1024x1024 pixel) image tile directly from Google Earth via its propre internal engine.


Can that same internal engine in Sketchup be used to import ex: a GeoTiff (treated as a generic Tiff I would assume), or a BMP file ?


I suspect that internal engine in Sketchup ne peut pas be used by the end user to do this with any other image import process.


So, IIUC, if we import non-georectified imagery into ex: Sketchup to use as a background for tracing objects to be placed into FS on the ground (ex: custom ground polys), what is a freeware application that end users could utilize to "re-project" or 'warp' the imagery correctement before import ?


I think we'd all like to be certain about this for FS Development, especially considering such issues as I've seen addressed in the Global Mapper forum:

And regarding the making of large FS airport backgrounds / very longue RWYs etc, this was rather intriguing in its implications for would-be 3D modelers:

" We need to decide: what is the acceptable difference between distance measured
on the ground and distance measured on the map.
Local ‘flat earth’ grids only work over a very small area. If your work area extends
beyond a kilometre, you can no longer use a localised flat earth grid . However, you
could define your own local map projection with a small (negligible, but acceptable)
scale-factor. It is relatively easy to convert data between your local projection and the
national grid.
If you are designing a long linear feature, it will usually have its own reference
system – chainage (i.e. distance from a known starting point), along the feature and
offset perpendicular to the feature. You could introduce a variable scale-factor, so that
chainage corresponds with what you measure on the ground. An even neater solution
would be to use what is commonly known as ‘snake’ projection: this dynamically
converts between chainage and grid co-ordinates on large, generally linear projects.
"

" in Great Britain, Ordnance Survey (ordnancesurvey.co.uk) provides
geodata referenced to a map projection known as British National Grid.
In most places, there will be a difference between a distance measured
par terre. This difference is known as scale-error or scale-factor
distortion. It is variable depending upon location in the country, and can affect
measurements by an amount ranging between zero and 4cms every 100m. "

Thanks for your further clarification of what the proper name of that conversion / re-projection format is which we should look for in GIS apps, for subsequent 3D world modeling in GMAX / 3DSMAX and Sketchup with FS as the ultimate destination of such content.


PS: I see that Google Maps (AFAIK the same as Google Earth) says this about their projection and coordinate system:

There are several coordinate systems that the Google Maps API uses:

* Latitude and Longitude values which reference a point on the world uniquely. (Google uses the World Geodetic System WGS84 standard.)
* World coordinates which reference a point on the map uniquely
* Tile coordinates which reference a specific tile on the map at the specific zoom level

World Coordinates

Whenever the Maps API needs to translate a location in the world to a location on a map (the screen), it needs to first translate latitude and longitude values into a "world" coordinate. This translation is accomplished using a map projection. Google Maps uses the Mercator projection for this purpose. You may also define your own projection implementing the google.maps.Projection interface. (Note that interfaces in V3 are not classes you "subclass" but instead are simply specifications for classes you define yourself.)

For convenience in the calculation of pixel coordinates (see below) we assume a map at zoom level 0 is a single tile of the base tile size. We then define world coordinates relative to pixel coordinates at zoom level 0, using the projection to convert latitudes & longitudes to pixel positions on this base tile. This world coordinate is a floating point value measured from the origin of the map's projection to the specific location. Note that since this value is a floating point value, it may be much more precise than the current resolution of the map image being shown. A world coordinate is independent of the current zoom level, in other words.

World coordinates in Google Maps are measured from the Mercator projection's origin (the northwest corner of the map at 180 degrees longitude and approximately 85 degrees latitude) and increase in the x direction towards the east (right) and increase in the y direction towards the south (down). Because the basic Mercator Google Maps tile is 256 x 256 pixels, the usable world coordinate space is <0-256>, "


Hmmm. remarkably similar to the FS world quad matrix cell system of ex: LOD-13 tiles having 256 x256 Area Points each cell !


But then, I've seen rumors posted that a number of folks with the original KeyHole / Google Earth team used to work for ACES.

This was a rather interesting "MS" Document as well:

Introduction to Spatial Coordinate Systems: Flat Maps for a Round Planet


So, I guess we should also look at what MS 'Bing' Maps (aka MS Virtual Earth") uses too:

To make the map seamless, and to ensure that aerial images from different sources line up properly, we have to use a single projection for the entire world. We chose to use the Mercator projection, which looks like this:
Bb259689.150afcdc-99eb-4296-9948-19c0a65727a3(en-us,MSDN.10).jpg

Although the Mercator projection significantly distorts scale and area (particularly near the poles), it has two important properties that outweigh the scale distortion:

1. It’s a conformal projection, which means that it preserves the shape of relatively small objects. This is especially important when showing aerial imagery, because we want to avoid distorting the shape of buildings. Square buildings should appear square, not rectangular.

2. It’s a cylindrical projection, which means that north and south are always straight up and down, and west and east are always straight left and right.

Since the Mercator projection goes to infinity at the poles, it doesn’t actually show the entire world. Using a square aspect ratio for the map, the maximum latitude shown is approximately 85.05 degrees.

To simplify the calculations, we use the spherical form of this projection, not the ellipsoidal form. Since the projection is used only for map display, and not for displaying numeric coordinates, we don’t need the extra precision of an ellipsoidal projection. The spherical projection causes approximately 0.33% scale distortion in the Y direction, which is not visually noticeable. "


2 Answers 2

Your code tries to convert from EPSG:4326 to EPSG:4326.
But your data is actually in EPSG:3857 (wild guess on my end).
You need to use EPSG:3857 as "source" CRS.

GeoJSON had no specification until https://tools.ietf.org/html/rfc7946 so older GeoJSON might not be in WGS84 which rasterio might assume to be the case.

In the end, reprojecting the geoTIFF using QGIS (instead of trying to unsuccessfully reproject geoJSON) to WSG84 did the trick and made the rasterio example work.


Conversion and georeferencing of maps for ttMaps

Before using a map with ttMaps, make sure you have the right! Some publishers of paper maps prohibit you to scan them, and vendors of maps on CDROM impose very restrictive licenses, such as the ban on the use of maps with other software than that provided by the vendor.

Maps format

The format that was chosen for ttMaps is the ECW format, developed by the company Er Mapper, recently acquired by ERDAS.. Why use such a choice? This format has many advantages for mapping applications:

  • Excellent lossless compression ratio, better then that of JPEG
  • Rapid access to image file, at any zoom level
  • Ability to decompress a portion of the image without decompressing the entire file

It is a format chosen by the IGN to deliver maps and aerial photographs. It is therefore possible to directly use the files from IGN. See, for example, screenshots some of which were produced with samples of IGN maps.

Calibrating maps

Before you start converting the maps to ECW format, you must make sure they are calibrated (georeferenced), therefore you need to know:

  • the datum
  • the projection
  • the geographic extents, ie :
    • either the coordinates of the upper left corner and the pixels sizes
    • either the coordinates of the upper left corner and of the lower right corner

    In general, the datum and projection are shown in the legend text of paper maps. If you do not know the geographical extents, you must first calibrate the map. I will not explain in detail how to do this, because there are many programs available, such as:

    Tools for conversion to ECW format

    • Er Mapper software (for Microsoft Windows), of which you can download a trial version (30 day trial) at http://www.erdas.com/Products/ERDASProductInformation/tabid/84/currentid/1052/default.aspx .
    • GDAL free software, available for several systems :
      • Windows : Install the GDAL package. Binaries are available on the GisInternals web site.
      • Linux : GDAL is available in most distributions (choose a package with ECW support)
      • Mac OS X : On the page http://www.kyngchaos.com/software:frameworks, download the frameworks GDAL, UnixImageIO, GEOS 3, SQLite3 and PROJ 4.6, together with the ECW plugin. Install them. In a console, type: export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH

      Conversion of maps with GDAL

      Here are three examples of using GDAL to convert the maps to the ECW format.

      Georeferenced maps

      Some file formats, such as GeoTIFF, include georeferencing data. For example, download the Sample GeoTiff file. After unpacking the archive, you can use the gdalinfo command to display information about the file:

      This information allows us to know the reference system (here WGS84) and the projection used (here, UTM Zone 10 North). Unfortunately, Er Mapper uses its own names to designate the reference systems and projections. Before making the conversion, you must find these names.

      • For France, there is a paper (in French): Syst mes de projections dans Er Mapper et Mapinfo.
      • For the rest of the world, you can read this document: Supported map projections and datums
      • If you installed the Er Mapper software, you can also browse the files in the GDT_DATA folder, especially the files datum.dat and project.dat .

      For the image above, the datum is WGS84 and the projection is NUTM10. The command to use to convert the file to ECW format is:

      The parameter TARGET=90 indicates that you wish to obtain a compression ratio of the image equal to 90% (Change the rate depending on the type of image).

      Note : Avoid using the same base name for input file and output file, otherwise you may encounter some problems with the current version of gdal_translate.

      Conversion of maps shipped with a world file

      Some files, such as TIFF, JPG or PNG may be accompanied by a world file containing the geographical extents of the map. These files have respective extensions .TFW, .JGW and .PGW .

      One example is the SCAN1000 map of IGN, which is freely available on the IGN site.

      Download for example SCAN 1000 Corsica (Version Lambert II extended) and unzip the archive. The useful files are COR_0000.TIF (TIFF image) and COR_0000.TFW, which contains the coordinates and size of pixels. The gdalinfo command shows that this file uses a color palette:

      The gdal_translate utility is not capable of converting a file based on a color palette to a file in 16 million of colors. Fortunately, a conversion utility ( pct2rgb.py ) comes with the GDAL package. You must therefore begin doing this conversion:

      The website of the IGN tells us that this map uses the NTF reference system (New Triangulation French), with the Lambert projection II extended. The conversion command will be as follows:

      However, we note that the map is surrounded by an unnecessary white border. It is possible to convert only a part of the map, using the -srcwin option of the gdal_translate command:

      For the Linux users (or Windows with CygWin), here's a script that converts the SCAN100 map.

      Conversion of a map without georeferencing data

      Download for example the image BlueMarble. The coordinates of the upper left corner are -180 , 90 and those of the lower right corner are 180 ,-90 . The datum is WGS84, and the projection is GEODETIC.

      The command you can use for the conversion is:

      Georeferencing of an ECW file

      If you have a file in ECW format not containing the whole georeferencing information (datum, projection and extents), it is possible to add this information without having to decompress / recompress the file with gdal_translate. To do this, you may use the tool ECW Header Editor , that you can download at ERDAS.

      To georeference the image, you must proceed as follows:

      • Install the software "ECW Header Editor"
      • Launch the software "ECW Header Editor"
      • Open the ECW file through File menu - Open
      • Fill in the Datum and Projection
      • Give the coordinates of the upper left corner of the image and size of pixels
      • Save the image through the File menu - Save.

      It is also possible to use gdal_edit , that works on Windows, MacOS and Linux.

      Georeferencing a file already calibrated for OziExplorer

      There are lots of maps for OziExplorer. Thus, we are going to show step by step, on three practical examples, how to convert your maps to ECW format.

      How to check if the calibration file is correct

      Often it occurs that calibration files for Oziexplorer do not contain the correct information on the projection and the datum. This is because some people use OziExplorer without sufficient knowledge, sometimes without even taking the time to try to understand the datums and projections.

      There is a simple test to be carried out on a calibration file. The principle is to verify that the edges of the map are parallel to the axes of the projection.

      To start, display the .map file contents (it contains plain text). In general, the file contains the coordinates of four control points, most often the four corners of the map. To ensure that the map is aligned with the axes of the projection, you just have to check that:

      • The left corner have the same eastings
      • The right corners have the same eastings
      • The upper corners have the same northings
      • The bottom corners have the same northings

      Example 1

      Point Easting Northing
      Coordinates of the upper left corner 0 2700000
      Coordinates of the upper right corner 1100000 2700000
      Coordinates of the lower right corner 1100000 1600000
      Coordinates of the lower left corner 0 1600000

      In this first example, we see that the map is well positioned. Indeed:

      • The easting of the left edge is 0
      • The easting of the right edge is 11000000
      • The northing of the upper edge is 2700000
      • The northing of the lower edge is 1600000

      Example 2

      Point Easting Northing
      Coordinates of the upper left corner 4.269040 45.882684
      Coordinates of the upper right corner 4.912734 45.870023
      Coordinates of the lower right corner 4.891744 45.422152
      Coordinates of the lower left corner 4.253713 45.434813

      In this second example, we see that the map is misaligned. Indeed:

      • The left edge of the map passes through points of eastings 4.269040 and 4.253713: it is tilted
      • The right egde of the map passes through points of eastings 4.912734 and 4.891744: it is tilted
      • The upper edge of the map passes through points of northings 45.882684 and 45.870023: it is tilted
      • The lower edge of the map passes through points of northings 45.422152 and 45.434813: it is tilted

      Note : This simple method is able to prove that the projection indicated in the .map file is incorrect. It does not check that the projection, datum and coordinates are correct, this should be done by displaying known points in ttMaps.

      Tutorial n 1 : Conversion of a french map coming with a correct calibration file

      Here is an example, with a .map file found on a forum:

      .map file lines La description
      Map Datum,NTF France Datum
      Map Projection,(II) France Zone II,PolyCal,No,AutoCalOnly,No,BSBUseWPX,No Projection
      Image Width,44000 Image width, in pixels
      Image Height,44000 Image height, en pixels
      Point01,xy, 0, 0,in, deg, , ,N, , ,E, grid, , 0, 2700000,N Coordinates of the upper left corner
      Point02,xy,44000, 0,in, deg, , ,N, , ,E, grid, , 1100000, 2700000,N Coordinates of the upper right corner
      Point03,xy,44000,44000,in, deg, , ,N, , ,E, grid, , 1100000, 1600000,N Coordinates of the lower right corner
      Point04,xy, 0,44000,in, deg, , ,N, , ,E, grid, , 0, 1600000,N Coordinates of the lower left corner

      This table contains all the information necessary to create a georeferenced ECW file.

      • The first line gives the datum: New Triangulation French (ErMapper code is NTF)
      • The second line gives the projection: Lambert France Zone II (Ermapper code is LM2FRANC)
      • The fifth line gives the coordinates of the upper left corner: (0, 2700000)
      • With the coordinates of the lower right corner (1100000, 1600000) and the number of pixels, you can deduce the size of pixels:
        • Horizontal (1100000 - 0) / 44000 = 25 m
        • Vertical: (2700000 - 1600000) / 44000 = 25 m

        The easiest way to make the conversion is to use GDAL:

        Note : The LARGE_OK=YES option here is essential, because the size of the map (44000 x 44000 x 3 = 5808000000 bytes) exceeds 500 Mb.

        Tutorial n 2 : Conversion of a french map coming with an incorrect calibration file

        .map file lines La description
        WGS 84,WGS 84, 0.0000, 0.0000,WGS 84 Datum
        Map Projection,Latitude/Longitude,PolyCal,No,AutoCalOnly,No,BSBUseWPX,Yes Projection
        Image Width,5000 Image width, in pixels
        Image Height,5000 Image height, en pixels
        Point01,xy, 0, 0,in, deg, 45, 52.96104,N, 4, 16.1424,E, grid, , , ,N Coordinates of the upper left corner (in degrees + decimale minutes)
        Point02,xy, 5000, 0,in, deg, 45, 52.20138,N, 4, 54.76404,E, grid, , , ,N Coordinates of the upper right corner (in degrees + decimale minutes)
        Point03,xy, 5000, 5000,in, deg, 45, 25.32912,N, 4, 53.50464,E, grid, , , ,N Coordinates of the lower right corner (in degrees + decimale minutes)
        Point04,xy, 0, 5000,in, deg, 45, 26.08878,N, 4, 15.22278,E, grid, , , ,N Coordinates of the lower left corner (in degrees + decimale minutes)

        We see immediately that the map is not aligned with the WGS84 projection. The main challenge will be to find the correct projection and datum.

        As it is a French map, it's not very difficult: Most french maps use the NTF datum and the extented Lambert II projection (at the exception of some maps of Corsica that use the Lambert IV projection).

        We'll check that by converting the coordinates of the four corners of the map into NTF/extented Lambert II. There are many softwares to convert coordinates, use whatever you prefer.

        Points Longitude (WGS84) Latitude (WGS84) X (Lambert IIe) Y (Lambert IIe)
        Upper left corner 4.269040 45.882684 750000.06 2099889.00
        Upper right corner 4.912734 45.870023 799993.98 2099915.04
        Lower right corner 4.891744 45.422152 799982.76 2050105.46
        Lower left corner 4.253713 45.434813 750024.37 2050092.08

        This map has probably been calibrated manually, and there are errors of a few tens of meters. After rounding results, we get:

        Points X (Lambert IIe) Y (Lambert IIe)
        Upper left corner 750000 2100000
        Upper right corner 800000 2100000
        Lower right corner 800000 2050000
        Lower left corner 750000 2050000

        So we see that the map is well aligned with respect to the axes of the projection NTF/extended Lambert II, and we can now convert it:

        Note : The LARGE_OK=YES option is not useful, because the size of the map (5000 x 5000 x 3 = 75000000 bytes) does not exceed 500 MB

        Tutorial n 3 : Conversion of an english map coming with a correct calibration file

        To understand this example, it is necessary to know the reference system used in the United Kingdom. Read for example http://en.wikipedia.org/wiki/British_national_grid_reference_system.

        .map file lines La description
        Ord Srvy Grt Britn,WGS 84, 0.0000, 0.0000,WGS 84 Datum
        Map Projection,(BNG) British National Grid,PolyCal,No,AutoCalOnly,No,BSBUseWPX,No Projection
        IWH,Map Image Width/Height,32704,32576 Taille de l'image, en pixels
        Point01,xy, 200, 128,in, deg, , ,N, , ,W, grid, HV, 01000, 00000,N
        Point02,xy,32400, 128,in, deg, , ,N, , ,W, grid, HW, 62000, 00000,N
        Point03,xy,32400,32328,in, deg, , ,N, , ,W, grid, NG, 62000, 39000,N
        Point04,xy, 200,32328,in, deg, , ,N, , ,W, grid, NF, 01000, 39000,N

        This conversion is more complex than previous ones:

        • First difficulty: the coordinates use a system of letters + numbers, useless with most software.
        • Second difficulty: the coordinates given by the .map file do not correspond to the corners of the map.

        We will begin by converting the coordinates from letters + numbers to fully numerical coordinates. For that, it is enough to know that the squares measure 100km aside and that the reference is the southwest corner of the SV square (see map on Wikipedia). We deduce the following table:

        Carré X Offset Y Offset
        HV 0 1000 km
        HW 100 km 1000 km
        NG 100 km 800 km
        NF 0 km 800 km

        that allows us to calculate the coordinates of the four control points:

        Control Point of reference Easting Northing
        Point01 1000 1000000
        Point02 162000 1000000
        Point03 162000 839000
        Point04 1000 839000

        We note in passing that the map is well aligned with respect to the axes of the projection.

        Now calculate the coordinates of the upper left corner and lower right corner of the map.


        2 Answers 2

        A workaround for this issue is to first use gdalwarp to transform the input image into EPSG:3857 so that no re-projection is done by gdal2tiles . Par example:

        It looks like the image has been reprojected but the areas outside the original image (black pixel areas) have not been assigned as nodata. Also, in your last gdal2tiles run it looks like you're assigning nodata to color white (255).

        Perhaps you could try running gdal2tiles again with the " –srcnodata=0 " option,

        or run gdal_translate on your already-reprojected images to specify black pixels (0) as noData, then re-run gdal2tiles with either the srcnodata option unused or set to black " –srcnodata=0 ".


        How to disable geolocation in browsers?

        Google Chrome

        1. Click the Chrome menu button on the browser toolbar (with the 3 dots).
        2. Click on Settings .
        3. Scroll down and click on Advanced .
        4. In the &lsquoPrivacy and security&rsquo section, click Site settings .
        5. Click &lsquoLocation&rsquo and toggle &lsquoAsk before accessing&rsquo to &lsquoBlocked&rsquo.

        For further information see Google&rsquos location sharing page.

        Firefox

        1. In the URL bar, type about:config .
        2. In the search bar type geo.enabled .
        3. Double click on the geo.enabled preference. Location-Aware Browsing should now be disabled.

        For further information see the Firefox Location-Aware Browsing page.

        Internet Explorer

        1. Open the Tools menu by clicking on the gear icon in the upper-right corner of the browser window.
        2. Open the Privacy tab.
        3. Under Location, select the option Never Allow Websites To Request Your Physical Location .

        Microsoft Edge

        1. Hit the Windows button & select Settings
        2. Navigate to Privacy -> Location and toggle location to Off

        Apple Safari

        1. Choose System Preferences from the Apple () menu.
        2. Click the Security & Privacy icon in the System Preferences window.
        3. Click the Privacy tab.
        4. If the padlock icon in the lower left is locked

          , click it and enter an admin name and password to unlock it
        5. Select Location Services .
        6. Uncheck &lsquoSafari&rsquo to disable geolocation.

        Opj_stream_get_number_byte_left: Assertion `p_stream->m_byte_offset >= 0' failed. #1151

        The text was updated successfully, but these errors were encountered:

        We are unable to convert the task to an issue at this time. Please try again.

        The issue was successfully created but we are unable to update the comment at this time.

        KermMartian commented Sep 25, 2019 •

        Edit: in fact, I replicate this issue with the exact same input file, trying to perform exactly the same conversion. Impressionnant. Makes me slightly doubt the integrity of the file.

        Rouault commented Sep 26, 2019

        The image is fine. This works with GDAL 3.0 (presumably 2.4 as well) with openjpeg >= 2.3.0. But you need a lot of memory (9 GB)

        KermMartian commented Oct 3, 2019

        Using a machine with 256GB of RAM:

        And GDAL 2.4.0 compiled against apt-installed libopenjp2-7-dev (that's OpenJPEG 2.3.0):

        I grabbed a fresh copy of the image to perform this test, then unzipped it, but perhaps that's still tripping me up in some way (e.g., a problem unzipping a 3+GB zip)? Does your sha256sum match?

        Rouault commented Oct 3, 2019

        OK, I finally reproduced the issue. A debug build of openjpeg was needed (at least one without -DNDEBUG), so that assertions are compiled in. The bug was actually in the GDAL OpenJPEG I/O layer, that wasn't ready to deal with reading single-tile codestreams larger than 2GB.

        There was also an error in the MCT stage of openjpeg triggered by opj_decompress, which caused the "Tiles don't all have the same dimension" error. couldn't test the fix since I don't have actually enough RAM to reach that point, but there were a few 'obvious' uint32 overflows I've fixed. Anyway that part wasn't triggered by GDAL since it reads by much smaller chunks, whereas opj_decompress is brute force and decompresses the whole image at once.

        You can’t perform that action at this time.

        You signed in with another tab or window. Reload to refresh your session. You signed out in another tab or window. Reload to refresh your session.


        Help - Information Point on Telecommunications

        The Help tab contains detailed information on using the Information Point on Telecommunications (PIT) and the Electronic Services Platform (PUE).

        The Help tab is arranged by topic and provides information on the method for setting up an account and logging in to the services, as well as instructions, FAQ and instructional videos.

        In the case of questions please contact us using the submission form.

        To access PIT, proceed as follows:

        1. Registration/log-in (the manual is available in Polish version only).

        Registering and logging in to UKE systems takes place via PUE (https://pue.uke.gov.pl) or PIT (https://pit.uke.gov.pl) websites.

        We recommend registration at the PUE level, which enables to fill in and submit “Wniosek o dostęp do PIT” (“Application for access to PIT”).

        2. Setting up an account for an entity (the manual is available in Polish version only).

        To correctly submit “Wniosek o dostęp do PIT” (“Application for access to PIT”) and to continue work in the PIT system, you need to set up an account for an entity (organisation, company) and assign a representative.

        3. Submission of “Wniosek o dostęp do PIT” (“Application for access to PIT” (the manual is available in Polish version only).

        An entity’s representative, acting on behalf of the entity, submits the application. Correct submission of the application is confirmed by an Official Confirmation of Receipt (UPO), which will be delivered to the applicant’s mailbox.

        4. The PIT system administrator immediately grants authorisation based on a submitted application for access to PIT.

        5. Once authorisation is granted pursuant to the application, the user is notified via e-mail.

        6. The user logs in to PIT (https://pit.uke.gov.pl), selects the context of the entity and works in the system.

        Detailed information on registering, logging in, setting up an entity account, changing the operation context and working on a document on the UKE Electronic Services Platform can be found at https://pit.uke.gov.pl/en-us/help/.

        Any questions may be submitted using the submission form.

        Information Point on Telecommunications (PIT)

        As part of the Information Point on Telecommunications (PIT), the Office of Electronic Communications provides access to view services WMS1 and WMTS2 for the different categories of spatial data:

        Service layers after logging in:

        1. Existing linear infrastructure
        2. Planned linear infrastructure
        3. Existing surface infrastructure
        4. Planned surface infrastructure
        5. Existing point infrastructure
        6. Planned point infrastructure
        7. GESUT data

        Service layers without logging in:

        1. Rates for occupying traffic lanes
        2. Forest districts
        3. Access to real estate
        4. SIIS: hotspots
        5. SIIS: co-locations
        6. SIIS: lines
        7. SIIS: planned expansion
        8. SIIS: nodes

        A key is required for layers available after logging in. The key can be generated by authorised entities in the PIT administrative panel.

        1) Web Map Service (WMS) is an international standard for sharing spatial data over the Internet in raster form. Technical standards are available at the Open Geospatial Consortium (OGC) website: http://www.opengeospatial.org/standards/wms.

        2) Web Map Tile Service (WMTS) is an international standard for sharing spatial data over the Internet in raster form, as pre-defined sections of a map, so-called tiles. The process of generating tiles is launched upon the update of a particular product, while the files are appropriately structured and stored on servers. Such a solution speeds up the response of the service to a user’s enquiry concerning a section of a map, since they receive graphics already prepared in advance, as opposed to the WMS service, which generates a graphic file upon each such request.

        For operators of electronic communications networks, in particular new entities on the market, it may be much more efficient to use already existing technical infrastructure, including the one belonging to other public utilities, for the deployment of electronic communications networks, in particular in areas where an adequate electronic communications network is not available or where the construction of a new teletechnical infrastructure may be unprofitable. In addition, sectoral synergies can significantly reduce the need for construction works related to the deployment of electronic communications networks, and thus can also reduce the associated social and environmental costs, such as pollution, nuisance and traffic congestion. For this reason, locating elements of telecommunications networks on a wide and ubiquitous technical infrastructure, such as technical networks used to provide electricity, gas, water supply, sewage and rainwater drainage, heating and transport networks can bring significant savings to investors.

        The basic objectives of running PIT are:

        - to ensure the most efficient planning and deployment of high-speed telecommunications networks,

        - to facilitate the efficient use of technical infrastructure for the purpose of building high-speed telecommunications networks, by providing access to information useful from the point of view of a telecommunications undertaking.

        Entity performing tasks in the field of public utilities - a natural person, a legal person or an organizational unit without legal personality which is granted legal capacity through specific legal provisions, that provides technical infrastructure for the purpose of:

        a) generating, transmitting or distributing gas, electricity or heat,
        b) providing lighting in places referred to in Article 18 (1) item 2 of the Act of 10 April 1997 - Energy Law,
        c) supplying people with water, collecting, transferring, treating or draining sewage, heating, drainage systems, including drainage ducts,
        d) transport, including railways, roads, ports and airports.

        Technical infrastructure - any element of the infrastructure or network that can be used to place in it or on it infrastructure or telecommunications network elements, without becoming an active element of this telecommunications network, such as pipelines, sewers, masts, ducts, chambers, manholes, cabinets, buildings and entrances to buildings, aerial installations, towers and columns, excluding:

        a) cables, including optical fibers,
        b) elements of the network used for the supply of water intended for human consumption,
        c) service ducts within the meaning of Article 4 (15a) of the Act of 21 March 1985 on public roads.

        Service duct – an arrangement of protective housing elements, cable wells and other facilities or devices used for placing or operating:

        a) technical infrastructure devices related to the road management or traffic needs,
        b) telecommunications lines with power supply and power lines, not related to the road management or traffic needs.

        Telecommunications undertaking – an undertaking or another entity authorised to pursue business activities under separate provisions and which conducts business activities consisting in the provision of telecommunications networks, associated facilities or in the provision of telecommunications services, whereby the telecommunications undertaking authorised to provide:

        a) telecommunications services shall be referred to as a “service provider”,
        b) public telecommunications networks or associated facilities shall be referred to as an “operator”.

        Network operator - a telecommunications undertaking or an entity performing public utility tasks, including a local government unit.

        Provisions of the Act of 7 May 2010 on supporting the development of telecommunications services and networks (Mega-law) define groups of information provided to the President of UKE via the PIT portal.

        Table 1. Groups of information to be provided to the President of UKE via the PIT portal.

        On the existing technical infrastructure, other than the infrastructure covered by the inventory, as referred to in Article 29 (1) of the Mega-law, as well as on service ducts, specifying:

        a) their location and arrangement,
        b) type and current status and method of use,
        c) contact details on access matters.

        On investment plans in the scope of the performed or planned construction works, financed in whole or in part from public funds, regarding technical infrastructure or service ducts, specifying:

        a) location and type of works,
        b) element of the technical infrastructure or of the service duct which is subject of the works,
        c) expected date of launching works and their duration,
        d) contact details for coordination of construction works.

        On addresses of the internet websites with the description of the conditions of access, as referred to in Article 30 (1) and (3) of the Mega-law, and of placing objects and devices on the property, as referred to in Article 33 (1) of the Mega-law.

        The following entities are required by the provisions of the Mega-law to provide the President of UKE with the specified information via the PIT portal:

        Table 2. Entities obligated to provide the President of UKE with information via the PIT portal.

        Way of providing
        informations

        Timeframe for providing and updating information

        The PIT portal is constantly connected to the publishing system.

        Time of submission agreed with the President of UKE.

        On the date specified in the application of the President of UKE.

        Information from group 1. in Table 1., at the request of the President of UKE.

        Direct input of information to the PIT system by the entity.

        On the date specified in the application of the President of UKE.

        Not later than within 30 days from the date of completion of the construction of the service duct.

        Within 30 days from the date of issuing the permit.

        Article 40 (8) of the Act of 21 March 1985 on public roads.

        Submission immediately after receiving the application.

        The network operator can provide information to PIT, which will exempt it from the obligation to provide such information at individual requests of telecommunications undertakings. The provided information should be consistent with the actual status and be updated annually by 31 March reflecting the status valid on 31 December of the previous year.

        In accordance with the applicable provisions of law, the following information can be provided:

        1. On the existing technical infrastructure, other than the infrastructure covered by the inventory, as referred to in Article 29 (1) of the Mega-law, as well as on service ducts, specifying:

        a) their location and arrangement,
        b) type and current status and method of use,
        c) contact details on access matters.

        2. On investment plans in the scope of the performed or planned construction works, financed in whole or in part from public funds, regarding technical infrastructure or service ducts, specifying:

        a) location and type of works,
        b) element of the technical infrastructure or of the service duct which is subject of the works,
        c) expected date of launching works and their duration,
        d) contact details for coordination of construction works.


        Voir la vidéo: Pragmatic GDAL (Octobre 2021).