Géotransformation pour la stéréographie polaire?

16

Je travaille actuellement pour importer des données climatiques CANGRID (fournies sous forme de fichiers Surfer Grid ascii, ".grd") dans ArcGIS. La grille est de 95 lignes par 125 colonnes. Les métadonnées fournissent le lat / lon d'origine (coin inférieur gauche), la taille des cellules (50 km) ainsi que la projection des notes en stéréographie polaire avec méridien central (110 degrés W) et latitude d'origine (60 degrés N).

Après avoir tenté de convertir le .grd en rasters en .ascii et .flt sans succès, j'ai réussi à utiliser GDAL pour définir l'étendue et la projection, mais l'ensemble de données ne s'aligne pas correctement avec les limites de la zone prévue. Voir l'image ci-dessous.

Existe-t-il une géotransformation acceptée pour la stéréographie polaire qui pourrait expliquer ce manque d'alignement?

Par exemple, y a-t-il un facteur de conversion ou une rotation spécifique que je devrais utiliser?

Un exemple de fichier de l'ensemble de données est ici: "t201113.grd"

Voici le code que j'utilise actuellement dans GDAL

ds = gdal.Open("t201113.grd")
array = ds.ReadAsArray()

x_rotation = 0
y_rotation = 0
xres = 1
yres = -1

llx = -129.8530
lly = 40.0451
ulx = -175.144
uly = 71.385

input_osr = osr.SpatialReference()
input_osr.ImportFromWkt(ds.GetProjection())

wgs84_osr = osr.SpatialReference()
wgs84_osr.ImportFromEPSG(4326)

wgs_to_nps_trans = osr.CoordinateTransformation(wgs84_osr, input_osr)
x, y, z = wgs_to_nps_trans.TransformPoint(llx,lly)

geo_transform = [ x, xres, x_rotation, y, y_rotation, yres ]

ncol = ds.RasterXSize
nrow = ds.RasterYSize

out_driver = gdal.GetDriverByName("HFA")
out_ds = out_driver.Create(t201113.img", ncol, nrow, 1, gdal.GDT_Float32)

out_ds.SetGeoTransform(geo_transform)

out_prj = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-110.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["50_Kilometers",50000.0]]'

out_ds.SetProjection(out_prj)

out_ds.GetRasterBand(1).WriteArray(array)
out_ds.GetRasterBand(1).SetNoDataValue(1.70141e+038)
out_ds.FlushCache()
out_ds = None
`

Voici également les informations de projection, telles que définies par l'entrée, c'est-à-dire de "GetProjection ()":

'PROJCS ["North_Pole_Stereographic", GEOGCS ["GCS_WGS_1984", DATUM ["WGS_1984", SPHEROID ["WGS_1984", 6378137.0,298.257223563]], PRIMEM ["Greenwich", 0,0], UNIT ["Degree", 0,099432] PROJECTION ["Stereographic"], PARAMETER ["False_Easting", 0.0], PARAMETER ["False_Northing", 0.0], PARAMETER ["Central_Meridian", 0.0], PARAMETER ["Scale_Factor", 1.0], PARAMETER ["Latitude_Of_Origin", 90.0 ], UNIT ["50_Kilomètres", 50000.0]] »

Et l'entrée GeoTransform:

(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)

Lat, les coordonnées longues de la grille sont également fournies, et lorsqu'elles sont affichées dans le système de coordonnées projeté, elles sont comme ci-dessous. Lorsque la géotransformation est définie par les coordonnées de la coordonnée inférieure gauche (jaune) ou supérieure droite (rose), je peux effectivement définir l'étendue, mais il reste des problèmes d'alignement tout au long du raster.

entrez la description de l'image ici

jsnider
la source
Si vous utilisez ArcGIS, passez au pôle Nord stéréographique et définissez le parallèle standard sur 60,0. L'implémentation stéréographique ArcGIS utilise un facteur d'échelle plutôt qu'un parallèle standard car le projet peut être centré n'importe où.
mkennedy
Merci @mkennedy - parlez-vous du projet "North Pole Stereographic" (WKID 102018)? J'ai défini la latitude d'origine et les valeurs du méridien central à l'aide de cette projection et j'ai toujours le même problème. Vous faites peut-être référence à une autre projection?
jsnider
Non, vous en avez besoin où la projection (méthode) est Stereographic_North_Pole. Je ne pense pas que nous ayons le PCS exact; essayez de modifier à partir de 3995 ou 3413.
mkennedy
1
Les métadonnées notent que "Le fichier grid_pnt_lls.txt répertorie les lat / longs pour chaque x / y (0,0 = coin SW de la grille)." Avec ce fichier en main, vous pouvez reprojeter cette grille à n'importe quel système de coordonnées que vous souhaitez et continuer à partir de là.
whuber
1
Où pouvons-nous télécharger la couche vectorielle à tester?
Farid Cheraghi

Réponses:

2

Trop de temps pour un commentaire, c'est pour accompagner la réponse de @ Matej .

  1. Ajoutez les données «.grd» dans ArcGIS.
  2. Utilisez la fonction «Raster vers un autre format» et convertissez votre fichier .grd en un format ESRII GRID. Ceci est important car la plupart des fonctions raster d'ArcGIS sont accessibles uniquement pour ce format, soit cela, soit cela devient généralement trop lent lorsque vous l'utilisez sur d'autres formats.

  3. Comme il a déjà le fichier de projection associé au fichier. Plutôt que de projeter les nouvelles données converties, définissez sa projection. ArcToolbox> Outils de gestion des données> Projections et transformations> Définir la projection. Vous pouvez naviguer vers la projection ESRII graphique stéréo polaire prédéfinie et voir si ses paramètres correspondent à ceux donnés dans les métadonnées (ce n'est pas le cas), vous pouvez donc la modifier selon @Matej. Seulement ici - plutôt que de modifier, créez-en une nouvelle basée sur la projection NPS avec le méridien central et la latitude d'origine modifiés et enregistrez-le sur le disque, puis accédez à la nouvelle projection et utilisez-la lors de la définition de votre projection. En effet, votre modification à la volée ne sera pas disponible ultérieurement lorsque vous souhaitez l'utiliser pour définir le système de coordonnées de votre bloc de données,

yanes
la source
1

Je ne pense pas que vous ayez besoin de reprojeter l'image. Procédez uniquement comme suit:

  1. définir (modifier) ​​la projection du fond de carte
  2. géoréférencer (déplacer) l'image vers l'emplacement spécifié

Notez que l'image (grd) est déjà dans la projection StereoGraphic du pôle Nord, ce qui ne vous donne que des indications sur la façon d'ajuster le fond de carte qui sera aligné avec l'image.

Étape 1 :

Modifiez la projection stéréographique originale du pôle Nord (WKID: 102018) pour ajuster la latitude d'origine et le méridien central:

entrez la description de l'image ici

Étape 2:

Géoréférencer le fichier grd en définissant le coin inférieur gauche sur les coordonnées spécifiées (lat, lon). Lorsque vous mettez à jour le géoréférencement, le fichier .gdwx est créé dans le même dossier. Lors de l'attribution d'un coin SW à (40.0451, -129.853), le contenu du fichier se présente comme suit:

50000
0
0
-50000
-1730620.4315
2744092.9724

modifier: le fichier mondial ci-dessus a été modifié manuellement en fonction de la taille de la cellule et de l'emplacement du coin SW - la 5e et la 6e ligne représentent l'emplacement calculé du pixel supérieur gauche de l'image. La position de l'image a légèrement changé.

Les valeurs ci-dessus placent (décalent) l'image à l'emplacement spécifié et définissent l'échelle.

Et voici la sortie: entrez la description de l'image ici

Si cela ne semble pas être correctement aligné, je remettrais en question les coordonnées fournies pour le coin SW de l'image. Dans le cas où vous avez accès aux coordonnées du coin NE de l'image, vous pouvez recalculer les paramètres de transformation qui redimensionneraient et tourneraient l'image entre deux (ou plus) points.

Matej
la source
Le fichier .gdwx est-il un fichier mondial ? Si c'est le cas, alors l'article wiki lié dit que les lignes 5 et 6 font référence au pixel supérieur gauche. Suggérez-vous de le régler sur le pixel sud-ouest (en bas à gauche)?
Kirk Kuykendall
Non. J'ai seulement spécifié l'emplacement du coin SW comme indiqué dans le fichier Lisez-moi. La structure du fichier ressemble au fichier mondial. Il a été généré à l'aide de l'outil de géoréférencement d'ArcMap qui peut avoir calculé et stocké l'emplacement du coin nord-ouest - n'a pas encore été vérifié.
Matej
1
Oui, je l'ai vérifié maintenant. L'emplacement stocké dans le .gdwx est le coin supérieur gauche.
Matej