Actualmente estoy trabajando para importar datos climáticos CANGRID (proporcionados como Surfer Grid ascii, archivos ".grd") en ArcGIS. La cuadrícula tiene 95 filas por 125 columnas de tamaño. Los metadatos proporcionan lat / lon de origen (esquina inferior izquierda), tamaño de celda (50 km), así como proyección de notas como estereografía polar con meridiano central (110 grados W) y latitud de origen (60 grados N).
Después de intentar convertir el .grd a rásteres como .ascii y .flt sin éxito, logré usar GDAL para establecer la extensión y la proyección, sin embargo, el conjunto de datos no se alinea correctamente con los límites del área deseada. Ver la imagen de abajo.
¿Existe una geotransformación aceptada para la estereografía polar que pueda explicar esta falta de alineación?
Por ejemplo, ¿hay un factor de conversión específico o rotación que debería estar usando?
Un archivo de ejemplo del conjunto de datos está aquí: "t201113.grd"
Aquí está el código que estoy usando actualmente en 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
`Además, aquí está la información de proyección, tal como se define por la entrada, es decir, desde "GetProjection ()":
'PROJCS ["North_Pole_Stereographic", GEOGCS ["GCS_WGS_1984", DATUM ["WGS_1984", SPHEROID ["WGS_1984", 6378137.0,298.257223563]], PRIMEM ["Greenwich", 0.0], UNIT ["Grado", 0.0174532, 0.0174532, 0.0174532 PROYECCIÓN ["Stereographic"], PARÁMETRO ["False_Easting", 0.0], PARAMETER ["False_Northing", 0.0], PARAMETER ["Central_Meridian", 0.0], PARAMETER ["Scale_Factor", 1.0], PARAMETER ["Latitude_Of_Origin", 90 ], UNIDAD ["50_Kilómetros", 50000.0]] '
Y la entrada GeoTransform:
(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)Lat, también se proporcionan longitudes de las coordenadas de la cuadrícula, y cuando se ve en el sistema de coordenadas proyectado, se ve a continuación. Cuando la geotransformada se define mediante coordenadas del cordinado inferior izquierdo (amarillo) o superior derecho (rosa), puedo establecer de manera efectiva la extensión, pero persisten problemas de alineación en todo el ráster.
fuente
Respuestas:
Demasiado tiempo para un comentario, esto es para acompañar la respuesta de @ Matej .
Use la función "Ráster a otro formato" y convierta su archivo .grd a un formato ESRII GRID. Esto es importante porque la mayoría de las funciones ráster en ArcGIS son accesibles solo para este formato, ya sea eso o generalmente se vuelve demasiado lento cuando lo usa en otros formatos.
Como ya tiene el archivo de proyección asociado al archivo. En lugar de proyectar los nuevos datos convertidos, defina su proyección. ArcToolbox> Herramientas de administración de datos> Proyecciones y transformaciones> Definir proyección. Puede navegar a la proyección ESRII gráfica estéreo polar predefinida y ver si sus parámetros coinciden con los que figuran en los metadatos (no lo hace), por lo que puede modificarlo según @Matej. Solo aquí: en lugar de modificar, cree uno nuevo basado en la proyección NPS con el meridiano central y la latitud de origen cambiados y guárdelo en el disco, luego navegue a la nueva proyección y úsela cuando defina su proyección. Esto se debe a que su modificación sobre la marcha no estará disponible más adelante cuando desee usarla para establecer el sistema de coordenadas para su marco de datos,
fuente
La información de proyección adecuada para estos archivos se puede encontrar aquí: Definir cadena Proj4 para el conjunto de datos CANGRD
fuente
No creo que necesites reproyectar la imagen. Solo haga lo siguiente:
Tenga en cuenta que la imagen (grd) ya está en la proyección StereoGraphic del Polo Norte, que solo le indica cómo ajustar el mapa base que se alineará con la imagen.
Paso 1 :
Modifique la proyección estereográfica del Polo Norte original (WKID: 102018) para ajustar la latitud de origen y el meridiano central:
Paso 2:
Georreferencia el archivo grd configurando la esquina inferior izquierda a la coordenada especificada (lat, lon). Cuando actualiza la georreferenciación, el archivo .gdwx se crea en la misma carpeta. Al asignar la esquina SW a (40.0451, -129.853), el contenido del archivo tiene el siguiente aspecto:
editar: el archivo del mundo anterior se ha modificado manualmente en función del tamaño de la celda y la ubicación proporcionada de la esquina SW: la quinta y sexta línea representan la ubicación calculada del píxel superior izquierdo de la imagen. La posición de la imagen cambió ligeramente.
Los valores anteriores colocan (desplazan) la imagen a la ubicación especificada y definen la escala.
Y esta es la salida:
Si esto no parece estar alineado correctamente, cuestionaría las coordenadas proporcionadas para la esquina SW de la imagen. En caso de que tenga acceso a las coordenadas, es decir, la esquina NE de la imagen, puede volver a calcular los parámetros de transformación que escalarían y rotarían la imagen entre dos (o más) puntos.
fuente