OSGeo Planet

GeoSolutions: Developer’s Corner: dynamic raster palettes and other spatiotemporal extensions for GeoServer

OSGeo Planet - Mon, 2016-07-25 14:23


Dear Reader,

today we would like to introduce a new community module that we have been working on: ncWMS like extensions to perform easy dynamic mapping over raster data.

The first change introduced by the extension is the ability to easily add dynamic color palettes that will be applied on single banded raster data based on their statistics, but with the ability to control range and palette application from request parameters too. For example, let's say we want to create a simple red to blue progression, with the ncWMS extension installed we can simply create a new style, choose the "dynamic palette" format, and enumerate the two colors:

redblue-editor Then, let's make sure the layer band configuration has a range consistent with the data, in this case, a very narrow one:


Finally, we can make a simple request that will use the dynamic palette against the configured range:


You can also try dynamically on this link: http://cloudsdi.geo-solutions.it/geoserver/wms?STYLES=,redblue&LAYERS=countries,flexpart&FORMAT=image%2Fpng8&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=-137,15,-110,42&WIDTH=600&HEIGHT=600

The dynamics of the data appears to favor the low values, it's possible to enhance the high portion by using a logarithmic palette instead, by adding "&logscale=true" to the request, this way: http://cloudsdi.geo-solutions.it/geoserver/wms?STYLES=,redblue&LAYERS=countries,flexpart&FORMAT=image%2Fpng8&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=-137,15,-110,42&WIDTH=600&HEIGHT=600&logscale=true


It's also possible to concentrate on a specific range using the "&colorscalerange=0.00001,0.0002", as follows: http://cloudsdi.geo-solutions.it/geoserver/wms?STYLES=,redblue&LAYERS=countries,flexpart&FORMAT=image%2Fpng8&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=-137,20,-110,42&WIDTH=600&HEIGHT=500&colorscalerange=0.00001,0.0002


Finally, one can generated an animation providing a time or elevation range, asking for the GIF format, and adding the "&animation=true" (here, with a different, more complex palette): http://cloudsdi.geo-solutions.it/geoserver/wms?STYLES=,x-Sst&LAYERS=countries,flexpart&FORMAT=image%2Fgif&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&SRS=EPSG%3A4326&BBOX=-137,20,-110,42&WIDTH=600&HEIGHT=500&time=2016-02-01/2016-03-01&animation=true&format_options=gif_loop_continuosly:true;gif_frames_delay:200


The official documentation provides more insight into these new features, while the module can be already downloaded as part of the GeoServer 2.10.x nightly builds (it's a development series):

If you want to know more about how we can help your organization, don't hesitate and get in touch with us! Make sure to check our Enterprise Services Offer which has been used to complete the work described in this post.

The GeoSolutions Team,

Categories: OSGeo Planet

Fernando Quadro: Análise e visualização de dados LiDAR – Parte 5

OSGeo Planet - Mon, 2016-07-25 10:30

Para o nosso exemplo de hoje vamos dar um zoom e encontrar um edifício para analisar (se você aumenta o zoom, os edifícios se tornam clicáveis).


Podemos ver que o edifício #81719 (identificador da chave primária do banco de dados) tem um número limitado de atributos, incluindo uma elevação de 1.446,43! Nossos dados LIDAR variam cerca de 450 metros, de modo que a elevação nos edifícios é provavelmente medida em pés.

Como podemos calcular a elevação para o edifício #81719 utilizando a tabela LIDAR? Veja como a lógica funciona visualmente:

1. Comece com o edifício.


2. Localize todos os patches que se cruzam naquele edifício.


3. Transforme essas manchas em pontos individuais.


4. Filtre esses pontos usando o contorno do edifício.


5. Finalmente a média dos pontos de elevações dará um valor para o edifício.

Veja o procedimento acima em SQL:

-- We run a set of subqueries sequentially using
-- the "with" keyword
-- Get the one building we are interested in
building AS (
  SELECT geom FROM buildings
  WHERE buildings.gid = 81719
-- All the patches that intersect that building
patches AS (
  SELECT pa FROM medford
  JOIN building ON PC_Intersects(pa, geom)
-- All the points in that patch
pa_pts AS (
  SELECT PC_Explode(pa) AS pts FROM patches
-- All the points in our one building
building_pts AS (
  SELECT pts FROM pa_pts JOIN building
  ON ST_Intersects(geom, pts::geometry)
-- Summarize those points by elevation
  Avg(PC_Get(pts, 'z')) AS lidar_meters
FROM building_pts;

O resultado do SQL é 441.075 metros, que em pés é 1447.1, quase exatamente o mesmo valor do arquivo de edifícios, 1.446,43!

Ok! Porém podemos acrescentar a elevação a partir da nossa imagem LIDAR para todos os nossos edifícios? Sim, mas vai demorar algum tempo de processamento. Primeiro vamos adicionar uma coluna para aceitar o valor, e então executar uma atualização na nossa tabela do banco de dados.

-- Add column for our calculate Z value
ALTER TABLE buildings ADD COLUMN z real;
-- Update into the column
UPDATE buildings SET z = elevs.z
  -- For every building, all intersecting patches
  WITH patches AS (
      buildings.gid AS buildings_gid,
      medford.id AS medford_id,
      medford.pa AS pa
    FROM medford
    JOIN buildings
    ON PC_Intersects(pa, geom)
  -- Explode those patches into points, remembering
  -- which building they were associated with
  pa_pts AS (
    SELECT buildings_gid, PC_Explode(pa) AS pts FROM patches
  -- Use the building associations to efficiently
  -- spatially test the points against the building footprints
  -- Summarize per building
    Avg(PC_Get(pts, 'z')) AS z
  FROM pa_pts
  JOIN buildings
  ON buildings.gid = buildings_gid
  WHERE ST_Intersects(buildings.geom, pts::geometry)
  GROUP BY buildings_gid
) AS elevs
-- Join calculated elevations to original buildings table
WHERE elevs.buildings_gid = gid;

Nossa tabela está atualizada! Confira as elevações originais e calculadas (temos que converter a coluna z de metros para pés para compará-la com a coluna elevação):

SELECT z AS z_m, z*3.28084 AS z_ft, elevation AS elevation_ft
FROM buildings;

Os valores estão bem próximos, veja:

   z_m | z_ft | elevation_ft 
--------- + ------------------ + --------------- 
 438,128 | 1.437,42663560303 | 1434.43000000 
 433,556 | 1.422,4291678418 | 1413.63200000 
 435,489 | 1.428,76987573853 | 1406.25400000 
 439,244 | 1.441,09014682129 | 1422.58200000 
 433,738 | 1.423,02460105347 | 1416.93600000 
 429,648 | 1.409,60687857788 | 1403.92400000 
 437,264 | 1.434,59264585083 | 1425.84000000 
 430,607 | 1.412,75115040894 | 1404.03675300

Vamos colocar a diferença de elevação (em metros) na tabela, veja como:

-- Add column for our calculated height
ALTER TABLE buildings ADD COLUMN height real;
-- Update the building heights by subtracting elevation of
-- the nearest street from the elevation of the building
UPDATE buildings SET height = heights.height
  WITH candidates AS (
      b.gid AS building_gid,
      s.gid AS street_gid,
      s.z AS streets_z,
      b.z as buildings_z
    FROM buildings b, streets s
    WHERE ST_DWithin(b.geom, s.geom, 0.001)
      ST_Distance(b.geom, s.geom)
  SELECT DISTINCT ON (building_gid)
    building_gid, street_gid,
    buildings_z - streets_z AS height
  FROM candidates
) AS heights
WHERE heights.building_gid = buildings.gid;

No próximo post, finalizaremos nossa série de posts sobre LiDAR, não perca!

Posts RelacionadosZemanta
Categories: OSGeo Planet

Free and Open Source GIS Ramblings: OSM turn restriction QA with QGIS

OSGeo Planet - Sat, 2016-07-23 16:10

Wrong navigation instructions can be annoying and sometimes even dangerous, but they happen. No dataset is free of errors. That’s why it’s important to assess the quality of datasets. One specific use case I previously presented at FOSS4G 2013 is the quality assessment of turn restrictions in OSM, which influence vehicle routing results.

The main idea is to compare OSM to another data source. For this example, I used turn restriction data from the City of Toronto. Of the more than 70,000 features in this dataset, I extracted a sample of about 500 turn restrictions around Ryerson University, which I had the pleasure of visiting in 2014.

As you can see from the following screenshot, OSM and the city’s dataset agree on 420 of 504 restrictions (83%), while 36 cases (7%) are in clear disagreement. The remaining cases require further visual inspection.


The following two examples show one case where the turn restriction is modelled in both datasets (on the left) and one case where OSM does not agree with the city data (on the right).
In the first case, the turn restriction (short green arrow) tells us that cars are not allowed to turn right at this location. An OSM-based router (here I used OpenRouteService.org) therefore finds a route (blue dashed arrow) which avoids the forbidden turn. In the second case, the router does not avoid the forbidden turn. We have to conclude that one of the two datasets is wrong.

turn restriction in both datasets missing restriction in OSM?

If you want to learn more about the methodology, please check Graser, A., Straub, M., & Dragaschnig, M. (2014). Towards an open source analysis toolbox for street network comparison: indicators, tools and results of a comparison of OSM and the official Austrian reference graph. Transactions in GIS, 18(4), 510-526. doi:10.1111/tgis.12061.

Interestingly, the disagreement in the second example has been fixed by a recent edit (only 14 hours ago). We can see this in the OSM way history, which reveals that the line direction has been switched, but this change hasn’t made it into the routing databases yet:

now before

This leads to the funny situation that the oneway is correctly displayed on the map but seemingly ignored by the routers:


To evaluate the results of the automatic analysis, I wrote a QGIS script, which allows me to step through the results and visually compare turn restrictions and routing results. It provides a function called next() which updates a project variable called myvar. This project variable controls which features (i.e. turn restriction and associated route) are rendered. Finally, the script zooms to the route feature:

def next(): f = features.next() id = f['TURN_ID'] print "Going to %s" % (id) QgsExpressionContextUtils.setProjectVariable('myvar',id) iface.mapCanvas().zoomToFeatureExtent(f.geometry().boundingBox()) if iface.mapCanvas().scale() < 500: iface.mapCanvas().zoomScale(500) layer = iface.activeLayer() features = layer.getFeatures() next()

You can see it in action here:

I’d love to see this as an interactive web map where users can have a look at all results, compare with other routing services – or ideally the real world – and finally fix OSM where necessary.

This work has been in the making for a while. I’d like to thank the team of OpenRouteService.org who’s routing service I used (and who recently added support for North America) as well as my colleagues at Ryerson University in Toronto, who pointed me towards Toronto’s open data.

Categories: OSGeo Planet

MapProxy: New MapProxy 1.9.0 release

OSGeo Planet - Fri, 2016-07-22 19:00

We are pleased to announce the release of MapProxy 1.9.0. It contains a lot of major and minor improvements.

The latest release is available at: https://pypi.python.org/pypi/MapProxy

To upgrade within your virtualenv:

$ pip install --upgrade --no-deps MapProxy

Updated documentation is available at: https://mapproxy.org/docs/1.9.0/

Support for ArcGIS REST services

You can now directly integrate ArcGIS services without the need to enable WMS support.

See: https://mapproxy.org/docs/1.9.0/sources.html#arcgis-label

Band merging

The new band merge feature allows you to create false-color or grayscale images on the fly, by selecting different sources for each (color) band.

See https://mapproxy.org/docs/1.9.0/configuration.html#band-merging and https://mapproxy.org/docs/1.9.0/configuration_examples.html#create-grayscale-images

Other fixes and changes

There are many more changes and improvements. For a complete list of see: http://github.com/mapproxy/mapproxy/blob/1.9.0/CHANGES.txt

Categories: OSGeo Planet

QGIS Blog: QGIS 2.16 ‘Nødebo’ is released!

OSGeo Planet - Fri, 2016-07-22 14:41

We’re happy to announce the release of QGIS 2.16.0 ‘Nødebo’. The University of Copenhagen’s Department of Geoscience and Natural Resource Management Forest and Landscape College in Nødebo were hosts to the First International QGIS conference and developer meeting in May 2015. For all of us who are not fluent in Danish, Lene Fischer has prepared the following video teaching us how to pronounce the release name:

QGIS 2.16 is not designated as a Long Term Release (LTR). Users wishing to have a version of QGIS which does not change and receives bug fixes for at least 1 year are invited to use the current LTR release 2.14.
If you are upgrading from QGIS 2.14 you will find a great many new features in this release.
Whenever new features are added to software they introduce the possibility of new bugs – if you encounter any problems with this release, please file a ticket on the QGIS Bug Tracker.

We would like to thank the developers, documenters, testers and all the many folks out there who volunteer their time and effort (or fund people to do so). From the QGIS community we hope you enjoy this release! If you wish to donate time, money or otherwise get involved in making QGIS more awesome, please wander along to qgis.org and lend a hand!

QGIS is supported by donors and sponsors. A current list of donors who have made financial contributions large and small to the project can be seen on our donors list. If you would like to become and official project sponsor, please visit our sponsorship page for details. Sponsoring QGIS helps us to fund our six monthly developer meetings, maintain project infrastructure and fund bug fixing efforts.

QGIS is Free software and you are under no obligation to pay anything to use it – in fact we want to encourage people far and wide to use it regardless of what your financial or social status is – we believe empowering people with spatial decision making tools will result in a better society for all of humanity. If you are able to support QGIS, you can donate here.

Categories: OSGeo Planet

Fernando Quadro: Fontes de dados geográficos gratuitos

OSGeo Planet - Fri, 2016-07-22 10:30

Quem trabalha com Geotecnologia sabe a dificuldade de encontrar dados geográficos. Apesar da INDE estar em “pleno funcionamento”, ainda estamos engatinhando nesse sentido. Temos poucos portais que disponibilizam informações geográfica se comparado ao tamanho do nosso país.

O LabGIS, após uma vasta pesquisa, compilou uma extensa base de 587 links com dados geográficos gratuitos para consulta e download. Essa base está agora aberta ao público, podendo ser acessada por qualquer pessoa.

Caso você conheça alguma fonte de dados que não consta nesta listagem você pode informá-los por meio desse link, ajudando assim a aumentar essa base pública e disponível a toda a comunidade.

Fonte: Sistema LabGIS

Posts RelacionadosZemanta
Categories: OSGeo Planet

Jackie Ng: React-ing to the need for a modern MapGuide viewer (Part 2): The baseline

OSGeo Planet - Thu, 2016-07-21 14:17
I've decided that cataloging my adventure in developing this viewer is enough to warrant its own blog mini-series. So here's part 2 of my adventure. I have no idea how many parts this will take :)

So to continue where we left off, we had a basic map viewer component built on the main pillars of
And it is all written in glorious TypeScript. I've mentioned once on a previous post that writing web frontend code in a language with static typing and a compilation phase (gasp!), with a clean and simple-to-understand component model offered by React is just heavenly! If I could go back in time, I'd have written the Fusion framework in TypeScript, but sadly it didn't exist around 2008 (when Fusion made its debut in MapGuide) and my involvement in the MapGuide project was nowhere near the level it is now.
So before I continue, I should also probably mention what I aim to achieve with this viewer:
  • It will be built on modern, best-of-breed web front-end technologies (React, OpenLayers 3, etc)
  • I will not be shackled with the burden of supporting legacy browsers. This viewer will demand a modern web browser. If one has to use Internet Explorer, it must be IE11
  • It will require MapGuide Open Source 3.0 as the minimal version as this viewer will leverage capabilities present in this particular version onwards.
  • It will at least have some level of feature parity with the AJAX/Fusion viewer offerings, that is to say:
    • It will have a Legend component for toggling layer visibility and comprehending layer symbology/thematics.
    • It will have a Task Pane component that functions as a generic container for content and/or contextual UI
    • It will have a Map Selection component for easy viewing of selected features
    • It will include a stock set of commands, with extension points for developers to add their own.
    • It will have opt-in integration with external base layers (OpenStreetMap, etc), provided your Map Definition meets the requirements (ie. It is in EPSG:3857/WGS84.PseudoMercator).
    • I will not at this time consider any kind of integration with Google Maps as their APIs are not only a constantly moving target, but a TOS minefield I just do not want to wander in at this point in time.
    • I do hope to allow for some level of free-form component layout (ala. templates) so one is not tied to a layout I choose by default.
So with that being said, for this post I want to focus on the last major bullet point: Getting to parity with our current viewer offerings.

So let's introduce our new Legend component.

Implementing this component wasn't difficult, I mainly copypasta'd most of the presentation logic from the existing feature-complete legend used in my various mapguide-rest sample applications.

The difference here is that because this is now a React component, there is zero jQuery or DOM manipulation. Instead, we are just render Group/Layer child components for each layer and group we encounter in our runtime map. It is all pure UI as a function of props and state, and letting React's virtual DOM figure out how to actually update the DOM in the most efficient manner. What this means is that our component model cleanly matches how our layers and groups are conceptually structured.

As for the other components, I'm still working on them (currently focused on the Task Pane) so we'll leave that for another post. But visually speaking, when you put these 3 items together:

It's already starting to look like our existing AJAX viewer!

And before I close out this post, In terms of our current JS bundle size, here's the results of the current weigh-in (using our current Fusion production bundle for comparison).

We've got plenty of breathing room.
Categories: OSGeo Planet

Fernando Quadro: Problema na internacionalização do GeoServer

OSGeo Planet - Thu, 2016-07-21 10:30

A partir da versão 2.8.3 foi adicionado ao GeoServer o idioma português (pt-BR), porém, ao que tudo indica ele foi criado a partir da versão em espanhol e disponibilizado sem que a tradução estivesse completa, ou seja, existem termos em português, espanhol e alguns até em inglês.

Se você está passando por esse problema, saiba que é possível retirar esses arquivos da sua versão do GeoServer e fazer com que ela volte a exibir as informações em inglês.

Para isso, basta seguir o passo-a-passo abaixo:

1. Renomeie o arquivo geoserver.war para geoserver.zip
2. Entre na pasta WEB-INF\lib
3. Encontre o arquivo gs-web-core-X.X.X.jar e renomei-o para .zip
4. Procure o arquivo de properties (GeoServerApplication_pt_BR.properties).
5. Exclua o arquivo
6. Renomeie o arquivo gs-web-core-X.X.X.zip para gs-web-core-X.X.X.jar
6.1. Repita os passos 3, 4, 5 e 6 para os arquivos:

  • gs-web-demo-X.X.X.jar
  • gs-web-sec-core-X.X.X.jar
  • gs-web-gwc-X.X.X
  • gs-web-wcs-X.X.X.jar
  • gs-web-wfs-X.X.X.jar
  • gs-web-wms-X.X.X.jar

7. Renomeie o arquivo geoserver.zip para geoserver.war.
8. Agora está pronto, basta subir sua aplicação e verificar se está tudo em inglês novamente.

Caso esteja utilizando a versão binária, executável (Windows) ou DMG (Mac), desconsidere os passos 1 e 7, e lembre-se que onde está indicado X.X.X você deve substituir pelo número da versão do GeoServer que você está utilizando no momento.

Qualquer problema, não deixe de comentar nesse post.

Posts RelacionadosZemanta
Categories: OSGeo Planet

Free and Open Source GIS Ramblings: One “add” button to rule them all

OSGeo Planet - Wed, 2016-07-20 19:46

Reducing the number of “Add layer” buttons in the QGIS GUI is a commonly voiced wish. Multiple approaches have been discussed but no decision has been made so far. One idea is to use the existing browser functionality to replace the “Add layer” dialogs. Others are envisioning completely novel approaches.

Since the topic came up again today on Twitter, I decided to implement a quick & dirty version of a unified Add layer button. This way, I can comfortably reduce my Layer toolbar to three buttons using Settings | Customization …



I pretty much just kept the “Create new layer” button and the “Add delimited text layer” button because, as far as I know, there is no way to call the dialog from the browser. (Instead, CSVs are opened with OGR, which doesn’t have nearly as many nice features.)

And here it is in action:

(I recommend to undock the Browser panel to get the dialog-like behavior that you see in the video.)

To install the plugin: download it and unzip it into your QGIS plugin folder, then activate it in the plugin manager.

I would love to hear what you think about this UX experiment.

Categories: OSGeo Planet

gvSIG Team: “Introduction to gvSIG” free online workshops

OSGeo Planet - Wed, 2016-07-20 11:20

The gvSIG Association is now pleased to inform the “Introduction to gvSIG” summer webinars at the beginning of August.

Summer has arrived in a lot of places around the world. This is a good time to learn to use gvSIG, the open source GIS.

Do you want to create a thematic map with gvSIG or a new vector layer with your parcel, a road…? You can do it now at these open and free webinars.

There will be two webinars, that will last about 30 minutes:

At the first webinar we’ll see how to manage a gvSIG project. We’ll create new views with cartography, and we’ll apply symbology and labelling. We’ll introduce the reference systems, and finally we will create a thematic map to be printed or exported to a PDF file.

At the second webinar, we’ll show how to create a new vector layer, where we’ll digitalize several elements. We also will see the gvSIG toolbox, where we can find all the geoprocesses. We will apply one of them. We finally will see how the 3D extension works.

Registration for the “Getting started with gvSIG” webinar

Registration for the “Vector editing, 3D view and other gvSIG tools” webinar

We expect your participation!

Filed under: community, english, events, gvSIG Desktop, training Tagged: introduction to gvSIG, webinar
Categories: OSGeo Planet

gvSIG Team: 3as Jornadas gvSIG México: geomática y software libre. Propuestas para comunicaciones

OSGeo Planet - Tue, 2016-07-19 12:36

Se ha ampliado el plazo para la recepción de propuestas de comunicaciones para las 3as Jornadas gvSIG México, que se celebrarán en Ciudad de México del 7 al 9 septiembre de este año. La nueva fecha límite es el próximo 5 de agosto.

Para ello no tenéis más que enviar un resumen siguiendo la plantilla facilitada en el apartado de Comunicaciones de la web a la dirección 3asjornadasgvsig@igg.unam.mx. El resumen será valorado por el Comité Científico de cara a su inclusión en el programa de las Jornadas. Existen dos modalidades de comunicación: ponencia y póster.

banner-3as_J_Mex_gvSIGPor otro lado, las inscripciones, que son gratuitas, siguen abiertas, y pueden realizarse desde el formulario disponible en la propia web.

Durante las jornadas se realizarán un buen número de talleres para usuarios y desarrolladores. Tened en cuenta que las plazas son limitadas, no esperéis a inscribiros a última hora.

¡Esperamos vuestra participación!

Filed under: community, events, spanish, training Tagged: jornadas, México
Categories: OSGeo Planet

Even Rouault: Speeding up computation of raster statistics using SSE-2/AVX-2

OSGeo Planet - Tue, 2016-07-19 11:25
GDAL offers a method ComputeStatistics() that given a raster band returns the minimum and maximum values of pixels, the mean value and the standard deviation.

For those not remembering how to compute mean and standard deviations, the basic formulas for values indexed from 0 to N-1 are :mean = sum(value(i) for i = 0 to N-1) / Nstd_dev = square root of the mean of the square of the differences of values to the mean
std_dev = sqrt(sum(i = 0 to N-1, (value(i) - mean)^2)) / N)A very naive version would first compute the mean, and in a second pass compute the standard deviation.
But it can be easily proven (by expanding the (value(i) - mean)^2 term),that it is also equivalent to :
std_dev = sqrt(sum(i = 0 to N-1, value(i)^2)/N - mean^2)
std_dev = sqrt(mean_of_square_values - square_of_mean)

std_dev = sqrt(sum(i = 0 to N-1, value(i)^2)/N - (sum_of_values/N)^2)
std_dev = sqrt(N^2 *(sum(i = 0 to N-1, value(i)^2)/N - (sum_of_values/N)^2)) / N
std_dev = sqrt(N * sum_of_square_values - sum_of_values^2) / NA less naive implementation would compute the sum of values and the sum of square values in a single pass. However the standard deviation computed like that might be subject to numeric instability given that even if the result is small, sum_of_square_values and sum_of_values can be very big for a big number of pixels, and thus if represented with floating point numbers, the difference between both terms can be wrong.
Welford algorithmSo in recent GDAL versions, the computation of the mean and standard deviation is done in a progressive and numerically stable way, thanks to the Welford algorithm
The generic code is:
pixel_counter = 0
mean = 0
M2 = 0
foreach(value in all pixels):
    if value < minimum or pixel_counter == 0: minimum = value
    if value > maximum or pixel_counter == 0: maximum = value
    pixel_counter = pixel_counter + 1
    delta = value - mean
    mean = mean + delta / pixel_counter
    M2 = M2 + delta * (value - mean);

std_dev = sqrt( M2 / pixel_counter )
Proof of Welford algorithm(You can skip this paragraph and still follow the rest of this article)

The magic of Welford algorithm lies in the following recurrence relations.

For the mean, it is rather obvious :

N*mean(N) = sum(i = 0 to N-1, value(i))
N*mean(N) = sum(i = 0 to N-2, value(i)) + value(N-1)
N*mean(N) = (N-1) * mean(N-1) + value(N-1)
mean(N) = (N-1)/N * mean(N-1) + value(N-1)/N
mean(N) = mean(N-1) + (value(N-1) - mean(N-1)) / N
Hence mean = mean + delta / pixel_counter

For the standard deviation, the proof is a little bit more lengthy :

N*stddev(N)^2 = sum(i=0 to N-1, (value(i) - mean(N))^2 )

N*stddev(N)^2 = sum(i=0 to N-1, (value(i) - (mean(N-1) + (value(N-1) - mean(N-1)) / N))^2 )

N*stddev(N)^2 = sum(i=0 to N-1, ((value(i) - mean(N-1)) - ((value(N-1) - mean(N-1)) / N))^2 )

N*stddev(N)^2 = sum(i=0 to N-1, (value(i) - mean(N-1))^2 + ((value(N-1) - mean(N-1)) / N)^2
             - 2 * (value(i) - mean(N-1))*((value(N-1) - mean(N-1)) / N)  )

N*stddev(N)^2 = sum(i=0 to N-1, (value(i) - mean(N-1))^2) + N * ((value(N-1) - mean(N-1)) / N)^2
              - 2 * sum(i=0 to N-1, (value(i) - mean(N-1)))*((value(N-1) - mean(N-1)) / N))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 +  (value(N-1) - mean(N-1)) ^2
                    +  N * ((value(N-1) - mean(N-1)) / N)^2
              - 2 * sum(i=0 to N-1, (value(i) - mean(N-1)))*((value(N-1) - mean(N-1)) / N))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 + (value(N-1) - mean(N-1))^2 * (1 + 1 / N)
              - 2 * N( mean(N) - mean(N-1)) *((value(N-1) - mean(N-1)) / N))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 + (value(N-1) - mean(N-1)) *
            ((1 + 1 / N) *  (value(N-1) - mean(N-1)) - 2 * N( mean(N) - mean(N-1)) / N))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 + (value(N-1) - mean(N-1)) *
            ((value(N-1) - mean(N-1) + (value(N-1) - mean(N-1) / N - 2 * N( mean(N) - mean(N-1)) / N))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 + (value(N-1) - mean(N-1)) *
            ((value(N-1) - mean(N-1) - (mean(N) - mean(N-1))))

N*stddev(N)^2 = (N-1)*stddev(N-1)^2 + (value(N-1) - mean(N-1)) * (value(N-1) - mean(N))

Hence M2 = M2 + delta * (value - mean)

Integer based computation of standard deviationThe Welford algorithm is good but it involves floating point operations for each pixel to compute the progressive mean and variance. Whereas fundamentally we would need those floating point operations only at the end if using the original formulas, and we could use integer arithmetics for the rest. Another drawback of Welford approach is that it prevents any direct parallelization (there might still be ways to reconcile partial computations, but I have not explored those), whereas if you have a set of pixels, you can conceptually divide it in as many subsets you want, and for each subset compute its local minimum, maximum, sum of values and sum of square values. Merging subsets is then trivial: take the minimum of minimums, maximum of maximums, sums of sum of values and sums of sum of square values.
Let us consider the case of pixels whose type is unsigned byte, ie with values in the range [0,255]. We want to computestd_dev = sqrt(N * sum_of_square_values - sum_of_values^2) / NFor practical reasons, we want N, sum_of_square_values and sum_of_values to fit on a 64bit unsigned integer (uint64), which is the largest natural integral type that can be easily and efficiently used on today's CPUs. The most limiting factor will be sum_of_square_values. Given that in the worse case, a square value is equal to 255*255, the maximum number of pixels N we can address is (2^64-1) / (255*255) = 283 686 952 306 183, which is large enough to represent a raster of 16 million pixels x 16 million pixels. Good enough.
We know need to be able to multiply two uint64 values and get the result as a uint128, and compute the difference of two uint128 values. The multiplication on Intel/AMD CPUs in 64bit mode natively yields to a 128 bit wide result. It is just that there is no standardized way in C/C++ how to get that result. For GCC compiler in 64 bit mode, the __uint128_t type can be used in a transparent wayto do that :
__uint128_t result = (__uint128_t)operand_64bit * other_operand_64bitFor Visual Studio compilers in 64 bit mode, a special instruction _umul128() is available.
What about non-Intel or non-64bit CPUs ? In that case, we have to do the multiplication at hand by decomposing each uint64 values into its lower uint32 and uint32 parts, doing 4 uint32*uint32->uint64 multiplications, summing the intermediary results, handling the carries and building the resulting number. Not very efficient, but we do not really care about that, since it is just a final operation.
To make it is easier, that partial 128 bit arithmetics is abstracted in a GDALUInt128 C++ class that has different implementations regarding the CPU and compiler support.
Now that we have solved the final part of the computation, we can then writethe computation loop as following :

    minimum = maximum = value[0]
    foreach value:
        if value < minimum: minimum = value
        else if value > maximum: maximum = value
        sum = sum + value
        sum_square = sum_square + value * value

Can we do better ? A bit of loop unrolling can help :

    minimum = maximum = value[0]
    foreach value pair (value1, value2):
        if value1 < minimum: minimum = value1
        else if value1 > maximum: maximum = value1
        sum = sum + value1
        sum_square = sum_square + value1 * value1
        if value < minimum: minimum = value2
        else if value > maximum: maximum = value2
        sum = sum + value2
        sum_square = sum_square + value2 * value2
    (deal with potential remaining pixel if odd number of pixels)

If we start with comparing value1 and value2, we can actually save a comparison (resulting in 3 comparisons for each pair of pixel, instead of 4) :
    minimum = maximum = value[0]
    foreach value pair (value1, value2):
        if value1 < value2:
            if value1 < minimum: minimum = value1
            if value2 > maximum: maximum = value2
            if value2 < minimum: minimum = value2
            if value1 > maximum: maximum = value1
        sum = sum + value1
        sum_square = sum_square + value1 * value1
        sum = sum + value2
        sum_square = sum_square + value2 * value2
    (deal with potential remaining pixel if odd number of pixels)

This improvement can already dramatically reduce the computation time from
1m10 to 7s, to compute 50 times the statistics on a 10000 x 10000 pixel raster.
Parallelization with SSE2We have not yet explored the parallelization of the algorithm. One way to do it would be to use multi-threading, but for Intel-compatible CPU, we can also explore the capabilities of the SIMD (Single Instruction/Multiple Data) instruction set. On 64bit Intel, the SSE2 instruction set, which offers vectorized operations on integers, is guaranteed to be always present. 16 registers (XMM0 to XMM15) are available, each 128 bit wide.
So each register is wide enough to hold 16 packed int8/uint8, 8 packed int16/uint16, 4 packed int32/uint32 or 2 packed int64/uint64, depending on the wished representation. A diverse set of operations are offered and generally operate on the sub-parts of each register independently. For example c=_mm_add_epi8(a,b) will add independently c[i]=a[i]+b[i] for i=0 to 15, and that in just one CPU cycle._mm_add_epi16() will work on packed uint16, etc. To add some salt, not all operators are available for all elementary subtypes however.
Compilers are supposed to be able to automatically vectorize some C code, but in practice they rarely manage to do so for real world code, hence requiring the programmer to use the SIMD instruction set at hand. All major compilers (gcc, clang, Visual Studio C/C++) offer access to the SSE2 instruction set through "intrinsics", which are C inline functions that wrap the corresponding assembly instructions, but while still being C/C++. This allows the compiler to do the register allocation and various other optimizations (such as re-ordering), which is a huge win over coding directly in assembly. The Intel intrinsics guide is a useful resource to find the appropriate intrinsics.
So a temptative vectorized version of our algorithm would be :
    v_minimum = vector_of_16_bytes[0]
    v_maximum = vector_of_16_bytes[0]
    v_sum = vector_of_16_zeros
    v_sum_square = vector_of_16_zeros

    foreach vector_of_16_bytes v:
        v_minimum = vector_minimum(v_minimum, v)
        v_maximum = vector_maximum(v_maximum, v)
        v_sum = vector_add(v_sum, v)
        v_sum_square = vector_sum(v_sum_square, vector_mul(v, v))

    minimum = minimum_of_16_values(v_minimum)
    maximum = maximum_of_16_values(v_minimum)
    sum = sum_of_X??_values(v_sum)
    sum_square = sum_of_X??_values(v_sum_square)
    (deal with potential remaining pixels if number of pixels is not multiple of 16)

vector_minimum and vector_maximum do exist as _mm_min_epu8 and _mm_max_epu8. But for vector_add, which variant to use _mm_add_epi8, _mm_add_epi16, _mm_add_epi32 or _mm_add_epi64 ? Well, none directly. We want to add uint8 values, but the result cannot fit on a uint8 (255+255=510). The same holds for sum_square. The result of each square multiplication requires at least a uint16, and we want to loop several times, so we need at least a width of uint32 to hold the accumulation. We designed the overall algorithm to be able to handle an accumulator of uint64, but this would decrease the performance of the vectorization if using that in the tigher loop. So we will decompose our loop into one upper loop and and one inner loop. The inner loop will do as many iterations as possible, while still not overflowing a uint32 accumulator. So (2^32-1)/(255*255) = 66051.xxxx iterations. Which we round down to the closest multiple of 16.
So what about v_sum = vector_add(v_sum, v) ?
The first idea would be to extract the 4 lowest order bytes of v, unpack them so that they fit each on a uint32 and then use _mm_add_epi32 to add them in the v_sum accumulator.
    v_sum = _mm_add_epi32(v_sum, _mm_unpacklo_epi16(_mm_unpacklo_epi8(v, zero), zero)
_mm_unpacklo_epi8(v, zero) expands the 8 lowest order bytes of v as 8 uint16. And similarly _mm_unpacklo_epi16(v, zero)  expands the 4 lowest order uint16 of v as 4 uint32.
And then repeat that with the 3 other groups of 4 bytes :

    v_sum = _mm_add_epi32(v_sum, _mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_shuffle_epi32(v, 1), zero), zero)
    v_sum = _mm_add_epi32(v_sum, _mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_shuffle_epi32(v, 2), zero), zero)
    v_sum = _mm_add_epi32(v_sum, _mm_unpacklo_epi16(_mm_unpacklo_epi8(_mm_shuffle_epi32(v, 3), zero), zero)

But we can do better thans to the _mm_sad_epu8 intrinsics. It is designed to "compute the absolute differences of packed unsigned 8-bit integers in a and b, then horizontally sum each consecutive 8 differences to produce two unsigned 16-bit integers, and pack these unsigned 16-bit integers in the low 16 bits of 64-bit elements in dst." If we notice that ABS(x-0) = x when x >= 0, then it does what we want.
    v_sum = _mm_add_epi64(v_sum, _mm_sad_epu8(v, zero))

Pedantic note: we can actually use _mm_add_epi32, since there is no risk of overflow : 8 * 66051 * 255 fits on a uint32. The advantage of using _mm_add_epi32 is that as we will use it elsewhere, the compiler can re-order additions to group them in pairs and benefit from their 0.5 throughput.
_mm_sad_epu8() has a relatively high latency (5 cycles), but it is still a big win since it replaces 14 intrinsics of our initial version.
What about the computation of the square value ? There is no mnemonics to directly multiply packed bytes and get the resulting packed uint16 (or even better uint32, since that is the type we want to operate on eventually to be able to do several iterations of our loop!). One approach would be to take the 8 lowest order bytes, un-pack them to uint16, use the  _mm_mullo_epi16() intrinsics that does uint16 x uint16->uint16. Then you would take the 4 lowest order uint16 of this intermediate result, un-pack them to uint32 and finally use _mm_add_epi32 to accumulate them in v_sum_square.
    v_low = _mm_unpacklo_epi8(v, zero)
    v_low_square = _mm_mullo_epi16(v_low, v_low)
    v_sum_square = _mm_add_epi32(v_sum_square, _mm_unpacklo_epi16(v_low_square, zero)

Then repeat the operation with the 4 upper order uint16 of the intermediate result.
    v_sum_square = _mm_add_epi32(v_sum_square,
        _mm_unpacklo_epi16(_mm_shuffle_epi32(v_low_square, 2 | (3 <<2)), zero) )

_mm_shuffle_epi32(v, 2 | (3 <<2) is a trick to replicate the high 64 bits of a XMM register into its low 64 bits. We don't care about the values of the resulting high 64 bits since they will be lost with the later unpack operations.
And then repeat the whole process with the 8 highest order bytes.

    v_high = _mm_unpacklo_epi8(_mm_shuffle_epi32(v, 2 | (3 <<2)), zero)
    v_high_square = _mm_mullo_epi16(v_high, v_high)
    v_sum_square = _mm_add_epi32(v_sum_square, _mm_unpacklo_epi16(v_high_square, zero)
    v_sum_square = _mm_add_epi32(v_sum_square,
        _mm_unpacklo_epi16(_mm_shuffle_epi32(v_high_square, 2 | (3 <<2)), zero) )

We can actually do much better with the _mm_madd_epi16() mnemonics that "Multiply packed signed 16-bit integers in a and b, producing intermediate signed 32-bit integers. Horizontally add adjacent pairs of intermediate 32-bit integers, and pack the results". This is really close to what we need. We just need to prepare uint16/int16 integers (the sign convention here does not matter since a uint8 zero-extended to 16 bit is both a uint16/int16)
    v_low_16bit = _mm_unpacklo_epi8(v, zero)
    v_sum_square = _mm_add_epi32(v_sum_square, _mm_madd_epi16(v_low_16bit, v_low_16bit))
    v_high_16bit = _mm_unpacklo_epi8(_mm_shuffle_epi32(v, 2 | (3 <<2)), zero)
    v_sum_square = _mm_add_epi32(v_sum_square, _mm_madd_epi16(v_high_16bit, v_high_16bit))

The latencies and throughput of _mm_mullo_epi16 and _mm_madd_epi16 are the same, so the second version is clearly a big win.
Use of AVX2We can tweak performance a bit by doing a 2x loop unrolling, which will enable the compiler to re-order some operations so that those who have a throughput of 0.5 cycle to be consecutive (such as _mm_add_epi32, _mm_unpacklo_epi8) and thus be able to executive 2 of them in a single cycle. When doing so, we can notice that we are operating on a virtual 256 bit register. But 256 bit registers do exist in the AVX2 instruction set, that was introduced in relatively recent hardware (2013 for Intel Haswell). AVX/AVX2 offer the YMM registers, equivalent of XMM registers but on a doubled bit width (the 128 bit low part of a YMM register is its corresponding XMM register). One particularity of the YMM register is that it operates on quite distinct "128 bit lanes", but you can still extract each lane.
The port to AVX2 is quite straightforward :

    v = _mm256_load_si256(data + i)
    v_sum = _mm256_add_epi32(v_sum, _mm256_sad_epu8(v, zero))
    v_low_16bit = _mm256_cvtepu8_epi16(_mm256_extracti128_si256(v, 0));
    v_sum_square = _mm256_add_epi32(v_sum_square, _mm256_madd_epi16(v_low_16bit, v_low_16bit))
    v_high_16bit = _mm256_cvtepu8_epi16(_mm256_extracti128_si256(v, 1));
    v_sum_square = _mm_add_epi32(v_sum_square, _mm_madd_epi16(v_high_16bit, v_high_16bit))

_mm256_extracti128_si256(v,0) extracts the 128 bit lower part of the register,
and _mm256_extracti128_si256(v,1) the 128 bit upper part.
The good news is that we can have a single code base for the SSE2 and AVX2 variants, by using the AVX2 code. In the case of SSE2, we in fact define the _mm256 functions with their corresponding _mm 128bit functions operating on the low and high 128 bit parts. For example:
static inline GDALm256i GDALmm256_add_epi32(GDALm256i r1, GDALm256i r2)
    GDALm256i reg;
    reg.low = _mm_add_epi32(r1.low, r2.low);
    reg.high = _mm_add_epi32(r1.high, r2.high);
    return reg;

The AVX2-with-SSE2 emulation can be found in :

Thanks to inlining and usual compiler optimizations, this will be equivalent to our hand 2x unrolled version ! The final code is here.

Regarding timings, our SSE2-emulated AVX2 version runs in 1.6s, so roughly a 4x time improvement with respect to the portable optimized C version. On a hardware capable of AVX2, the pure AVX2 version is 15% faster than the SSE2-emulated version. So definitely not enough to justify a dedicated code path, but here as we have a unified one, it comes almost for free. Provided that the code is explicitly compiled to enable AVX2.
Nodata valuesUp to now, we have ignored the potential existence of nodata values. When computing statistics, we do not want pixels that match the nodata value to be taken into account in the minimum, maximum, mean or standard deviation.
In the pure C approach, this is easy. Just ignore pixels that match the nodata value:
    minimum = maximum = value[0]
    foreach value:
        if value != nodata:
            valid_pixels = valid_pixels + 1
            minimum = min(minimum, value)
            maximum = max(minimum, value)
            sum = sum + value
            sum_square = sum_square + value * value

We cannot directly translate that with SSE2/AVX2 mnemonics since the result of the value != nodata test can be different for each of the 32 packed bytes of the (virtual) AVX2 register, and making tests for each components of the vector register would kill performance to a point where it would be worse than the pure C approach !
We can however rewrite the above in a vector friendly way with :
    minimum = maximum = first value that is not nodata
    neutral_value = minimum (or any value in final [min,max] range that is not nodata)
    foreach value:
        validity_flag = if (value != nodata) 0xFF else 0
        value_potentially_set_to_zero = value & validity_flag
        value_potentially_set_to_neutral = (value & validity_flag) | (neutral_value & ~validity_flag)
        valid_pixels = valid_pixels + validity_flag / 255
        minimum = min(minimum, value_potentially_set_to_neutral)
        maximum = max(minimum, value_potentially_set_to_neutral)
        sum = sum + value_potentially_set_to_zero
        sum_square = sum_square + value_potentially_set_to_zero * value_potentially_set_to_zero

(value & validity_flag) | (neutral_value & ~validity_flag) is a quite common pattern in SIMD code to implement a if/then/else pattern without branches (for classic scalar code, if/then/else branches are more efficient due to the CPU being able to do branch prediction)
The only complication is that there is no SSE2 intrinsics for non-equality testing, so we have to transform that a bit to use equality testing only. And we will also remove the need for division in the loop :
    foreach value:
        invalidity_flag = if (value == nodata) 0xFF else 0
        value_potentially_set_to_zero = value & ~invalidity_flag
        value_potentially_set_to_neutral = (value & ~invalidity_flag) | (neutral_value & invalidity_flag)
        invalid_pixels_mul_255 = invalid_pixels_mul_255 + invalidity_flag
        minimum = min(minimum, value_potentially_set_to_neutral)
        maximum = max(minimum, value_potentially_set_to_neutral)
        sum = sum + value_potentially_set_to_zero
        sum_square = sum_square + value_potentially_set_to_zero * value_potentially_set_to_zero

    valid_pixels = total_pixels - invalid_pixels_mul_255 / 255

The computation of invalid_pixels_mul_255 in a vectorized way is the same as
v_sum, using the _mm_sad_epu8() trick. The resulting SSE2 code is :
    foreach vector_of_16_bytes v:
        v_invalidity_flag = _mm_cmpeq_epi8(v, v_nodata)
        v_value_potentially_set_to_zero = _mm_andnot_si128(v_invalidity_flag, v)
        v_value_potentially_set_to_neutral = _mm_or_si128(
            v_value_potentially_set_to_zero, _mm_and_si128(v_invalidity_flag, v_neutral))
        v_invalid_pixels_mul_255 = _mm_add_epi32(invalid_pixels_mul_255,
                                        _mm_sad_epu8(v_invalidity_flag, zero))
        [ code for min, max operating on v_value_potentially_set_to_neutral ]
        [ code for sum and sum_square operating on v_value_potentially_set_to_zero ]

The transposition to AVX2 is straightforward.
We can notice that this version that takes into account nodata value can only be used once we have hit a pixel that is not the nodata value, to be able to initialize the neutral value.
What about uint16 rasters ?
The same general principles apply. If we still want to limit ourselves to operate with at most uint64 accumulators, given that the maximum square value of a uint16 is 65535*65535, this limits to rasters of 2^64/(65535*65535) ~= 2 billion pixels, which remains acceptable for common use cases.
One oddity of the SSE-2 instruction set is that it includes only a _mm_min_epi16() / _mm_max_epi16() mnemonics, that is to say that operates on signed int16. The _mm_min_epu16() that operates on uint16 has been introduced in the later SSE 4.1 instruction set (that is quite commonly found in not so recent CPUs).
There are tricks to emulate _mm_min_epu16() in pure SSE2 using saturated subtraction and masking :
    // if x <= y, then mask bits will be set to 1.
    mask = _mm_cmpeq_epi16( _mm_subs_epu16(x, y), zero )

    // select bits from x when mask is 1, y otherwise
    min(x,y) = _mm_or_si128(_mm_and_si128(mask, x), _mm_andnot_si128(mask, y));

Another way is to shift the unsigned values by -32768, so as to operate on signed 16bit values.
This -32768 shift trick is also necessary since, like for the byte case, we want to still be able to use the _madd_epi16 intrinsics, which operates on signed int16, to compute the sum of square values. One subtelty to observe is that when you operate on 2 consecutive pixels at 0, _mm_madd_epi16 does :
 (0 - 32768) * (0 - 32768) + (0 - 32768) * (0 - 32768)
= 1073741824 + 1073741824
= 2147483648 = 2^31

Which actually overflows the range of signed int32 ( [-2^31, 2^31-1] ) ! The good news is that _mm_madd_epi16 does not saturate the result, so it will actually return 0x80000000 as a result. This should normally be interpreted as -2^31 in signed int32 convention, but as we know that the result of _madd_epi16(x,x) is necessary positive values, we can still correctly interpret the result as a uint32 value. This is where you feel lucky that Intel chose two's complement convention for signed integers.
To compute the sum of values, no nice trick equivalent to _mm_sad_epu8. So we just do it the boring way: unpack separately the 64bit low and high part of the value register from uint16 to uint32 and accumulate them with _mm_add_epi32.
Exactly as for the byte case, the uint16 case can be easily transposed to AVX2 oremulated-AVX2.

Conversion between integer and floating-point operations can be costly, so avoiding them as much as possible is a big win (provided that you make sure not to overflow your integer accumulators)
In theory, the gains offered by a SSE2/AVX2 optimized version are at most limited to a factor of with_of_vector_register / with_of_elementary_type, so, for bytes and SSE2, to 16. But often the gain is lesser, so do that only when you have come to an already optimized portable C version (or if the SIMD instruction set includes a dedicated intrinsics that just do what you want)
Categories: OSGeo Planet

Fernando Quadro: Análise e visualização de dados LiDAR – Parte 4

OSGeo Planet - Tue, 2016-07-19 10:30

Agora vamos adicionar ao mapa uma camada de edifícios. Para isso descompacte o arquivo BuildingFootprints.zip que baixado anteriormente.

Nós vamos carregar os edifícios no banco de dados como uma tabela espacial, mas antes disso, nós temos de descobrir o qual o “SRID” usar.

Para descobrir o SRID basta abrir o arquivo BuildingFootprints.prj, veja:

PROJCS [ "NAD_1983_StatePlane_Oregon_South_FIPS_3602_Feet_Intl" , 
  GEOGCS [ "GCS_North_American_1983" , 
    DATUM [ "D_North_American_1983" , 
      SPHEROID [ "GRS_1980" , 6378137.0 , 298.257222101 ]], 
    PRIMEM [ "Greenwich" , 0.0 ], 
    UNIT [ "Degree" , 0.0174532925199433 ]], 
  PROJECTION [ "Lambert_Conformal_Conic" ], 
  PARAMETER [ "False_Easting" , 4921259.842519685 ], 
  PARAMETER [ "False_Northing" , 0.0 ], 
  PARAMETER [ "Central_Meridian" , - 120.5 ], 
  PARAMETER [ "Standard_Parallel_1" , 42.33333333333334 ], 
  PARAMETER [ "Standard_Parallel_2" , 44.0 ], 
  PARAMETER [ "Latitude_Of_Origin" , 41.66666666666666 ], 
  UNIT [ "Foot" , 0.3048 ]]

É muito claro que a informação que estamos procurando é esta “NAD_1983_StatePlane_Oregon_South_FIPS_3602_Feet_Intl”, mas qual número “SRID” devemos usar?

Para descobrir basta ir para ao site Prj2EPSG (http://prj2epsg.org), e procurar pela informação acima, e você vai obter a resposta que é 2270.

Agora, usando o shp2pgsql podemos carregar os dados em uma tabela chamada buildings:

shp2pgsql -s 2270 -D BuildingFootprints.shp buildings | psql -d lidar

Os nossos dados LIDAR estão todos em coordenadas geográficas (EPSG: 4326) e nós vamos ter de integrá-los com a camada buildings (EPSG:2270) para envitar problemas, sendo assim vamos utilizar a função ST_Transform do PostGIS pra isso:

-- Update SRID and transform all geoms
ALTER TABLE buildings
TYPE geometry(MultiPolygon,4326)
USING ST_Transform(geom, 4326);
-- Rename to area column to area
ALTER TABLE buildings
RENAME COLUMN shape_st_1
TO shape_area;
-- Index the table
CREATE INDEX buildings_gix ON buildings USING GIST (geom);

Para agilizar a nossa análise, vamos eliminar todos os edifícios que não estão contidos na área do nosso dado LiDAR.

-- Find the LIDAR extent
SELECT st_extent(pa::geometry) FROM medford;
-- BOX(-122.8874999 42.3125,-122.8749998 42.325)
-- Delete unneeded building polygons
DELETE FROM buildings
WHERE NOT ST_Contains(
  ST_MakeEnvelope(-122.8874999, 42.3125, -122.8749998, 42.325, 4326),

Agora, é publicar esta camada no GeoServer. Faça a publicação da forma tradicional, apenas na aba Cache desmarque a opção “Criar camada de cache para esta camada” e na aba Publishing altere as informações das propriedade do KML confirme a figura abaixo:


Salve, e então temos agora uma camada visível! Você pode visualizá-la no Google Earth usando o KML reflector do GeoServer.


Quando você ampliar a imagem no Google Earth, vai notar algo estranho sobre nossos edifícios: eles são “lisos”! Queremos edifícios em 3D, como podemos obtê-los? É isso que iremos ver nos próximos posts.

Posts RelacionadosZemanta
Categories: OSGeo Planet

GeoServer Team: Online GeoServer Bug Stomp

OSGeo Planet - Mon, 2016-07-18 18:56


Dear Readers,

a quick post to spread the word about the Online GeoServer Bug Stomp which will take place this Friday, the 22nd of July 2016.

Developers as well as users from the GeoServer community will gather online to spend up to a full day (in their timezone) on tasks like:

  • Reviewing JIRA Reports to make sure they are valid
  • Fix bugs as we come across them
  • Improved docs
  • Test and report new bugs or close existing reports

The rules of engagement as well as a first rough list of participants can be found in this document; the event is an online gathering, people will be working from their place coordinating using the GeoServer Gitter channel with whomever will be online in their timezone.

If you feel like helping don’t be scared, jump onboard, read the rules of engagement say hi on the gitter channel (github login required) and help us make GeoServer even better.

If you cannot, don’t worry, we are going to try and make this event a monthly event.

Happy GeoServer to everybody.


Categories: OSGeo Planet

Fernando Quadro: Análise e visualização de dados LiDAR – Parte 3

OSGeo Planet - Mon, 2016-07-18 10:30

O problema com o LIDAR é o seu tamanho! Este conjunto de dados de exemplo é muito, muito pequeno e inclui 11 milhões de pontos. Desenhá-los em um mapa seria uma tarefa muito lenta, e visualmente não é muito instrutivo.

Para exibir um resumo dos nossos dados, vamos colocar os patches de ponto do banco de dados em um mapa.

O GeoServer só pode mapear geometrias do PostGIS, por isso vamos usar a geometria do pcpatch para criar uma visão dos dados (polígono).

CREATE VIEW medford_patches AS
  pa::geometry(Polygon, 4326) AS geom,
  PC_PatchAvg(pa, 'Z') AS elevation
FROM medford;

Agora vamos configurar uma camada no GeoServer a partir da visão acima (medford_patches).


Quando estiver configurando a camada, na aba Tile/Cache desmarque a opção “Criar camada de cache para esta camada”.

Após salvar a camada, vá em Layer Preview e veja como ficou:


Você deve estar vendo um grande quadrado preto, como demonstra a figura acima, mas se ampliar o zoom em poucos passos as coisas se tornam mais claras.


Existem 28547 manchas e, quando tudo é desenhado em uma pequena prévia, eles parecem com uma massa escura. Mas ampliando, podemos ver o detalhe dos pequenos blocos de pontos, cada um com cerca de 400 pontos dentro.

No entanto, como são manchas escuras, deixam muito a desejar! Seria melhor se tivessem cor de acordo com sua elevação, por isso precisamos adicionar um estilo no GeoServer com as cores que desejamos.

Vamos então usar o arquivo elevation_ramp.xml e criar um novo estilo no GeoServer com o nome elevation_ramp e vinculá-lo a nossa camada.

Este não é um modelo padrão de SLD, ele não tem regras que definem as quebras de cor. Em vez disso, ele usa um recurso de estilo interpolado para criar uma efeito de cores contínuo, utilizando as quebras de cores sugeridas neste post.

Após a aplicação do estilo na camada, ela ficará da seguinte forma:


Agora a (pequena) variação da elevação pode ser vista. Na verdade, a variação é tão pequena temos dificuldade até de identificar grandes edifícios comerciais.

No próximo post iremos inserir as edificações (Buildings) no mapa. Não Perca.

Posts RelacionadosZemanta
Categories: OSGeo Planet

Free and Open Source GIS Ramblings: Traveltime routing & catchment extension for the QGIS IDF router

OSGeo Planet - Sun, 2016-07-17 20:26

As announced in Salzburg a few days ago, I’m happy to present the lastest enhancement to my IDF router for QGIS: travel time routing and catchment computation.

Travel times for pedestrians and cyclists are computed using constant average speeds, while car travel times depend on the speed values provided by the road network data.

Catchment computations return the links that can be traversed completely within the given time (or distance limit). The current implementation does not deal with links at the edge of the catchment area, which can only be traversed partially.

Loading the whole network (2.7GB unzipped IDF) currently requires around 10GB of memory. One of the next plans therefore is to add a way to only load features within a specified bounding box.

Plans to turn this into a full-blown plugin will most likely have to wait for QGIS 3, which will ship with Python 3 and other updated libraries.

Screenshot 2016-07-17 22.04.54

Categories: OSGeo Planet

GIScussions: Local Government Shared Services – wrestling with data, PostGIS and stuff

OSGeo Planet - Fri, 2016-07-15 18:20

This making maps business is certainly not as easy as it seems! I have just burnt an inordinate amount of time trying to make a half decent map and publish it to the web, but I am learning stuff on the way even if I did disappear up some blind alleys.


2 years ago, I wrote a piece “Could you make a better map than this?” which questioned the way in which the Local Government Association were presenting shared services partnerships on a map. I challenged people to make a better map with the offer of a small prize but no one took up my offer. A few weeks ago, the LGA published an updated version of the map and also provided a link to download the data, well the map is pretty similar to the original one that I criticised, have a look:

For me this map is cluttered with too many markers and colours and numbers on the markers that are confusing. So I thought this time I will try to make something, hopefully better. I am not sure that I have succeeded as representing this data was more difficult than I imagined and I think I have stretched the limits of my technical skills and of the QGIS2Web plugin. So here is a summary of what I tried, maybe it will help someone else.

Wrestling with the data

The data is in a csv file with almost 1700 rows, 1 row per organisation per service/partnership (so if 20 organisations are in a partnership there are 20 rows with pretty much identical data). Some simple metadata about the columns and their significance would have been helpful (e.g. there was a count field that didn’t seem to make any sense to me).

Lesson 1 when you are cleaning up data don’t delete stuff that you think you don’t need too early, you might discover why it was useful later on and have to restart!

I spent a while looking at the data and decided that it might be more effective to represent each service as a multi-polygon covering the Local Authorities. So the first task was to give each partnership a unique ID, how do you do that in Excel? A bit of googling and I cam up with this simple answer: Insert a column (A in this case) for the Service ID then sort the data by service name (column B), set A1=1 then use this formula =if(B2=B1, A1, A1+1) in A2 and copy down the column. Now you have a series of Service IDs, copy and paste as values so that they are fixed and don’t recalculate as you resort or delete some of the data.

SQL is hard

Teaching yourself SQL involves a heck of a lot of trial and error, if I had’t had a lot of help from my friend Aileen Heal (a database guru IMHO) I would never have got a result. Here is the sequence of stuff that I did:

  1. I imported the data into PostGIS and filtered out the services that had shut down or where there was no information to show that it was live to create a table of live services (ca 1100 rows)
  2. The data had an ONS Code for almost every row, so I grabbed the latest version of BoundaryLine from the OS OpenData site and loaded the Counties and Districts/Boroughs/Unitaries into PostGIS
  3. Then I appended the Counties data to the Districts so that I had one table with all of the boundaries and ONS codes
  4. Then I added a geometry column to the table of live services and updated the column with the geometry from the combined boundaries table using the ONS Code (setting the coordinate system and adding a spatial index)
  5. Then I created a table that created a multi polygon for each row within a group of a service ID. I wouldn’t have had a clue how to do this and couldn’t find an example on line so it was good that I could “phone a friend”.
  6. Then I could join the attribute data from the live table to the union table that had the multi polygons to create a view with 256 live services displayed as multi polygons.
  7. The final step was to work out how to split these 256 rows into a manageable number of categories to portray on the map. The original data had 27 different categories which was way too much (see the massive rainbow legend at the bottom of the LGA map above). I worked out how to run a Select, Count, Group and Order query to list all of the categories and the number of rows in each, pretty pleased with that (google is fantastic for working this stuff out). I ended up creating 9 higher level groupings that made some sense to me (e.g combining “Shared Management” with “Shared Leadership”) using select queries.

Job done! Well not first time, there were several iterations before I got it all right

Categories: OSGeo Planet

Jackie Ng: React-ing to the need for a modern MapGuide viewer

OSGeo Planet - Fri, 2016-07-15 17:41
A few months ago, I mused about a hypothetical future map viewer for MapGuide.
Having a good solid 6 months of working with React and TypeScript (and dealing with Webpack in anger) under my belt, I felt confident enough to give this idea a crack.
This is what the consumption story looks like so far:

The HTML is just:
  • One script tag (to the webpack bundle, containing the map viewer component and all of its dependencies
  • Making an instance of the entry point class (currently called MapGuide.Application)
  • Mounting the class (using React terminology here) at the desired DOM element with properties to pass into the root React component
And voila! Instant OpenLayers3-powered map viewer. Pretty simple stuff.
When we look at the current React component composition:

It's one root component (I'm intending for this to be something conceptually similar to the ApplicationDefinition/WebLayout for our existing viewers), and a child component that is wrapped by a react-dimensions higher-order component to automatically resize the map viewer according to parent element/window size changes (so we can hook in and call updateSize on the internal ol.Map to prevent image distortion)

The major deviation from my original thoughts is that in the above HTML code, there is no references to React or having to mount the viewer component using the React APIs. I've come to realize that this would be cumbersome as your ultimate entry point will always be HTML/JavaScript, way outside the zone of React/TypeScript/(T|J)SX, and although this is a "component", our final bundle is ultimately more an application instead of a library. So having a wrapper (which is what MapGuide.Application does) to do all of this leg work simplifies the whole setup process.

That's not to say we can't componentize this map viewer (as a library/package you can npm install), but my primary goal is to have an experience similar to our existing viewer offerings, which are applications whose features and behaviour is driven by configuration files. So in that respect, including a script tag and writing some JavaScript to point to your configuration and kickstart the application is much more easier to comprehend than installing node, npm, Webpack and messing with frontend development stacks.

So far so good. We're just beginning to scratch the surface.
Categories: OSGeo Planet
Syndicate content