Mais

Por que gdalwarp retorna coordenadas de canto incorretas para meu geotiff?

Por que gdalwarp retorna coordenadas de canto incorretas para meu geotiff?


Eu tenho uma imagem de origem png que está na área igual de Lambert Azimuthal e tenho a extensão longa / lat. Então, executei gdal_translate para criar pontos de empate e converter para gtif. Em seguida, executei gdalwarp para reprojetar para Web Mercator. gdalinfo mostra que as coordenadas delimitadoras estão ok para a saída gdal_translate, mas não para a saída gdalwarp. Que pasa?

Esta é minha operação gdal_translate no arquivo de origem:

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 ">

A razão pela qual você parece estar obtendo valores muito pequenos é porque sua imagem original é projetada, então ela usa unidades lineares - que neste caso PROJ4 eu acho que o padrão é metros - ao invés de graus. Então você está especificando uma parte de aproximadamente 680 metros quadrados do oeste dos EUA.

Agora, como você pode determinar as coordenadas certas para usar é uma questão complicada, porque sua imagem cobre uma área tão grande, você obterá muita distorção nas bordas.

Desculpe, isso não é uma solução - alguém mais se importa em assumir? Se eu tiver um lampejo de inspiração, vou postar aqui.


Figura 1. Exemplo de área de cobertura e extensão espacial do domínio norte-americano para Previsões de Fumaça RAP. Em breve no RealEarth.

Muitos pesquisadores em ciências da Terra ficaram frustrados com a dificuldade de mapear a saída do modelo Rapid Refresh (RAP) dos Centros Nacionais de Predição Ambiental (NCEP) do NOAA e # 8217s em GIS (Sistemas de Informação Geográfica). Aqui está a boa notícia: agora temos uma maneira de fazer isso para o domínio norte-americano com resolução de 13,5 km (Figura 1) usando as bibliotecas padrão recentemente atualizadas: PROJ v6.x, GDAL v3.x. Por que isso importa? Porque à medida que a pesquisa em ciências atmosféricas começa a convergir com a pesquisa em ciências da terra, ainda existem barreiras técnicas para o compartilhamento de dados que inibem a inovação e a colaboração. Esta solução remove uma dessas barreiras. Ele oferece aos praticantes de GIS uma maneira de acessar a saída do modelo do domínio norte-americano do RAP usando suas ferramentas em seu ambiente de computação familiar. & # 8212 Sam Batzli, CIMSS / SSEC, University of Wisconsin-Madison

Reconhecimento

Eu tenho trabalhado nesse problema ocasionalmente por alguns meses. Eu estendi a mão para Kevin sampson na UCAR porque me deparei com um código que ele escreveu para tentar resolver esse problema. Trocamos exemplos, testes e ideias em uma série de e-mails, cada vez mais perto de uma solução. Não poderíamos ter resolvido esse problema sem trabalhar juntos. O mesmo é verdade para Eric James, uma afiliada da NOAA no Earth System Research Lab (ESRL) envolvida em colocar previsões de fumaça na saída do RAP (os dados brutos ainda são experimentais e ainda não foram distribuídos publicamente). Ele me apresentou os pacotes de software wgrib2 e NCL (NCAR Command Language) e forneceu parâmetros de domínio não publicados que confirmaram nossa solução definitiva.

Fundo

Minha entrada neste quebra-cabeça veio através de uma solicitação do JPSS-PGRR (Joint Polar Satellite Proving Ground e Risk Reduction), Fire and Smoke Initiative para colocar previsões de fumaça por hora do domínio total RAP norte-americano em nossa plataforma de visualização & # 8220RealEarth. & # 8221 & # 8220 Sem problema, & # 8221 eu pensei. & # 8220Será tão fácil quanto foi inserir as previsões de fumaça do HRRR CONUS que estamos fazendo há mais de um ano. & # 8221 Eu estava errado. O problema originou-se de uma incompatibilidade básica na descrição de projeções geográficas entre disciplinas acadêmicas.

A forma como a extensão geográfica e a forma do domínio norte-americano completo RAP (também conhecido como uma & # 8220Arakawa escalonada E-grade na grade girada de latitude / longitude & # 8221) é definida no ambiente de modelagem é como uma projeção cilíndrica equidistante (ou Placa Carrée), em que a Terra (representada como uma esfera perfeita) foi girada para o sul de forma que o Pólo Norte e todo o círculo ártico são visíveis, bem como partes da Sibéria, Norte da Europa, Havaí e até mesmo partes da América do Sul abaixo do equador (veja a Figura 1). Para modelagem de clima, existem vantagens distintas em usar uma grade de Arakawa. Conforme observado em sua entrada na Wikipedia, ele oferece uma maneira & # 8220 de representar e calcular quantidades físicas ortogonais (especialmente quantidades relacionadas à velocidade e massa) em grades retangulares usadas para modelos de sistemas terrestres para meteorologia e oceanografia. & # 8221 Embora não haja nada inerentemente errado em representar a Terra como uma grade escalonada em uma esfera girada, não é convencional no mundo da geodésia, geografia e GIS. Geógrafos nas ciências da terra geralmente constroem mapas como este com uma projeção azimutal equidistante. A abordagem de pólo girado (referida pelos geógrafos como uma transformação oblíqua geral), funciona bem em um ambiente de modelagem, mas raramente é usada por geógrafos porque depende de um modelo esférico da Terra (não o elipsóide preferido e mais preciso, WGS84) e envolve duas transformações (uma projeção cilíndrica e uma rotação). Simplesmente não era, e ainda não é, totalmente suportado na maioria dos softwares GIS. Além disso, o conjunto completo de parâmetros necessários ao GIS para descrever a & # 8220Arakawa escalonada E-grid na grade girada de latitude / longitude & # 8221 com resolução de 13,5 km na América do Norte não foi publicado em um só lugar (até agora). Pule para as etapas 4, 5 e 6 se quiser ir direto às respostas.

Pacotes de software e ferramentas mais comumente usados ​​nas ciências atmosféricas, como MATLAB, o módulo GEOGRID dentro do Sistema de Pré-processamento de Previsão de Pesquisa do Tempo (WRF) (WPS), Operadores de Dados Climáticos (CDO), Ferramentas de Mapeamento Genérico (GMT), wgrib2 e NCL oferecem suporte parcial ou interno para projeções de pólo girado. Com exceção talvez do MATLAB, eles não são comumente usados ​​em ambientes GIS. Para complicar ainda mais o problema, está o uso de formatos de agrupamento de dados hierárquicos e multidimensionais supereficientes, como GRIB2, HDF e netCDF, como formatos de saída de dados do RAP e modelos semelhantes. Embora o software GIS esteja cada vez melhor na análise de metadados para esses formatos, ele fica em branco ao encontrar projeções e descrições da Terra fora do padrão, como um pólo girado. Portanto, se você deseja sobrepor dados da grade norte-americana do RAP & # 8217s em dados da mesma localização geográfica em um GIS sem instalar e aprender um novo conjunto de ferramentas, não é tão fácil.

Felizmente, a combinação de software de transformação de PROJ e GDAL agora pode acomodar a noção de um pólo girado & # 8220. & # 8221 Infelizmente, etapas adicionais são necessárias para traduzir os parâmetros do domínio norte-americano RAP do pólo girado em termos compreensíveis e esperados por PROJ ou GDAL. O restante desta postagem descreve as etapas que seguimos para fazer isso funcionar.

Descrevendo a projeção em termos amigáveis ​​ao GIS

Nosso objetivo aqui é produzir um GeoTIFF com uma projeção bem definida que possamos então reprojetar em algo que um software de mapeamento GIS como QGIS e ArcGIS (que atualmente são construídos em versões anteriores de GDAL e PROJ) possa lidar. Portanto, mesmo quando definimos corretamente a saída GeoTIFF (o que fazemos na Etapa 5 abaixo), ela não funcionará no GIS até que as versões futuras do software incorporem os recursos mais recentes do GDAL e do PROJ. Portanto, aplicamos uma projeção adicional na Etapa 6 para tornar o GeoTIFF compatível com versões anteriores.

Etapa 1: Obtenha alguns dados de exemplo

Primeiro, precisamos de um conjunto de dados de exemplo para trabalhar. Os dados do modelo RAP de acesso público estão disponíveis aqui: ftp://ftpprd.ncep.noaa.gov/pub/data/nccf/com/rap/prod/rap.<yyyymmdd>/ (substitua a data de hoje & # 8217s por & ltyyyymmdd & gt ou navegue até um diretório da sua escolha). Recomendo fazer login com & # 8220user: anonymous & # 8221 e & # 8220password: & ltyour email & gt & # 8221 usando um cliente FTP como FileZilla ou utilitário de linha de comando como wget ou curl para obter um conjunto de dados de exemplo. Como os dados de exemplo que desejamos estão no formato binário GRIB2, certifique-se de fazer o download com o tipo de transferência definido como & # 8220binary. & # 8221 Se você estiver navegando no FileZilla, notará muitos arquivos no diretório & # 8220prod & # 8221. A maioria está associada a Grids AWIPS específicos (o sistema de software fechado usado pelo National Weather Service). Esses arquivos são subconjuntos do domínio norte-americano completo, geralmente em uma projeção Lambert Conformal Conic padrão que já funcionará em um GIS porque seus parâmetros de projeção são reconhecíveis por um GIS, mas o que queremos é um arquivo de exemplo do domínio RAP completo com o padrão & # 8220rap.t06z.wrfprsf02.grib2 & # 8221. Traduzindo este nome de arquivo & # 8230 & # 8220rap & # 8221 é o nome da implementação do modelo & # 8220t06z & # 8221 significa a saída da etapa de hora do modelo executado em UTC / Zulu (há 21 etapas de uma hora em cada model run & # 8220wrfprs & # 8221 significa saída do modelo WRF com níveis de pressão, essencialmente compreendendo uma matriz vertical de níveis de altitude, do domínio norte-americano para a maioria das variáveis ​​& # 8220f02 & # 8221 significa a hora em UTC da previsão de execução (o modelo é executado a cada hora) e & # 8220grib2 & # 8221 é o formato de saída.

Eu apenas testei essa metodologia nos arquivos públicos acima e nos arquivos experimentais GRIB2 ainda não públicos que contêm variáveis ​​adicionais, como fumaça de superfície, fumaça integrada de coluna. Em teoria, essa abordagem funcionará para todos os arquivos de saída GRIB2 para a grade E escalonada & # 8220Arakawa na grade girada de latitude / longitude & # 8221 de resolução de 13,5 km, domínio norte-americano completo com arquivos de saída GRIB2.

Etapa 2: E Pluribus Unum & # 8211 Lendo os metadados e extraindo um TIFF

Primeiro, vamos ver o que temos. Executar gdalinfo no arquivo recém-baixado nos dará um resumo. Certifique-se de que está executando a versão mais recente (ou 3.0 ou superior) das ferramentas GDAL. Versões anteriores não conseguirão ler esses arquivos GRIB2.

Agora execute gdalinfo com o sinalizador -proj4 (proj ficou preso na v4 por tanto tempo que se tornou parte do nome por um tempo). A parte superior da saída ficará assim:

Portanto, pense nisso como um cubo de dados com várias camadas de dados representando uma etapa de tempo de saída de previsão: quatro dimensões (latitude, longitude, altitude e tempo). Essa é a beleza do GRIB2. É uma forma altamente eficiente de organizar, armazenar e mover muitas informações. Como outros formatos de dados científicos hierárquicos, como HDF e NetCDF, é essencialmente um banco de dados de arquivos autocontido e autoexplicativo. A dificuldade com GRIB2, no caso de & # 8220wrfprs & # 8221 e do domínio norte-americano completo, é que até agora você precisava de um software especializado comum apenas em ciências atmosféricas para abri-lo. Mesmo assim, da minha perspectiva de novato, é necessário um verdadeiro cinturão negro em wgrib2 ou software NCL para obter algo significativo. Meu objetivo é permitir que as pessoas nas ciências da Terra usem ferramentas conhecidas e que fazem parte de seu fluxo de trabalho para acessar esses conjuntos de dados importantes.

A partir do comando & # 8220gdalinfo & # 8221 acima, agora temos uma definição da esfera que é usada nesta versão do modelo: & # 8220 + proj = longlat + a = 6371229 + b = 6371229 + no_defs & # 8221 (& # 8220 + a & # 8221 é o raio equatorial e & # 8220 + b & # 8221 é o raio polar & # 8212 visto que esta é uma esfera, é o mesmo + R = 6371229, & # 8220 + no_defs & # 8221 instrui PROJ a não definir valores padrão para variáveis ​​indefinidas). Também podemos ver que o & # 8220Size é 953, 834 & # 8221. Estas são as dimensões em pixels das camadas raster em largura e altura. Fora isso, gdalinfo realmente não nos dá nada que nos permita situar isso no espaço geográfico. Por exemplo, gdalinfo não reconhece a existência de coordenadas de canto e, portanto, não pode calcular o tamanho do pixel (resolução espacial) neste arquivo GRIB2 de pólo girado. Portanto, esses dados simplesmente não podem ser mapeados pelos softwares GDAL / PROJ e GIS como estão. Teremos que recorrer às ferramentas wgrib2 ou NCL para obter ajuda na coleta de mais parâmetros do arquivo GRIB2.

Curiosidade: UNIT [& # 8220degree & # 8221,0.0174532925199433]] é a definição padrão de um grau. É calculado dividindo Π (pi) por 180 radianos.

Se você executou o comando gdalinfo acima, não terá dúvidas de ter notado que existem 779 & # 8220Bands & # 8221 ou conjuntos de dados no arquivo GRIB2, e você teve que rolar de volta por todas as 11.730 linhas para chegar ao topo da saída apenas para ver as informações gerais acima. Isso é muita informação. Vamos simplificar. De muitos bandas, vamos & # 8217s extrair 1 como um TIFF para ver sua aparência. Usando uma sugestão de Kevin Sampson, vamos extrair a variável & # 8220TSOIL & # 8221 (temperatura do solo) porque, ao contrário de outras variáveis ​​como nuvens ou vento, TSOIL mostra características geográficas fixas e elementos costeiros suficientes para que possamos ver com o que estamos trabalhando enquanto tentamos encaixar esses dados na realidade geográfica. A partir do gdalinfo acima, descobrimos que TSOIL está localizado como Banda 628. Agora vamos & # 8217s extrair um TIFF (I & # 8217m dando a ele o nome & # 8220unnav & # 8221 porque é & # 8220un-navegado & # 8221 que é um termo atmosférico os cientistas e meteorologistas costumam significar que ele não tem um sistema de coordenadas geográficas ou informações definidoras, ele não é navegado ou está perdido. Os geógrafos diriam que ele não é & # 8220 geografado & # 8221).

Obtemos um TIFF de ponto flutuante de banda única de 64 bits. Depois de brincar com o contraste e a escala no Photoshop e converter para um PNG amigável para a web, fica assim (Figura 2):

Figura 2: Saída bruta do RAP North American Domain. Observe que o norte e o sul estão invertidos.

Talvez a primeira coisa que notamos é que esta é uma imagem espelhada do que esperávamos para o domínio norte-americano, conforme visto na Figura 1. A localização da Groenlândia e os contornos tênues da costa norte-americana mostram que os dados têm uma & # 8220sul -up & # 8221 em vez de orientação & # 8220north-up & # 8221. Não é convencional, mas não & # 8220 errado & # 8221 per se. Sem julgamento. Apenas anotamos e voltamos ao blog, com gentileza. Podemos encontrar alegria e gratidão nas dimensões e formas gerais. Parece bom. Vamos continuar.

Etapa 3: espremendo a água de uma pedra & # 8211 Obtenha parâmetros adicionais

Precisamos de mais informações. Para definir esta projeção em termos de PROJ, usaremos uma transformação oblíqua geral de uma projeção cilíndrica equidistante baseada em esfera. O comando de transformação que queremos usar se parece com este: & gtproj + proj = ob_tran + o_proj = eqc + o_lon_p =? + o_lat_p =? + lon_0 =? + R = 6371229. Portanto, para esta parte, precisamos saber a localização do novo pólo (+ o_lon_p =? + O_lat_p =?) E a longitude central do mapa resultante (+ lon_0 =?). Para fazer nosso raster de saída, também precisaremos da latitude do centro de projeção e dos pontos de canto. Portanto, para uma definição completa desta grade em formato amigável ao GIS, precisamos de um mínimo de dez parâmetros (incluindo uma segunda versão da longitude do centro de projeção, sobre a qual falaremos mais tarde). A partir deles, seremos capazes de fazer um GeoTIFF compatível com GIS:

  • Tipo de transformação: + proj = ob_tran
  • Tipo de projeção inicial. + o_proj = eqc
  • Raio da esfera de projeção + R = 6371229
  • Dimensões em pixels do raster: -outsize 953 834
  • Nova longitude do pólo: + o_lon_p =?
  • Nova latitude do pólo: + o_lat_p =?
  • Longitude do centro de projeção: + lon_0 =? (fonte)
  • Longitude do centro de projeção: + lon_0 =? (alvo)
  • Latitude do centro de projeção: + lat_0 =?
  • Pontos de canto (superior esquerdo, inferior direito) em unidades projetadas (metros): -a_ullr ulx? uly? lrx? lry?

Depois de muitas pesquisas na Internet e muitas consultas com Eric James, obtive pontos críticos que pareciam plausíveis. Mas eu queria confirmá-los e ter uma maneira mais confiável de obter essas informações. Por sugestão de Eric & # 8217s, recorri a NCL e uma ferramenta chamada & # 8220ncl_filedump & # 8221 (com a opção -c para informações de coordenadas), para despejar mais metadados sobre este domínio para obter os quatro parâmetros restantes que descrevem o modelo de grade de pólo girado . Fazendo isso, obtenho o seguinte:

Agora temos pontos de esquina oficiais em longitude-latitude. Como pares de longitude / latitude, temos quatro pontos: 1) -139,0858 -10,5906, 2) -72,91415 -10,5906, 3) 22,66102 46,59194, 4) 125,339 46,59194.

Etapa 4: Milagres e pensamento mágico & # 8211 Lógica e sorte & # 8211 Tentativa e erro

Com uma projeção oblíqua, é difícil saber qual ponto vai para onde, a menos que você os plote como na Figura 3. Aqui estamos usando um equidistante azimutal (definido por + proj = aeqd + lat_0 = 54 + lon_0 = -106 + x_0 = 0 + y_0 = 0 + a = 6371229 + b = 6371229 + unidades = m + no_defs) como uma aproximação do domínio norte-americano do pólo girado, apenas para testar o canto e os pontos centrais. Aqui & # 8217s o que vemos em QGIS: LL) -139.0858 -10.5906, LR) -72.91415 -10.5906, UR) 22.66102 46.59194, UL) 125.339 46.59194 Centro) -106 54

Figura 3. Mapear os cantos e pontos centrais em uma projeção aproximada (Azimutal Equidistante) ajuda a identificar qual ponto é qual, em termos de superior, inferior, esquerdo e direito.

Para preservar nossa sanidade enquanto mergulhamos na realidade abstrata, é útil usar o pensamento mágico e considerar esses parâmetros em termos de & # 8220 condição da fonte & # 8221 ou como os dados são mapeados nativamente (Figura 2) e & # 8220 condição de destino, & # 8221 ou como finalmente aparece no mapa projetado (Figura 1). Os pontos parecem ótimos no mapa da Figura 3, no entanto, sabemos que nossos dados estão invertidos, então aqui precisamos descrever a condição de origem em termos de destino e inverter as latitudes ou valores y (apenas). Superior esquerdo torna-se Inferior Esquerdo, e Superior Direito torna-se Inferior Direito, desta forma: UL) -139,0858 -10,5906, LR) 22.66102 46.59194

A longitude central é dada como 254 e usa um modelo de 360 ​​graus do mundo em vez do mais convencional -180 a 180 em torno do Meridiano Principal ou & # 8220Greenwich. & # 8221 Subtrair 360 de 254 nos dá o mais familiar -106 (oeste) , que junto com 54 (norte) faz sentido como um ponto central de nossa condição-alvo quando olhamos para a Figura 3. Para obter o centro de nossa projeção a 54 graus norte em vez do equador (que seria a latitude central de projeção cilíndrica não girada), temos que girar a Terra dentro desse cilindro para o sul em 36 graus (90-36 = 54). Continuando nosso pensamento mágico, suspeitamos que o modelo de latitude está usando 180 em vez do mais convencional -90 a 90 e em nossa condição de fonte invertida, estamos reconhecendo que a latitude do pólo norte neste modelo é 0, o equador é 90 e o pólo sul é 180. Como a latitude do pólo é girada em 36 graus, vamos & # 8217s subtrair 36 de 180, girando 36 graus ao norte do pólo sul (180-36 = 144). Visto que não queremos girar o pólo leste ou oeste e não sabemos o que mais fazer com ele, colocamos nossa fé em nosso Poder Superior e deixamos a fonte & # 8220Novo pólo longitude & # 8221 em 180 genéricos, com uma vontade de alterá-lo para algum outro genérico como 0 ou 360, se necessário. Não foi necessário. Este é o único parâmetro misterioso que não consigo explicar. 180 é a única coisa que parece funcionar em terreno invertido. Talvez o pólo sul original tenha sido definido como 180, 180 e como estamos apenas fazendo uma rotação oblíqua, essa longitude arbitrária (já que qualquer longitude é válida para o pólo de uma esfera) permanece a mesma. Eu agradeço as explicações sugeridas. Pela minha correspondência com Eric James, posso confirmar que é realmente o valor atribuído a esse parâmetro no arquivo de inicialização de entrada do modelo que inicia o RAP. [Alerta de spoiler para longitude da fonte do centro de projeção: é 74! Tentei usar a longitude central como -106, mas não funcionou. Tentei usar o 254 e também não funcionou. Percebi que em uma de nossas trocas de e-mail meu colega Keven Sampson estava tentando usar 74 e presumiu que talvez o modelo de longitude contasse até 360 na outra direção (& # 8230 visto que os dados são invertidos, isso realmente faz muito sentido: 254-180 = 74). Mas se eu não tivesse percebido o que Kevin tentou, duvido que teria considerado isso e esta postagem do blog não existiria.]

Portanto, agora temos todos os dez dos nossos parâmetros, juntos em um só lugar:

  • Tipo de transformação: + proj = ob_tran (fonte)
  • Tipo de projeção inicial. + o_proj = eqc (fonte)
  • Raio da esfera de projeção + R = 6371229 (fonte)
  • Dimensões em pixels do raster: -outsize 953 834 (alvo)
  • Novo pólo longitude: + o_lon_p = 180 (fonte)
  • Nova latitude do pólo: + o_lat_p = 144 (fonte)
  • Longitude do centro de projeção: + lon_0 = 74 (fonte)
  • Longitude do centro de projeção: + lon_0 = -106 (alvo)
  • Latitude do centro de projeção: + lat_0 = 54 (alvo)
  • Pontos de canto (superior esquerdo, inferior direito): -a_ullr -6448701.88 -5642620.27 6448707.45 5642616.94 (fonte)

& # 8220Mas espere, & # 8221 você pergunta, & # 8220como você conseguiu os pontos de canto de lon / lat em x / y metros? & # 8221 Boa pergunta! Transformamos os pontos usando PROJ com os outros parâmetros de projeção & # 8230 e um pouco de inicialização & # 8230 e algumas suposições & # 8230 e um pouco mais de pensamento mágico. Em suas próprias palavras, & # 8220PROJ é um software de transformação de coordenadas genérico que transforma coordenadas geoespaciais de um sistema de referência de coordenadas (CRS) para outro. Isso inclui projeções cartográficas, bem como transformações geodésicas. PROJ inclui aplicativos de linha de comando para fácil conversão de coordenadas de arquivos de texto ou diretamente da entrada do usuário. & # 8221 É em parte que o GDAL e o funcionamento interno de muitos sistemas GIS usam para manipular rasters, pixel por pixel. Podemos usar o aplicativo de linha de comando para nossos pontos de canto individuais.

Também descobri que poderia usar isso como uma forma de testar e ajustar vários de nossos parâmetros de transformação descritos acima. No processo de encontrar uma solução para o ponto central, eu essencialmente fiz uma engenharia reversa da lógica que produziu os parâmetros de latitude e longitude do pólo e do centro acima. Testei várias combinações de parâmetros usando a transformação oblíqua geral. Meu entendimento é que uma projeção cilíndrica equidistante geralmente mapeia um globo que está centrado em lon / lat 0,0 em um quadrado bidimensional em unidades de metros com uma origem cartesiana de 0,0 e se estendendo para fora com +/- norte no eixo y e +/- leste no eixo x. Então, eu deduzi que um ponto central correto para essa projeção girada também seria 0,0. Isso acabou sendo um importante palpite lógico-feliz. Depois de muita tentativa e erro, aqui estão os comandos e resultados:

É realmente ótimo que GDAL 3.x reconheça e processe a transformação oblíqua geral que PROJ descreve. Acredito que, com o tempo, isso terá um impacto significativo na interoperabilidade e colaboração de dados em toda a Terra e nas ciências atmosféricas.

Etapa 5: Extraia um GeoTIFF do GRIB2 e aplique a definição de projeção

Então, tudo isso é condensado em um comando & # 8220gdal_translate & # 8221. Talvez eu deva colocar isso no topo.

Etapa 6: Reprojetando para compatibilidade com versões anteriores com GIS

Esta etapa é necessária apenas até que o software GIS incorpore GDAL 3.xe PROJ 6.x. I & # 8217m escolhendo uma projeção azimutal equidistante centrada em -106, 54 como na Figura 3 para que possamos admirar uma projeção de aparência semelhante no QGIS.

O resultado final é semelhante a este (Figuras 4 e 5) e mostra um excelente ajuste ao Natural Earth 1: 10m & # 8220Admin 0 - Países. & # 8221 É hora de soltar os balões!


Ajuda geoTIFF

Geotiff_Information:
Versão 1
Key_Revision: 1.0
Tagged_Information:
ModelTiepointTag (2,3):
0 0 0
-87.1435563 1.47207446 0
ModelPixelScaleTag (1,3):
8.97580174e-006 8.97580174e-006 0
End_Of_Tags.
Keyed_Information:
GTModelTypeGeoKey (Short, 1): ModelTypeGeographic
GTRasterTypeGeoKey (Short, 1): RasterPixelIsArea
GeographicTypeGeoKey (curto, 1): GCS_WGS_84
GeogCitationGeoKey (Ascii, 7): & quotWGS 84 & quot
GeogAngularUnitsGeoKey (Short, 1): Angular_Degree
End_Of_Keys.
End_Of_Geotiff.

GCS: 4326 / WGS 84
Datum: 6326 / World Geodetic System 1984
Elipsóide: 7030 / WGS 84 (6378137,00, 6356752,31)
Meridiano principal: 8901 / Greenwich (0,000000 / 0d 0 '0,00 & quotE)

Coordenadas de canto:
Superior esquerdo (87d 8'36,80 & quotW, 1d28'19,47 & quotN)
Inferior esquerdo (87d 8'36,80 & quotW, 1d22 '2,05 & quotN)
Superior direito (87d 3'25,34 & quotW, 1d28'19,47 & quotN)
Inferior direito (87d 3'25,34 & quotW, 1d22 '2,05 & quotN)
Centro (87d 6 '1,07 & quotW, 1d25'10,76 & quotN)

não está onde precisa ser convertido para wgs84 e reamostrado. Deve ser em torno de 34N e 93W. O que eu fiz errado ?


FSXA Exportando GeoTIFF lat / long

Talvez você esteja certo. Fiz como sugerido, e a nova amostra está compilando as células (6437 delas). Certamente era muito mais simples do que eu pensava.

No entanto, ainda estou um pouco perplexo porque a reprojeção no Global Mapper parece distorcida, porque SBX não reconhece os dados no arquivo .prj e porque SBX pensa que eu não selecionei nenhum item para exportação.

Vou dar uma checada nessas respostas, presumindo que seu método muito mais simples faça o trabalho. Obrigado por me ajudar.

Isso vale para todos que tiveram a gentileza de responder aqui também: Obrigado!

Kack911

Compilou ok e está ativado no FSX, mas. Está tudo vermelho. Parece Marte. Alguma ideia do que se trata?

Scott967

No GM você tem que salvar como 24 bits, acho que o padrão é a paleta de 8 bits, você verificou isso?

Hcornea

A partir da memória, a compilação do photoscenery no SBX é implementada apenas para arquivos BMP.

Eu costumo usá-lo como pano de fundo apenas para dados vetoriais e compilar os dados manualmente em um inf.

Os metadados são removidos pelo Photoshop e similares, então isso é importante para o meu fluxo de trabalho.

Não precisei converter Datum, pois todas as minhas imagens estavam em WGS84. e, como tal, nunca teve dificuldade com o SBX lendo os metadados dos arquivos exportados do Globalmapper.

Acho normal que ele não compile o photoscenery para você.

Pegue as informações e recorte / cole em um inf.

Kack911

Lembro-me de ver essa opção e acho que a deixei na configuração padrão. Vou voltar e tentar a opção de 24 bits. Obrigado pela sugestão.

A documentação SBX faz várias referências ao formato geoTIFF, mas suponho que isso possa ser algo legado. Para ser totalmente honesto, o método .inf usando Resample parece mais fácil de qualquer maneira. Provavelmente ainda usarei o SBX para fazer o tráfego e alguns outros recursos de vetor, no entanto.

Kack911

Só queria que você soubesse que a opção de exportação de 24 bits no Global Mapper corrigiu o problema do canal RGB. Tudo está no sim e está lindo.

Scott967

Bom ouvir. O GM também pode obter um DEM que funciona com a mesma facilidade. Basta exportar como elevação geotiff de 16 bits.

Kack911

Tenho um novo problema, mas não achei que deveria começar um novo post.

Resample digere meu GeoTiff, e eu tenho uma bela fotocenografia com resolução de 1 pé para meu aeroporto em FSX. Mas.

Meu 3ds Max criou pistas e pistas de taxiamento não se alinham corretamente. Sei que a resposta óbvia seria que simplesmente fiz algo errado no Max, mas verifiquei várias vezes e não encontro erros.

Eu criei os polis de solo usando as mesmas imagens que eventualmente 'reamostrou' no FSX como meu plano de fundo. No Max 8, ele se sobrepõe perfeitamente, mas no FSX, as coisas não se alinham muito bem. Se eu consertar uma das pontas da pista, a outra sairá quase 3 metros. Nenhuma quantidade de rotação do meu modelo vai consertar isso, ao que parece.

É como se a imagem (no FSX) fosse ligeiramente esticada, então as linhas paralelas não são mais paralelas. É muito sutil, mas as extremidades opostas de uma pista de 11.000 pés estão a quase 3 metros de distância.

Tenho a nítida sensação de que isso está relacionado à reprojeção que o Mapeador Global realizou em minhas imagens. Eu olhei no SDK e tentei mexer nas configurações de xDim e yDim, mas não parecia fazer nenhuma diferença.

Como posso vencer isso? Alterar os comprimentos de minhas pistas do Max 8 para valores diferentes dos comprimentos do mundo real não parece uma solução para mim, então isso me deixa com a correção ou alteração das imagens.

Devo retirar as informações de georreferenciamento do TIF e especificar manualmente minhas próprias coordenadas no .inf? Qual seria a aparência de um .inf?

Eu esperava usar o recurso Gridding do Global Mapper para tentar eliminar algumas das imagens extras em torno das bordas da minha área de exportação. Caso contrário, terei muito espaço em branco sobrando depois de fazer minha máscara alfa. Sem mencionar que quero reduzir a exclusão de autógenos.

Mas se eu mexer nos dados dos cantos do & quotmain & quot Tiff, como gerencio todos os blocos que desejo criar?

Desculpe por tantas perguntas e agradeço qualquer insight que vocês possam ter.

Scott967

Eu meio que duvido que tirar os dados georreferenciados do Geotiff vá mudar alguma coisa, a menos que haja um erro de arredondamento sendo introduzido. Existe um programa freeware, geotiff header examiner geotiffe.exe, que pode extrair os dados de um geotiff, mas é igualmente fácil criar um arquivo tfw no GM e obter os mesmos dados.

AFAIK GM usa algoritmos de reprojeção & quot padrão da indústria & quot. Uma coisa, eu verificaria se seus parâmetros de projeção / datum estão corretos (por exemplo, verifique a unidade & quotUS Pé de pesquisa & quot), embora eu não ache que isso causaria um erro de 10 pés no comprimento de uma pista.

Uma coisa no GM é que há uma opção & quotcriar pixels quadrados & quot na exportação que está habilitada por padrão. Eu sempre desmarco isso.

Obtive bons resultados pintando minha imagem em um programa de pintura com roxo na área que não queria e, em seguida, definindo roxo como meu nullValue no .inf. Meu programa de pintura destrói as tags geotiff, mas eu apenas executo a imagem de volta ao GM para colocá-la de volta.

Eu nunca tentei sobrepor polígonos de solo gerados por GMAX, então não sei quais problemas de posicionamento podem haver. Não sei como você determina a "verdade fundamental" quando está trabalhando com a precisão desejada.

Rhumbaflappy

Administrador

Não sei como você está posicionando as imagens da passarela, mas esteja ciente de que as direções em FS são magnéticas. não são títulos verdadeiros. Você precisa usar um programa como TCalcX para obter um título verdadeiro no sim.

Além disso, a imagem reamostrada, agora em projeção geográfica, não é mais uma representação do comprimento ou largura reais.

As projeções geográficas estão em graus de latitude e longitude. Max ou GMax usa medidas lineares de metros ou pés. Portanto, você realmente não pode usar uma imagem desse tipo como plano de fundo máximo.

Você pode ser capaz de usar uma imagem de fundo dessa projeção UTM, que é uma medida linear (eu acho). Talvez outra pessoa possa entrar aqui.

Kack911

@Dick, esqueci de mencionar que tudo o que projetei no Max está exatamente de acordo com a documentação do Plano Diretor do Aeroporto. Rolamentos reais, dimensões de pista, larguras de pista de taxiamento, distâncias de linha de centro a linha de centro e raios de filete são especificados nos documentos, e criei meus polis de solo fielmente a essas dimensões.

Para fazer backup desses dados, importei as imagens aéreas que mencionei no início deste tópico (ainda em projeção State Plane e a 1 pé / px), colocando-as em um plano e usando-as como pano de fundo ao projetar o pavimento. O resultado foi uma combinação perfeita com as dimensões publicadas.

Em suma, é por isso que acho que as imagens reprojetadas no FSX são as culpadas.

@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 ALGUM FS utility, whether FS SDK Resample, SHP2VEC, AFCAD, ADE, AFX /Airport Studio, SBuilder, SBuilderX, Slartibartfast. and/or GMAX, 3dsMAX, Sketchup etc. should primeiro 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.


Para sua informação: 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 deve "look like" but FS will 'warp' it back into "normal" looking 'square' quad matrix cell tiles at run time, so no worries !


Fazer não 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 não do any processing of imagery with 'arbitrary' parameters via GDALWarp or other such tools to "warp" an image into what you pensar "looks right". isto devo be projected correctly as WGS84. or IMHO you are better off NOT attempting to use it.

Administrador

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ão era flat ? <. just kidding, of course ! & gt


Apparently Sketchup performs the required re-projection when one geo-locates a (max 1024x1024 pixel) image tile directly from Google Earth via its own 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 não pode be used by the end user to do this with any other image import process.


So, IIUC, if we import não-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 devidamente 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 longo RWYs etc, this was rather intriguing in its implications for would-be 3D modelers:

& quot 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 . No entanto, você
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.
& quot

& quot 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
on the ground. 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. & quot

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. & quot


2 respostas 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

      Exemplo 1

      Apontar 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

      Exemplo 2

      Apontar 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 Descrição
      Map Datum,NTF France Datum
      Map Projection,(II) France Zone II,PolyCal,No,AutoCalOnly,No,BSBUseWPX,No Projeção
      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 Descrição
        WGS 84,WGS 84, 0.0000, 0.0000,WGS 84 Datum
        Map Projection,Latitude/Longitude,PolyCal,No,AutoCalOnly,No,BSBUseWPX,Yes Projeção
        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.

        Pontos 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:

        Pontos 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 Descrição
        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 Projeção
        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:

        Quadrado 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 respostas 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 . Por exemplo:

        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.

        Raposa de fogo

        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 borda

        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. Impressive. 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) Rede 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
        em formação

        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.


        Assista o vídeo: QGIS + GDAL: Conversão de GeoTIFF para KMZ com sobreposição de Raster no Google Earth