Atualmente, estou trabalhando para importar dados climáticos do CANGRID (fornecidos como Surfer Grid ascii, arquivos ".grd") para o ArcGIS. A grade é de 95 linhas por 125 colunas de tamanho. Os metadados fornecem lat / lon de origem (canto inferior esquerdo), tamanho da célula (50 km), bem como projeção de notas como estereográfica polar com meridiano central (110 graus W) e latitude de origem (60 graus N).
Após tentar converter o .grd em rasters como .ascii e .flt, sem êxito, consegui usar o GDAL para definir a extensão e a projeção, no entanto, o conjunto de dados não se alinha corretamente aos limites da área pretendida. Veja a imagem abaixo.
Existe uma geotransformação aceita para estereografia polar que possa explicar essa falta de alinhamento?
Por exemplo, existe um fator de conversão ou rotação específico que eu deveria usar?
Um arquivo de exemplo do conjunto de dados está aqui: "t201113.grd"
Aqui está o código que estou usando atualmente no 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
`Além disso, aqui estão as informações da projeção, conforme definidas pela entrada, ou seja, de "GetProjection ()":
'PROJCS ["North_Pole_Stereographic", GEOGCS ["GCS_WGS_1984", DATUM ["WGS_1984", SPHEROID ["WGS_1984", 6378137.0,298.257223563]], PRIMEM ["Greenwich", 0.0], UNIDADE ["3232", 1995], 33] 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 ], UNIDADE ["50_Kilômetros", 50000.0]] '
E a entrada GeoTransform:
(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)Lat, comprimentos das coordenadas da grade também são fornecidos e, quando visualizados no sistema de coordenadas projetado, são parecidos abaixo. Quando a geotransformação é definida por coordenadas do cordado inferior esquerdo (amarelo) ou superior direito (rosa), posso definir efetivamente a extensão, mas ainda há problemas de alinhamento durante a varredura.
fonte
Respostas:
Muito tempo para um comentário, é para acompanhar a resposta de @ Matej .
Use a função “Raster para outro formato” e converta seu arquivo .grd em um formato ESRII GRID. Isso é importante porque a maioria das funções de varredura no ArcGIS são acessíveis apenas para este formato, isso ou geralmente fica muito lento quando você o usa em outros formatos.
Como ele já possui o arquivo de projeção associado ao arquivo. Em vez de projetar os novos dados convertidos, defina sua projeção. ArcToolbox> Ferramentas de Gerenciamento de Dados> Projeções e Transformações> Definir Projeção. Você pode navegar para a projeção ESRII polar estéreo polar predefinida e ver se seus parâmetros correspondem aos dados nos metadados (não), para modificá-lo conforme @Matej. Somente aqui - em vez de modificar, crie um novo com base na projeção NPS com o meridiano central e Latitude de origem alterados e salve-o em disco, navegue até a nova projeção e use-a ao definir sua projeção. Isso ocorre porque sua modificação imediata não estará disponível mais tarde quando você desejar usá-la para definir o sistema de coordenadas para seu quadro de dados,
fonte
As informações de projeção adequadas para esses arquivos podem ser encontradas aqui: Define a seqüência de caracteres do Proj4 para o conjunto de dados CANGRD
fonte
Eu não acho que você precise reprojetar a imagem. Faça apenas o seguinte:
Observe que a imagem (GRD) já está na projeção estéreo do Polo Norte, que fornece apenas a indicação de como ajustar o mapa base que será alinhado com a imagem.
Etapa 1 :
Modifique a projeção estereográfica do Polo Norte original (WKID: 102018) para ajustar a latitude de origem e o meridiano central:
Passo 2:
Georeferencie o arquivo grd configurando o canto inferior esquerdo para a coordenada especificada (lat, lon). Quando você atualiza o georreferenciamento, o arquivo .gdwx é criado na mesma pasta. Ao atribuir o canto SW a (40.0451, -129.853), o conteúdo do arquivo é o seguinte:
editar: o arquivo mundial acima foi modificado manualmente com base no tamanho da célula e no local fornecido do canto SW - a 5ª e a 6ª linhas representam o local calculado do pixel superior esquerdo da imagem. A posição da imagem mudou ligeiramente.
Os valores acima colocam (deslocam) a imagem no local especificado e definem a escala.
E esta é a saída:
Se isso não parecer alinhado corretamente, eu questionaria as coordenadas fornecidas para o canto SW da imagem. Caso você tenha acesso às coordenadas do canto NE da imagem, é possível recalcular os parâmetros de transformação que seriam dimensionados e girariam a imagem entre dois (ou mais) pontos.
fonte