Mais cache

No post anterior intitulado cache ficou por referir o projecto MapProxy (erro lastimável da minha parte). Tratando-se também de um projecto da família Python, a provocação final continua a valer.

Sobre o MapProxy ver "TileCache,GeowebCache and MapProxy - a technical and usability comparison", apresentação do FOSS4G 2010.

Fica também a referência ao TileSeeder, interface gráfico o qual pode ajudar bastante na pré-geração de tiles. Ver também a apresentação do FOSS4G 2010, "TILESEEDER; A NEW TILE MANAGEMENT TOOL"

Bom fds

GIS Cloud no Blogger com iFrame

A integração do GIS Cloud num blog (Blogger neste caso) é tão simples como inserir um iFrame no corpo da mensagem. São ainda disponibilizadas outras formas de integração nomeadamente GIS Cloud API (Javascript), Google Maps (Javascript), Google Maps (Adobe Flash) e OpenLayers.

cache

Deixando de parte o TileCache e o GeoWebCache, dois clássicos consagrados no que diz respeito à cache de tiles resultantes de pedidos WMS, fica a referência a alguns projectos os quais jogam, ou tentam jogar, no mesmo campeonato.

JTileCache - Mencionado apenas por razões históricas, implementa a proposta WMS-C (como de resto os restantes projectos referidos) tendo sido escrito integralmente em JAVA. É o antepassado directo do GeoWebCache.

Por esta altura, e mesmo antes, era também frequente tentar correr o TileCache sobre Jython. A este propósito ver Jython + TileCache/FeatureServer: It just Works do Christopher Schmidt ou Adapting TileCache to work as a java servlet do Jon Blower, estando este igualmente por trás do gae-wms.

SpatialCache - Escrito também em Python, pode-se considerar o irmão mais novo do TileCache. Pode ser utilizado como mecanismo de cache para fins genéricos.

PIRI - Este é particularmente interessante. Uma implementação de TileCache para App Engine. Comigo não funciona mas ainda não desisti.

geoCache - Para quem procura uma abordagem .NET (neste caso 3.5) esta é a solução, sendo escrita em C#. Resta dizer que teve por base, em grande parte, a migração de código do TileCache.

Desta pequena enumeração nota-se um certo ascendente da família TileCache/Python. Alguma explicação?

Transformação de coordenadas, programaticamente

Quem queira utilizar programaticamente o conversor de coordenadas disponível aqui, pode sempre fazer algo do género:

Exemplo funcional disponível aqui.

De notar a necessidade de incluir não só o script proj4js.js mas também o tmerc.js.

Para os sistemas de coordenadas suportados ver a listagem do conversor.

A utilização no OpenLayers aplicada a features é possível através da classe OpenLayers.Projection, mas isso é outra história.

Códigos EPSG revisitados

Visto que, desde a última vez que este assunto foi abordado aqui no blog, terem surgido algumas alterações nos códigos EPSG relativos a Portugal, nomeadamente no que diz respeito aos identificadores das regiões insulares (a questão insofismável do ITRF93, importantíssima para quem faz a distinção relativamente ao WGS84) bem como o Datum 73/ Hayford-Gauss (o 27492 foi promovido a 27493) para não falar do ETRS89 (esse lobisomem que se confunde também com o WGS84, e que aparece agora em 4 sabores distintos), volta-se a publicar a lista actualizada à data do presente post:

Códigos EPSG dos sistemas de referência utilizados em Portugal

Portugal Continental - Sistemas Globais

• EPSG: 4936 (ETRS89/ Coordenadas Geocêntricas)
• EPSG: 4937 (ETRS89/ Coordenadas Geográficas 3D)
• EPSG: 4258 (ETRS89/ Coordenadas Geográficas 2D)
• EPSG: 3763 (ETRS89/ PT-TM06)

Portugal Continental - Sistemas Locais

• EPSG: 4274 (Datum 73/ Coordenadas Geográficas 2D)
• EPSG: 27493 (Datum 73/ Hayford-Gauss)
• EPSG: 4207 (Datum Lisboa/ Coordenadas Geográficas 2D)
• EPSG: 5018 (Datum Lisboa/ Hayford-Gauss)
• EPSG: 20790 (Datum Lisboa/ Hayford-Gauss com falsa origem - Coordenadas Militares)

Arquipélagos dos Açores e da Madeira - Sistemas Globais

• EPSG: 5011 (ITRF93/ Coordenadas Geocêntricas)
• EPSG: 5012 (ITRF93/ Coordenadas Geográficas 3D)
• EPSG: 5013 (ITRF93/ Coordenadas Geográficas 2D)
• EPSG: 5014 (ITRF93/ PTRA08 - UTM zona 25N) - Grupo Ocidental do Arquipélago dos Açores
• EPSG: 5015 (ITRF93/ PTRA08 - UTM zona 26N) - Grupo Central e Oriental do Arquipélago dos Açores
• EPSG: 5016 (ITRF93/ PTRA08 - UTM zona 28N) – Madeira, Porto Santo, Desertas e Selvagens

Arquipélagos dos Açores e da Madeira - Sistemas Locais

• EPSG: 2188 (Datum Observatório - Flores (Grupo Ocidental do Arquipélago dos Açores) / UTM zona 25N)
• EPSG: 2189 (Datum Base SW - Graciosa (Grupo Central do Arquipélago dos Açores) / UTM zona 26N
• EPSG: 2190 (Datum S. Braz - S. Miguel (Grupo Oriental do Arquipélago dos Açores) / UTM zona 26N)
• EPSG: 2942 (Datum Base SE - Porto Santo (Madeira) / UTM zona 28N)

Mais uma vez a fonte é o site do IGP (m@pas online) ao qual recorremos amiúde.

cp4 georeferenciados (centróides)

A forma mais comum de código postal, o de 4 dígitos, vulgo cp4, baseia-se na ideia de ter o território coberto por polígonos não sobrepostos sendo a cada um desses polígonos atribuído um identificador único: o cp4.

No caso de Portugal Continental temos 463 polígonos podendo ser calculado para cada um deles o respectivo centróide. Na posse dessa informação é possível publicar um webservice que devolva as coordenadas do centróide para um determinado cp4. Sendo assim,

--> Saber as coordenadas em WGS84 do centróide relativo ao polígono de um determinado código postal de 4 dígitos.

http://codigospostais.appspot.com/cp4? Devolve a latitude e a longitude para um determinado código postal de 4 dígitos.

Exemplo: http://codigospostais.appspot.com/cp4?codigo=2675

Parâmetros: codigo código postal de 4 dígitos

Parâmetro opcional: callback para a chamada de uma função em javascript

Resposta: JSON

Nota: serviço só disponível para códigos de Portugal Continental.

Este e outro serviços em http://geodivagar.appspot.com

Consumindo json web services com Python (exemplos)

Há já bastante tempo que o JSON (JavaScript Object Notation) deixou de ser um exclusivo do javascript. Actualmente, a maioria das linguagens de programação têm suporte a este formato, seja nativamente, seja através de uma biblioteca qualquer. Isto torna-o ideal para intercâmbio de dados.

O Python tem suporte nativo a JSON, na sua Standard Library, desde a versão 2.6. Para versões anteriores será necessário utilizar bibliotecas externas, sendo a mais conhecida a simplejson.

Como os exemplos seguintes foram testados em Python 2.5, também irei utilizar a simplejson a qual terá que estar instalada no sistema. Note-se que o Python 2.5 é a versão suportada pelo Google App Engine, logo quem também quiser desenvolver nesta plataforma e em Python terá que utilizar esta versão, por enquanto.

Exemplos:

- Transformar coordenadas geográficas WGS84 em PTTM06.

import simplejson
import urllib

lat = 38

lng = -8

url = 'http://geodivagar.appspot.com/geogauss?lat=' + str(lat) + '&lng=' + str(lng)

response = simplejson.load(urllib.urlopen(url))

x = response['x']
y = response['y']

print x , y

- Na posse de coordenadas PTTM06 podemos calcular a altitude SRTM do ponto em metros.

url = 'http://geo-pt.appspot.com/srtmPT?x=' + str(x) + '&y=' + str(y) + '&interpol=bilinear'

response = simplejson.load(urllib.urlopen(url))

altitude = response['altitude']

print altitude

- Para as coordenadas geográficas podemos também calcular a ondulação do geóide em metros.

url = 'http://geodivagar.appspot.com/geoidePT?lat=' + str(lat) + '&lng=' + str(lng) + '&interpol=bilinear'

response = simplejson.load(urllib.urlopen(url))

n = response['N']

print n

Serviços disponíveis em http://geodivagar.appspot.com/

(Continua brevemente)

Web Service: Altitudes SRTM para Portugal Continental

Mais um serviço:

http://geo-pt.appspot.com/srtmPT? Devolve a altitude SRTM em metros para um determinado ponto em coordenadas x,y PT-TM06. Suporta 2 métodos de interpolação: vizinho mais próximo e bilinear.

- Exemplo: http://geo-pt.appspot.com/srtmPT?x=117006.46&y=-11722.69&interpol=bilinear
- Exemplo: http://geo-pt.appspot.com/srtmPT?x=117006.46&y=-11722.69&interpol=nn
- Parâmetros: x y coordenadas PT-TM06, interpol interpolador "nn" (nearest neighbor) ou "bilinear"
- Parâmetro opcional: callback para a chamada de uma função em javascript
- Resposta: JSON

Nota: na implementação deste serviço foi utilizado o modelo digital de terreno SRTM disponibilizado para download pelo Professor José Alberto Gonçalves. Fica a referência para quem quiser comparar os resultados do serviço.

Web Service: Códigos Postais de 7 dígitos

Saber a Localidade, Artéria, Local/Zona e Troço para um determinado código postal de 7 dígitos.

http://codigospostais.appspot.com/cp7? Devolve, quando existem, a Localidade, Artéria, Local/Zona e Troço para um determinado código postal de 7 dígitos.

• Exemplo: http://codigospostais.appspot.com/cp7?codigo=2620182
• Parâmetro: codigo - código postal de 7 dígitos
• Parâmetro opcional: callback para a chamada de uma função em javascript
• Resposta: JSON

Web Service: GeodPT08

O Instituto Geográfico Português (IGP) disponibiliza para download, o GeodPT08, um modelo do geóide para Portugal Continental elaborado pela FCUL em parceria com o IGP.

Este modelo é disponibilizado sob a forma de um ficheiro ASCII, o qual descreve uma grelha regular encontrando-se organizada no ficheiro segundo um esquema XYZ.

Do facto deste modelo requerer um programa com ferramentas de interpolação para ser utilizado, surgiu a ideia de publicar um web service que aceita pedidos GET e que devolve uma resposta em JSON. Assim foi feito.

/geoidePT? Devolve a ondulação do geóide (N) (GRS80) para uma determinada posição em WGS84/ETRS89 (lng e lat)

- Exemplo: http://geodivagar.appspot.com/geoidePT?lng=-10.0&lat=36.525&interpol=bilinear

- Exemplo: http://geodivagar.appspot.com/geoidePT?lng=-10.0125&lat=42.2375&interpol=nn

- Parâmetro: lat latitude em graus decimais

- Parâmetro: lng longitude em graus decimais

- Parâmetro: interpol. Valores possíveis: "bilinear" ou "nn" (nearest neighbor = vizinho mais próximo). Algoritmo a utilizar na interpolação na grelha.

- Resposta: JSON

Para mais considerações sobre o modelo consultar a página de download e também:
http://www.igeo.pt/produtos/Geodesia/GeodPT08/GeodPT08.pdf

Nota: assumiu-se que o espaçamento da grelha (0º,0250) é suficientemente pequeno para que possamos admitir que a ondulação do geóide tem um comportamento linear. Para a implementação do interpolador bilinear, ver por exemplo:
http://www.geocomputation.org/1999/082/gc_082.htm


Este e outros web services disponíveis em http://geodivagar.appspot.com/

NavPT foi à vida

Pois é. Segundo o próprio site:

"Por razões de ordem superior e completamente alheias à equipa do Projecto NavPT, somos obrigados a terminar o projecto."

Enquanto não soubermos mais pormenores acerca das "razões de ordem superior", resta-nos lembrar que se tratava de um projecto que permitia visualizar trajectórias, em tempo real, relativas à RIV (Região de Informação de Voo) de Lisboa. Permitia também, entre outras coisas, ter acesso às comunicações entre as aeronaves e a torre de controlo.

Ou seja, nada a que não possamos ter acesso através de sites internacionais, o que poderá levantar alguma suspeição sobre a legitimidade das alegadas "razões de ordem superior". Para mais considerações sobre este e outros aspectos, consultar por exemplo este blog.

Tile Map Service

Depois de termos os nossos tiles gerados, precisamos mesmo de um servidor de mapas? Acho que não (proof of concept).

Proj4js - transformação de coordenadas

Relativamente à aplicação para transformação pontual de coordenadas que disponibilizo aqui, desde o ano passado, foram feitas algumas (pequenas) alterações:

• Actualização dos parâmetros de Bursa-Wolf, sendo agora utilizados os parâmetros mais recentes do IGP;
• Formatação das coordenadas geográficas (graus decimais) de saída para 8 casas decimais;
• Formatação das coordenadas cartográficas (metros) de saída para 2 casas decimais.

Da página do IGP referida constam também pontos transformados os quais poderão servir para comparar resultados.

Link:
Transformação de coordenadas com Proj4js (implementação em javascript da biblioteca PROJ)

Serviços: parâmetro opcional

Todos os serviços disponibilizados em geodivagar.appspot.com admitem um parâmetro opcional callback para a chamada de uma função em javascript. Pode ser utilizado, por exemplo, o JSONScriptRequest.

Para isso, incluir o script na página.

Definir o callback na chamada ao serviço:

   function onMapClick(overlay,latlng) {

  if (latlng) {
   coords=latlng;
   request = 'http://geodivagar.appspot.com/geogauss?lat=' + latlng.y + '&lng=' + latlng.x + '&callback=onComplete';   
   aObj = new JSONscriptRequest(request);
   aObj.buildScriptTag();
   aObj.addScriptTag();   
       }
   }

Definir a função que recebe o objecto JSON, por exemplo:

 function onComplete(jData) {
   var myHtml = 'X: ' + jData.x.toFixed(2)+ '
 Y: ' + jData.y.toFixed(2);
           map.openInfoWindow(coords, myHtml);
 }

Demo completa disponível aqui.

Dica simples e eficaz

Directamente do blog do Darren Cope:

Merge A Directory of Shapefiles Using OGR

mkdir merged
for %f in (*.shp) do (
if not exist merged\merged.shp (
ogr2ogr -f “esri shapefile” merged\merged.shp %f) else (
ogr2ogr -f “esri shapefile” -update -append merged\merged.shp %f -nln Merged )
)

Mais 2 serviços JSON

Mais 2 serviços com resposta JSON:

/wmsdetails Informação detalhada acerca de um determinado serviço WMS

Exemplo:

http://geodivagar.appspot.com/wmsdetails?url=http://mapas.igeo.pt/wms/caop/continente

Parâmetros: url do serviço WMS

Resposta:

{
     layers: [
          "Distritos",
          "NUT3",
          "NUT2",
          "Freguesias",
          "CAOP",
          "Concelhos",
     ],
     operations: [
          "GetCapabilities",
          "GetMap",
          "GetFeatureInfo",
          "DescribeLayer",
          "GetLegendGraphic",
          "GetStyles",
     ],
     title: "Carta Administrativa Oficial de Portugal (CAOP - Continente) - Versão 6.0",
     abstract: "Serviço de mapas WMS: Continente - Limites Administrativos Oficiais (Limites de País, Limites de NUT 2, Limites de NUT3, Limites de Distrito, Limites de Concelho e Limites de Freguesia) para o Continente e Ilhas. A esta informação está associada a toponímia, bem como outra informação descritiva como seja a área oficial de cada circunscrição administrativa.",
     version: "1.1.1",
     type: "OGC:WMS"
}

/layerdetails Informação detalhada acerca de uma determinada layer de um serviço WMS

Exemplo:

http://geodivagar.appspot.com/layerdetails?url=http://mapas.igeo.pt/wms/caop/continente&layer=Distritos

Parâmetros: url do serviço WMS e layer

Resposta:

{
     boundingBoxWGS84: [
         -10.1913,
         36.8987,
         -5.7141,
         42.1887,
     ],
     styles: {
           default: {
                legend: http://mapas.igeo.pt/wms/caop/continente?version=1.1.1&service=WMS&request=GetLegendGraphic&layer=Distritos&format=image/png,
                title: "default"
           }
     },
     crsOptions: [
          "EPSG:27492",
          "EPSG:20790",
          "EPSG:4326",
          "EPSG:4258",
          "EPSG:102160",
          "EPSG:102161",
          "EPSG:102164",
          "EPSG:102165"
     ],
    title: "Distritos"
}

Foi utilizada a biblioteca OwsLib

Restantes serviços disponíveis em http://geodivagar.appspot.com/

Serviços no App Engine

Aproveitei o facto de ter criado conta no Google App Engine para disponibilizar alguns serviços com parâmetros passados na URL e com resposta em JSON (quem quiser chamar REST a isto, está à vontade).

Os serviços são os seguintes (poucos por enquanto):

/geogauss Projecção directa de coordenadas geográficas ETRS89 em coordenadas cartográficas PTTM06

Exemplo: http://geodivagar.appspot.com/geogauss?lat=38.01132704&lng=-7.87167688

Parâmetros: "lat"=latitude "lng"=longitude

Resposta: JSON

/gaussgeo Projecção inversa de coordenadas cartográficas PTTM06 em coordenadas geográficas ETRS89

Exemplo: http://geodivagar.appspot.com/gaussgeo?x=22958.6550821&y=-183907.735462

Parâmetros: "x" "y"

Resposta: JSON

/arcmer Comprimento de arco de meridiano (metros) sobre o elipsóide associado ao sistema WGS-84 entre duas latitudes

Exemplo: http://geodivagar.appspot.com/arcmer?lat1=38.0&lat2=39.0

Parâmetros: "lat1"=latitude inicial "lat2"=latitude final

Resposta: JSON

/elevation Gera uma URL do GOOGLE CHART com o perfil de terreno entre dois pontos

Example: http://geodivagar.appspot.com/elevation?lat1=38.5&lat2=39.0&lng1=-8.0&lng2=-8.0

Parâmetros: "lat1" e "lng1" coordenadas do ponto de partida, "lat2" e "lng2" coordenadas do ponto de chegada

Nota: a utilização deste serviço pressupôe a aceitação dos Termos de Serviço da API do Google Maps bem como dos Termos das Google Chart Tools

Resposta: JSON

/dms2decimal Converte graus, minutos e segundos em graus decimais

Exemplo: http://geodivagar.appspot.com/dms2decimal?d=39&m=40&s=5.73

Parâmetros: "d"=graus "m"=minutos "s"=segundos

Resposta: JSON

/decimal2dms Converte graus decimais em graus, minutos e segundos

Exemplo: http://geodivagar.appspot.com/decimal2dms?decimal=39.6682583333

Parâmetro: "decimal"= graus decimais

Resposta: JSON

/geoutil Faz parsing de coordenadas geográficas

Exemplo 1: http://geodivagar.appspot.com/geoutil?str=39%2040m%205.73s%20N%208%207m%2059s%20W

Exemplo 2: http://geodivagar.appspot.com/geoutil?str=39.5%20N%209.0%20W

Parâmetro: "str"= aceita strings como "39 40m 5.73s N 8 7m 59s W" ou "39.5 N 9.0 W"

Resposta: JSON

Ainda não implementei serviços para processamento em batch mas poderei vir a fazê-lo no futuro. Adicionalmente, poderei vir a disponibilizar mais serviços dependendo da disponibilidade de tempo.

Na sequência do 10º aniversário do fim do SA

Na sequência do 10º aniversário do fim do Acesso Selectivo (SA) o qual ocorreu dia 1 de Maio, encontrei uma página com as primeiras edições da coluna Innovation da revista GPS World:

• GPS: A multipurpose system by D. Wells and A. Kleusberg
  January/February 1990, Vol. 1, No. 1, pp. 60-63.


• The limitations of GPS by A. Kleusberg and R.B. Langley
  March/April 1990, Vol. 1, No. 2, pp. 50-52.


• Why is the GPS signal so complex? by R.B. Langley
  May/June 1990, Vol. 1, No. 3, pp. 56-59.


• Electronic charts and GPS by M.J. Casey and P. Kielland
  July/August 1990, Vol. 1, No. 4, pp. 56-59.


• The issue of selective availability by Y. Georgiadou and K.D. Doucet
  September/October 1990, Vol. 1, No. 5, pp. 53-56.


• Comparing GPS and GLONASS by A. Kleusberg
  November/December 1990, Vol. 1, No. 6, pp. 52-54.


• The GPS receiver: An introduction by R.B. Langley
  January 1991, Vol. 2, No. 1, pp. 50-53.


• Precise, real-time dredge positioning by S.R. DeLoach
  February 1991, Vol. 2, No. 2, pp. 43-45.


• The orbits of GPS satellites by R.B. Langley
  March 1991, Vol. 2, No. 3, pp. 50-53.


• Ionospheric effects on GPS by J.A. Klobuchar
  April 1991, Vol. 2, No. 4, pp. 48-51.


• GPS vehicle location and navigation by E.J. Krakiwsky
  May 1991, Vol. 2, No. 5, pp. 50-53.


• Continuous monitoring of crustal deformation by Y. Bock
  June 1991, Vol. 2, No. 6, pp. 40-47.


• The mathematics of GPS by R.B. Langley
  July/August 1991, Vol. 2, No. 7, pp. 45-50.


• GPS in civil aviation by K.D. McDonald
  September 1991, Vol. 2, No. 8, pp. 52-59.


• GPS – satellites of opportunity for ionospheric monitoring by D. Coco
  October 1991, Vol. 2, No. 9, pp. 47-50.

ECW Compressor

A ERMapper sempre disponibilizou no seu site utilitários gratuitos nomeadamente a famosa ferramenta de compressão ECW Compressor (com GUI e tudo, mas também com o limite de 500 MB por imagem). Depois da empresa ter sido adquirida pela ERDAS deixei de ver o compressor na secção de downloads no site desta última, apesar de lá constarem outros produtos herdados da ERMapper, nomeadamente o visualizador gratuito ER Viewer.

Assim, para quem não queira utilizar o binómio GDAL/ECW ou directamente o ECW SDK (o qual neste momento também não se encontra disponível) ou outra ferramenta que utilize o SDK, disponibilizo para download a versão ECW Compressor 2.6. Actualmente não tenho bem a certeza acerca do licenciamento desta ferramenta por isso a utilização da mesma fica ao critério de cada um (por outras palavras, não culpem o mensageiro).

LINK

Pesquisa por códigos postais de 7 dígitos

Normalmente, um código postal de 7 dígitos (cp7) corresponde a uma frente de quarteirão (pelo menos em zonas urbanas). Consequentemente, a uma rua podem estar associados vários cp7 tanto do mesmo lado da rua como do lado oposto. Por isso, a pesquisa por cp7, na maioria dos casos, não chega para definir de forma inequívoca uma morada sendo por isso necessário completar a informação com o respectivo número de policia. Em princípio o nome do arruamento estará subentendido no cp7.

No entanto, existe quem insista na pesquisa por códigos postais de 7 dígitos.

Por exemplo, a pesquisa por códigos postais dos Mapas do Sapo deixa muito a desejar. Quem pesquisar por meia dúzia de códigos postais, constata rapidamente que o que é feito é uma comparação entre o código postal e a descrição de cada ponto de interesse (POI) constante da base de dados. Uma abordagem ad hoc, portanto. Se não houver nenhum POI cuja morada contenha o código postal pretendido, o resultado será errático (por exemplo, pesquisando pelo código postal do Fórum Picoas, 1050-996, vamos obter um resultado próximo da Fundação Calouste Gulbenkian, o que não tem nada a ver) .

Abordagem mais interessante encontramos em geocoder.pinguimcomfrio.net. O autor construiu uma base de dados com os códigos postais e respectivos arruamentos. Quando se pesquisa por um determinado código postal, a respectiva morada (rua e localidade) constante da base de dados é enviada para o geocoder do Google. Com esta solução ao menos temos a garantia que o resultado se encontra na rua pretendida, o que já é muito bom.

JSON no browser (Firefox)

Para quem queira visualizar documentos em JSON o ideal é utilizar a extensão JSONView do Firefox. O documento será visualizado de forma semelhante a um documento XML, formatado, com os campos em destaque e podendo ser expandido. Sempre é melhor do que ter que fazer o download do documento JSON e visualizá-lo num editor de texto, numa única linha.

Para adicionar a última versão JSONView 0.4 clicar aqui.

Marcadores como resposta de um serviço por bounding box

Quando temos uma aplicação em Google Maps (ou OpenLayers, ou outra) e fazemos uma chamada a um serviço dando como parâmetros uma bounding box, geralmente utilizamos o resultado para adicionar marcadores ao mapa (com ou sem informação associada) em que as coordenadas dos marcadores a adicionar fazem parte da resposta do serviço.

Se fizermos chamadas sucessivas a esse serviço, por exemplo na sequência de um evento "moveend" registado, teremos que garantir que não vamos adicionar ao mapa marcadores já existentes.

Assim, de forma a que não existam marcadores duplicados é necessário percorrer a resposta do serviço (ou uma cópia desta) e testar se já existe algum marcador sobre o mapa com as mesmas coordenadas. Isso pode ser feito recorrendo a uma função a qual recebe como parâmetro as coordenadas a testar:

 function isOnTheMap(latlng) {
  for (var i=0; i < markers.length ; ++i) {
   if ( markers[i].getLatLng().equals(latlng)) {
    return true;
   }   
  }
  return false;
 }

Outra questão consiste em garantir que os marcadores que já estavam sobre o mapa e que não façam parte da bounding box, sejam removidos do mapa. Ou seja, garantir que não existem marcadores adicionados ao mapa fora da vieweport. Mais uma vez podemos recorrer a uma função:

 function removeMarkers() {
  var tmpMarkers = markers;
  for (var i = 0; i < tmpMarkers.length; ++i) {
   var marker = tmpMarkers[i];
   if (bounds.containsLatLng(marker.getLatLng())) {
    continue;
   }
   map.removeOverlay(marker);
   markers.splice(i,1); 
  } 
 }

Neste caso poderá também ser necessário garantir a remoção de eventuais eventos registados.

Para que ambos exemplos resultem é ainda necessário que aquando da criação de cada marcador seja este adicionado a um array que vá guardando os nossos marcadores:

markers.push(marker);

O sermão do topógrafo

Where 2.0 2010: Paul Ramsey, "Why your data sucks"

Clip raster with polygon (GDAL)

#!/bin/bash
SHPFILE=$1
BASE=`basename $SHPFILE .shp`
EXTENT=`ogrinfo -so $SHPFILE $BASE | grep Extent \
| sed 's/Extent: //g' | sed 's/(//g' | sed 's/)//g' \
| sed 's/ - /, /g'`
EXTENT=`echo $EXTENT | awk -F ',' '{print $1 " " $4 " " $3 " " $2}'`
echo "Clipping to $EXTENT"
RASTERFILE=$2
gdal_translate -projwin $EXTENT -of GTiff $RASTERFILE /tmp/`basename $RASTERFILE .sid`_bbclip.tif
gdalwarp -co COMPRESS=DEFLATE -co TILED=YES -of GTiff \
-r lanczos -cutline $SHPFILE \
/tmp/`basename $RASTERFILE .sid`_bbclip.tif /tmp/`basename $RASTERFILE .sid`_shpclip.tif

Fonte

Geoservicios

O ICC (Institut Cartogràfico de Catalunya), numa atitude bastante pragmática, utiliza LizardTech Express Server e...pronto. Sem entrar em considerações de ordem financeira, só posso dizer que estes Geoservicios são estremamente rápidos, veja-se por exemplo: Ortofoto de Cataluña 1:5.000 ou Base Topográfica 1:5.000. Também interessante é a quantidade de serviços wms disponíveis. Ver em Servicios disponibles. Na página da LizardTech encontra-se mais informação sobre este case study.

Cartographer.js

Cartographer.js é uma pequena lib em javascript a qual permite adicionar Pie-Charts, Choropleth Maps e Clusters ao Google Maps. É portanto uma biblioteca para cartografia temática.

Esta lib vem já com a divisão estadual dos EUA, hardcoded. Podemos sempre fazer o mesmo para os distritos de Portugal. Os concelhos poderão ser excessivamente pesados para representar do lado do cliente. E as freguesias ainda mais.

Pessoalmente, acho que faz mais sentido apresentar Pie-Charts e Choropleths (coropletos em Português) sobre mapas estáticos. Quem quer ver dados estatísticos geralmente não está interessado em ver o telhado das casinhas ou o nome das ruas. É uma questão de escala.

Relativamente aos Clusters prefiro utilizar o ClusterMarker.

De qualquer forma fica a referência a este curioso projecto.

O repositório encontra-se aqui.

Ilhas mais próximas do Continente

Após fazer o download da shapefile referente aos códigos postais de 4 dígitos, no site dos Correios, reparei que os Arquipélagos da Madeira e dos Açores estavam representados exageradamente próximos do Continente.

O problema deve-se à arbitrariedade de tentar representar no sistema Hayford-Gauss Militar (Continente), dados provenientes de data geodésicos regionais.

Como só estava interessado nos dados do Continente purguei a informação restante. No entanto, convinha que os Correios corrigissem este erro.

Termo de Utilização: OpenLayers e Google Maps

Têm-me perguntado frequentemente o porquê da utilização da API do Google Maps preterindo a API do OpenLayers.

A resposta é que sempre que pretendo utilizar dados de base do Google Maps (e pretendo muitas vezes), prefiro utilizar directamente a sua API que ir recorrer ao OpenLayers adicionando uma OpenLayers.Layer.Google .

Ou seja, quando utilizo OpenLayers e adiciono uma Layer de Google Maps, não deixo de estar sujeito aos Termos de Utilização da API do Google Maps. Isto parece-me trivial.

De qualquer forma, as restrições de utilização da API do Google Maps não são nada de transcendente, basta ver em http://code.google.com/intl/pt-PT/apis/maps/

"The Maps API is a free service, available for any web site that is free to consumers"

"To use the Maps API on an intranet or in a non-publicly accessible application, please check out Google Maps API Premier"

Altitudes do Brasil

Clicar no mapa para saber a altitude

Altitudes em metros. Coordenadas Geográficas referentes ao sistema WGS84.
Especificação dos dados:
- Resolução dos dados: 30 * 30 metros.
- Exactidão vertical: 20 metros (nível de confiança de 95%).
- Exactidão horizontal: 30 metros (nível de confiança de 95%).

Fonte dos dados: Geonames

Raster Query com a GDAL

Problema: saber o valor do pixel de uma imagem georeferenciada para um determinado par de coordenadas e uma determinada banda.

Solução utilizando GDAL Python bindings:

Nota: no código que se segue substituir underscore(_) por espaço em branco.

import struct
import gdal
from gdalconst import *

filename = "image.tif"
dataset = gdal.Open( filename, GA_ReadOnly )

def CellValue(dataset, band, x, y):
_band = dataset.GetRasterBand(band)
_if (band.GetNoDataValue() == None):
___band.SetNoDataValue(-9999)
_ndv = band.GetNoDataValue()
_cols = band.XSize
_rows = band.YSize
_geotransform = dataset.GetGeoTransform()
_cellSizeX = geotransform[1]
_cellSizeY = -1 * geotransform[5]
_minx = geotransform[0]
_maxy = geotransform[3]
_maxx = minx + (cols * cellSizeX)
_miny = maxy - (rows * cellSizeY)
_print 'bbox(real-world coords):',minx,',',miny,',',maxx,',',maxy
_if ((x <> maxx) or (y <> maxy)):
___print 'given point does not fall within grid'
___return ndv
_xLoc = (x - minx) / cellSizeX
_xLoc = int(xLoc)
_yLoc = (maxy - y) / cellSizeY
_yLoc = int(yLoc)
_#print 'point (pixels): ',xLoc,',',yLoc
_if ((xLoc <> cols - 0.5)):
___#print 'xcoordinate out of bounds'
___return ndv
_if ((yLoc <> rows - 0.5)):
___#print 'ycoordinate out of bounds'
___return ndv
_strRaster = band.ReadRaster(xLoc, yLoc, 1, 1, 1, 1, band.DataType)
_sDT = gdal.GetDataTypeName(band.DataType)
_if (sDT == 'Int16'):
___dblValue = struct.unpack('h', strRaster)
_elif (sDT == 'Float32'):
___dblValue = struct.unpack('f', strRaster)
_elif (sDT == 'Byte'):
___dblValue = struct.unpack('B', strRaster)
_else:
___print 'unrecognized DataType:', gdal.GetDataTypeName(band.DataType)
___return ndv
_print sDT
_return dblValue[0]

Testar a função:

print CellValue(dataset, 1, 222960.0, 285040.0)

Fonte: thread da mailing-list gdal-dev.

Layer de SQL Server 2008 no MapServer

O MapServer tem a possibilidade de estabelecer conexão (em modo read-only) com uma base de dados espacial residente em Microsoft SQL Server 2008.

No caso do MS4W começar por confirmar a existência do respectivo plugin:

/ms4w/Apache/specialplugins/msplugin_mssql2008.dll

No mapfile, na respectiva layer, é necessário utilizar o parâmetro CONNECTIONTYPE PLUGIN bem como o parâmetro PLUGIN para definir a PATH do mesmo.

Utilizar o parâmetro CONNECTION para definir os parâmetros de ligação que permitem aceder ao SQL Server e o parâmetro DATA para especificar a tabela na qual residem os dados espaciais.

LAYER
...
CONNECTIONTYPE PLUGIN
PLUGIN "C:/ms4w/Apache/specialplugins/msplugin_mssql2008.dll"
CONNECTION "server=localhost/sqlexpress;uid=dbusername;
             pwd=dbpassword;database=Districts Database;
               Integrated Security=false"
DATA "the_geom from districts"
TYPE POLYGON
STATUS ON
CLASS
   ...
END
END

Curvas de nível no MapServer

A representação de curvas de nível no MapServer não oferece problemas de maior. Apenas duas questões requerem um cuidado especial:

• diferenciação entre curvas de nível mestras e secundárias através de simbologia diferente;
• colocação de índices (labels) nas curvas de nível mestras.

Consideremos uma representação de curvas de nível com uma equidistância de 25 metros. Segundo as boas práticas de representação cartográfica, as curvas mestras serão representadas de 100 em 100 metros. Assim, teremos 3 curvas secundárias entre cada 2 curvas mestras.

No ficheiro de configuração do MapServer poderemos definir duas classes para a mesma layer obtendo algo do género:

MapServer em modo FastCGI

Uma das formas de optimizar o MapServer do lado do servidor web consiste em corrê-lo em modo FastCGI.

Instalação
(Em Windows com MS4W)

1. Actualmente o MS4W vem com suporte a FastCGI sendo o módulo do Apache (mod_fcgid) necessário carregado por defeito. Instalar o MS4W em:

"c:/ms4w/"

Verificar se está instalado o Visual C++ 2008 Redistributable Package

2. Editar o ficheiro de configuração do Apache httpd.conf

• remover o comentário da linha:
  

  LoadModule fcgid_module modules/mod_fcgid.so

• remover os comentários das seguintes linhas e colocar os caminhos correctos:
  

  
     DefaultInitEnv PROJ_LIB    "c:/ms4w/proj/nad/"
     DefaultInitEnv PATH        "c:/ms4w/Apache/cgi-bin;
       c:/WINDOWS/system32;c:/WINDOWS;c:/WINDOWS/System32/Wbem;"
     DefaultInitEnv windir      "c:/WINDOWS"
     DefaultInitEnv SystemRoot  "c:/WINDOWS"
     DefaultInitEnv SystemDrive "c:"
     DefaultInitEnv GDAL_DATA   "c:/ms4w/gdaldata"
     DefaultInitEnv GDAL_DRIVER_PATH "c:/ms4w/gdalplugins"
     DefaultInitEnv TMP         "c:/ms4w/tmp"
     DefaultInitEnv TEMP        "c:/ms4w/tmp"
  

• adicionar parâmetros de configuração do FastCGI como por exemplo:
  

  
     IPCCommTimeout 60
     IdleTimeout 60
     DefaultMinClassProcessCount 2
     DefaultMaxClassProcessCount 20
     DefaultInitEnv PROJ_LIB    "c:/ms4w/proj/nad/"
     DefaultInitEnv PATH        "c:/ms4w/Apache/cgi-bin;
       c:/WINDOWS/system32;c:/WINDOWS;c:/WINDOWS/System32/Wbem;"
     DefaultInitEnv windir      "c:/WINDOWS"
     DefaultInitEnv SystemRoot  "c:/WINDOWS"
     DefaultInitEnv SystemDrive "c:"
     DefaultInitEnv GDAL_DATA   "c:/ms4w/gdaldata"
     DefaultInitEnv GDAL_DRIVER_PATH "c:/ms4w/gdalplugins"
     DefaultInitEnv TMP         "c:/ms4w/tmp"
     DefaultInitEnv TEMP        "c:/ms4w/tmp"
  

• salvar o ficheiro httpd.conf e reiniciar o Apache.

3. Modificar a aplicação CGI de forma a apontar para "/fcgi-bin/mapserv.exe" em vez de "/cgi-bin/mapserv.exe".

4. Adicionar o parâmetro seguinte às layers para as quais queremos uma ligação FastCGI:

         PROCESSING "CLOSE_CONNECTION=DEFER"

5. Testar a aplicação. Se a instalação tiver sido bem sucedida então deveremos ver o processo "mapserv.exe" no Windows Task Manager, mantendo-se aberto enquanto o utilizador interage com a aplicação.

Notas

Em modo FastCGI poderá não se verificar um aumento de performance embora haja um ganho em utilização de recursos e eficiência.

Com layers de shapefiles o parâmetro "CLOSE_CONNECTION=DEFER" é simplesmente ignorado, não havendo portanto, um ganho de performance assinalável, por maiores que sejam as shapefiles.

O FastCGI compensa em situações em que exista um "custo" de ligação elevado, por exemplo, em ligações a ArcSDE, Oracle ou PostGIS.

O melhor será sempre fazer testes de benchmarking.

Serviços geográficos disponibilizados em Portugal

O Capítulo Português do OSGeo teve a feliz ideia de, através do grupo de trabalho OGC Web Services, criar e manter uma lista de webservices nacionais ou contendo dados nacionais.

Os serviços listados são basicamente os do IGP, SNIRH, Entre Douro e Vouga Digital, Esri Portugal (MDT de 30 metros já mencionado neste blog) e Atlas Interactivo do Alentejo, Algarve e Andaluzia. Destes apenas o IGP dispõe de serviços WFS. Não é muito mas é o que temos.

Podemos, no entanto, ensaiar uma pequena listagem de "candidatos a felizes contemplados":

• Web Map Service do Atlas de Portugal , cuja página diz "O serviço WMS do Atlas de Portugal encontra-se temporariamente indisponível". Substituir "temporariamente" por "indefinidamente". Não me lembro de ver este serviço disponível, num passado recente;
  

• Web Map Service do SNIT da DGOTDU o qual necessita de acreditação, pelo menos "nesta fase de desenvolvimento do sistema". Não vejo que desenvolvimento visto utilizarem uma solução out of the box;

• Geo-WebServices do IGeoE. Dividem-se em serviços pagos (carta militar 1:25.000 e 1:50.000) e gratuitos (carta militar 1:250.000 e 1:500.000). O IGeoE disponibiliza não só serviços WMS mas também ArcIMS. Exelente oportunidade para testar a classe OpenLayers.Layer.ArcIMS.
  

E é tudo, tirando alguns serviços publicados por acidente (por exemplo, WebSIGs municipais que consomem serviços wms "internos" que de interno têm muito pouco).

WMS no Maps

Há uns tempos atrás andei a testar as possibilidades de integração de Web Map Services (WMS) com a API do Google Maps.

Deixo aqui um pequeno exemplo. Neste exemplo foi utilizado o serviço WMS arcgis_online/MDT30m_PT_WGS84 da ESRI.

Para quem esteja interessado na temática WMS/Google Maps deixo as seguintes referências obrigatórias:

• Google Maps API and WMS Servers no blog do Brian Flood;
• How To Overlay Custom Maps Using the Google Maps API no blog do Kyle Mulka;
  
• John Decks's Blog (os links neste blog estão quebrados mas fica a referência na mesma);
• http://www.geoskating.com/gmap/.
  

Google maps: curiosidades

Servidores

Mapa (G_NORMAL_MAP)

• No mapa normal o Google utiliza quatro servidores para fazer balanceamento de carga: mt0, mt1, mt2 e mt3 (http://mt_.google.com);
• Os servidores devolvem tiles de dimensão 256*256 no formato PNG;

Satélite (G_SATELLITE_MAP)

• Para servir imagens de satélite o Google utiliza quatro servidores para fazer balanceamento de carga: kh0, kh1, kh2 e kh3 (http://kh_.google.com).
  
• O "kh" significa Keyhole;
• Os servidores devolvem tiles de dimensão 256*256 no formato JPG;

No modo híbrido G_HYBRID_MAP os dados são servidos tanto por servidores http://mt_.google.com como por servidores http://kh_.google.com. É utilizada uma combinação de satélite e mapa em que os dados do mapa são servidos em tiles transparentes.

Projecção dos tiles

Os tiles são projectados segundo a versão esférica da projecção de Mercator em que o raio da esfera é igual ao semi-eixo maior do elipsóide associado ao sistema WGS-84.

Na sintaxe da PROJ4 vamos ter:

+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext
+no_defs

A projecção dos tiles é importante se quisermos trabalhar com Tile Overlays.

Nos restantes casos estaremos sempre a trabalhar com coordenadas geográficas (referentes ao sistema WGS-84), seja para fazer o upload de um ficheiro KML nos My Maps, seja na utilização da API.

Por exemplo, quando definimos o centro de um mapa e o nível de zoom, temos:

map.setCenter(new GLatLng(38.5, -9.5), 7);

GLatLng é um ponto definido em coordenadas geográficas (latitude e longitude) no sistema WGS-84.

Shapefile

Em 1998 a ESRI publicou um white paper, com a descrição técnica do formato ESRI Shapefile.

Nesse momento passou a ser possível, sem constrangimentos de qualquer ordem, a leitura e escrita neste formato, o qual se tornou (já era?) no standard de facto no que diz respeito a formatos vectoriais de SIG.

Aliás, nessa descrição técnica, na secção How Shapefiles Can Be Created, é referido explicitamente, depois de mencionar a linha de produtos ESRI, "Write directly to the shapefile specifications by creating a program".

Na realidade, este facto foi dos que mais contribuíram para o aparecimento de software open-source GIS. Repare-se que, existem ferramentas open-source que utilizam exclusivamente a shapefile como formato vectorial nativo. Exemplos? Por exemplo, o MapWindow GIS.

O mesmo podemos dizer da introdução do formato ARC/INFO ASCII GRID, o qual devido à sua natureza (ficheiro de texto com Header), não necessitou de ver a sua especificação publicada para que fosse adoptado por terceiros.

Enfim, a César o que é de César.

Clique no mapa para saber a altitude

Faça zoom sobre a área pretendida, clique no mapa e será aberto um popup com as coordenadas e a respectiva altitude.

Altitudes obtidas a partir do Aster Global Digital Elevation Model através de um Webservice do tipo REST.
Resolução dos dados: 30 * 30 metros.
Exactidão vertical: 20 metros (nível de confiança de 95%).
Exactidão horizontal: 30 metros (nível de confiança de 95%).

Personal Geodatabase sem recorrer ao ArcGIS

É frequente a utilização do formato ESRI @ Personal Geodatabase (PGeo) como formato de transferência de dados geográficos. Também é frequente encontrar sites que disponibilizam dados nesse mesmo formato para download.

Se estivermos interessados apenas em visualizar esses dados, isso pode ser feito, por exemplo, recorrendo ao AccuGlobe® Desktop 2007. Curiosamente, o TatukGIS Viewer não o permite.

Se por outro lado, estivermos interessados em editar esses dados, podemos sempre utilizar a biblioteca OGR Simple Feature Library, mais especificamente o utilitário ogr2ogr. Nota: o driver PGeo tem que estar instalado. Podem recorrer à instalação FWTools para Windows.

Por exemplo, o comando seguinte irá criar uma pasta de nome aaa contendo uma shapefile por cada feature class contida na geodatabase aaa.mdb,

>> ogr2ogr -f "ESRI Shapefile" bbb aaa.mdb



OpenLayers: Ortofotomapas do IGP

Da listagem de serviços WMS disponibilizados pelo IGP não constam os "famosos" Ortofotomapas os quais podemos ver, por exemplo, no Bing Maps.

No entanto, o IGP dispõe de um visualizador de dados escrito em OpenLayers, do qual, a partir do respectivo código fonte, podemos extrair o endereço desses mesmos Ortofotomapas.

A saber, http://mapas.igeo.pt/tilecache/tilecache.py? , Layer: Ortos

Actualização (30 de Março de 2010) - Constatei hoje que este visualizador de OpenLayers não está disponível (sabe-se lá porquê). Alternativamente, disponibilizo o link seguinte o qual funciona como cliente dos Ortos. Nota: os Ortos só estão disponíveis para níveis de zoom superiores a 8, inclusive.

http://geomatica.comoj.com/ortos.html

Importar ESRI ® Shape files no AutoCAD

Import Shapes é um aplicativo que corre sobre Autocad, o qual permite converter Shape Files em DWG (conversão tanto de geometrias, SHP, como de atributos, DBF). O resultado é um ficheiro Autocad normal, o qual é também 100% compatível com Autocad Map 2008.

Front-end gráfico para ogr2ogr

ogr2ogr é um utilitário para conversão e manipulação de dados vectoriais. Faz parte do projecto OGR (sub-projecto da GDAL, "Geospatial Data Abstraction Library”).

ogr2gui é um aplicativo com interface gráfico o qual permite converter dados sem recorrer à sintaxe complexa do ogr2ogr. É um projecto open-source.