Welcome to Planet OSGeo

May 19, 2024

May 18, 2024

May 17, 2024

May 16, 2024

The most fundamental and widely-used operations in the JTS Topology Suite are the ones that evaluate topological relationships between geometries.  JTS implements the Dimensionally-Extended 9 Intersection Model (DE-9IM), as defined in the OGC Simple Features specification, in the RelateOp API.  

DE-9IM matrix for overlapping polygons

The RelateOp algorithm was the very first one implemented during the initial JTS development, over 20 years ago.  At that time it was an appealing idea to implement a general-purpose topology framework (the GeometryGraph package), and use it to support topological predicates, overlay, and buffering.  However, some disadvantages of this approach have become evident over time:

  • the need to create a topological graph structure limits the ability to improve performance. This has led to the implementation of PreparedGeometry - but that adds further complexity to the codebase, and supports only a limited set of predicates.
  • a large number of code dependencies make it hard to fix problems and improve semantics
  • constructing a full topology graph increases exposure to geometric robustness errors
  • GeometryCollections are not supported (initially because the OGC did not define the semantics for this, and now because adding this capability is difficult)
There's an extensive list of issues relating (ha!) to RelateOp here.

The importance of this functionality is especially significant since the same algorithm is implemented in GEOS.  That codebase is used to evaluate spatial queries in popular spatial APIs such as Shapely and R-sf, and numerous systems such as PostGISDuckDBSpatialLiteQGIS, and GDAL (to name just a few).  It would not be surprising to learn that the RelateOp algorithm is executed billions of times per day across the world's CPUs.

During the subsequent years of working on JTS I realized that there was a better way to evaluate topological relationships.  It would required a ground-up rewrite, but would avoid the shortcomings of RelateOp and provide better performance and a more tractable codebase.  Thanks to my employer Crunchy Data I have finally been able to make this idea a reality.  Soon JTS will provide a new algorithm for topological relationships called RelateNG.

Key Features of RelateNG

The RelateNG algorithm incorporates a broad spectrum of improvements over RelateOp in the areas of functionality, robustness, and performance.  It provides the following features:

  • Efficient short-circuited evaluation of topological predicates (including matching custom DE-9IM matrix patterns)
  • Optimized repeated evaluation of predicates against a single geometry via cached spatial indexes (AKA "prepared mode")
  • Robust computation (only point-local geometry topology is computed, so invalid topology does not cause failures)
  • GeometryCollection inputs containing mixed types and overlapping polygons are supported, using union semantics.
  • Zero-length LineStrings are treated as being topologically identical to Points.
  • Support for BoundaryNodeRules.

Using the RelateNG API

The main entry point is the RelateNG class.  It supports evaluating topological relationships in three different ways:

  • Evaluating a standard OGC named boolean binary predicate, specified via a TopologyPredicate instance.  Standard predicates are obtained from the RelatePredicate factory functions intersects, contains, overlaps, etc.
  • Testing an arbitrary DE-9IM relationship by matching an intersection matrix pattern (e.g. "T**FF*FF*", which is the pattern for a relation called Contains Properly).
  • Computing the full value of a DE-9IM IntersectionMatrix.
RelateNG operations can be executed in two modes: stateless and prepared.  Stateless mode is provided by the relate static functions, which accept two input geometries and evaluate a predicate.  For prepared mode, a single geometry is provided to the prepare functions to create a RelateNG instance.  Then the evaluate instance methods can be used to evaluate spatial predicates on a sequence of geometries. The instance creates spatial indexes (lazily) to make topology evaluation much more efficient.  Note that even stateless mode can be significantly faster than the current implementation, due to short-circuiting and other internal heuristics.

Here is an example of matching an intersection matrix pattern, in stateless mode:
boolean isMatched = RelateNG.relate(geomA, geomB, "T**FF*FF*");
Here is an example of setting up a geometry in prepared mode, and evaluating a named predicate on it:
RelateNG rng = RelateNG.prepare(geomA);
for (Geometry geomB : geomSet) {
    boolean predValue = rng.evaluate(geomB, RelatePredicate.intersects());
}

Rolling It Out

It's exciting to launch a major improvement on such a core piece of spatial functionality.  The Crunchy spatial team will get busy on porting this algorithm to GEOS.  From there it should get extensive usage in downstream projects.  We're looking forward to hearing feedback from our own PostGIS clients as well as other users.  We're always happy to be able to reduce query times and equally importantly, carbon footprints.  

In further blog posts I'll describe the RelateNG algorithm design and provide some examples of performance metrics.

Future Ideas

The RelateNG implementation provides an excellent foundation to build out some interesting extensions to the fundamental DE-9IM concept.

Extended Patterns

The current DE-9IM pattern language is quite limited.  In fact, it's not even powerful enough to express the standard named predicates.  It could be improved by adding features like:

  • disjunctive combinations of patterns.  For example, touches is defined by  "FT******* | F**T***** | F***T****"
  • dimension guards to specify which dimensions a pattern applies to.  For example, overlaps is defined by "[0,2] T*T***T** | [1] 1*T***T**"
  • while we're at it, might as well support dotted notation and spaces for readability; e.g. "FT*.***.***"

Approximate Spatial Relationships

A requirement that comes up occasionally is to compute approximate spatial relationships between imprecise data using a distance tolerance.  This is useful in datasets where features don't match exactly due to data inaccuracy. Because RelateNG creates topology only in the neighbourhood of specified points, it should be possible to specify the neighbourhood size using a distance tolerance.  Essentially, vertices and line intersections will "snap together" to determine the effective topology.  Currently within-distance queries are often used to compute "approximate intersects".  Adding a distance tolerance capability to RelateNG will allow other kinds of approximate relationships to be evaluated as well.

Further Performance Optimizations?

A challenge with implementing algorithms over a wide variety of spatial types and use cases is how to provide general-purpose code which matches (or exceeds) the efficiency of more targeted implementations.  RelateNG analyzes the input geometries and the predicate under evaluation to tune strategies to reduce the amount of work needed to evaluate the DE-9IM.  It may be that profiling specific use cases reveals further hotspots in the code which can be improved by additional optimizations.  

Curve Support

GEOS has recently added support for representing geometries with curves.  The RelateNG design is modular enough that it should be possible to extend it to allow evaluating relationships for geometries with curves.  


by Dr JTS (noreply@blogger.com) at May 16, 2024 09:41 PM

A consulta geoespacial, ou SQL espacial, está revolucionando a maneira como conduzimos operações de Sistemas de Informações Geográficas (GIS). Ao aproveitar funções e recursos espaciais em bancos de dados SQL, podemos analisar e obter insights valiosos de dados espaciais de maneira transparente.

Uma das principais vantagens do SQL espacial é a sua capacidade de encontrar relações entre geometrias. Seja para determinar proximidade, sobreposição ou contenção, o SQL espacial nos permite desbloquear conexões significativas em conjuntos de dados espaciais. Esta funcionalidade é crucial para diversas aplicações, desde planejamento urbano até monitoramento ambiental e muito mais.

Além disso, a integração de SQL espacial em processos de back-end enriquece nosso código com poderosos recursos analíticos. Ao aplicar a análise espacial diretamente em nossas consultas de banco de dados, simplificamos os fluxos de trabalho e aumentamos a eficiência da tomada de decisões baseada em dados.

No mundo atual orientado por dados, dominar o SQL espacial abre portas para um mundo de possibilidades em GIS e muito mais. Você está pronto para aproveitar todo o potencial da análise de dados espaciais?

Gostou desse post? Conte nos comentários 👇

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 16, 2024 12:00 PM

May 15, 2024

Walter Schwartz sentg us this, he said “Northwestern Mexico, where Baja, Sonora, and Sinaloa provinces are found, has been home to a *lot* of missions since at least 1591. I counted around 100 established in the 1600s and 1700s alone. This picto-map is in the pictured church, La Mision de San Ignacio in the town of San Ignacio, Baja California Sur, Mexico. 

MapsintheWild Missions of Mexico

by Steven at May 15, 2024 09:00 AM

May 14, 2024

Um software capaz de identificar, a partir de imagens aéreas, caixas d’água sobre telhados ou lajes e piscinas em áreas abertas foi desenvolvido por pesquisadores brasileiros com o auxílio de ferramentas de Inteligência Artificial. A proposta é usar esse tipo de imagem como indicador de zonas especialmente vulneráveis a infestações do mosquito Aedes aegypti, transmissor de doenças como dengue, zika e chikungunya. Além disso, a estratégia desponta como potencial alternativa para um mapeamento socioeconômico dinâmico das cidades – um ganho para diferentes políticas públicas.

A pesquisa, apoiada pela Fapesp, foi conduzida por profissionais da USP, da Universidade Federal de Minas Gerais (UFMG) e da Superintendência de Controle de Endemias (Sucen) da Secretaria de Estado da Saúde de São Paulo.

Entre outras coisas, o grupo almeja incorporar outros elementos para serem detectados nas imagens e quantificar as taxas reais de infestação do Aedes aegypti em uma dada região para refinar e validar o modelo. “Nós esperamos criar um fluxograma que possa ser aplicado em diferentes cidades para encontrar áreas de risco sem a necessidade de visitas domiciliares, prática que gasta muito tempo e dinheiro público”.

Apesar de as fotos aéreas de Campinas terem sido obtidas com um drone, espera-se que, no futuro, a estratégia testada nessa pesquisa recorra apenas às imagens de satélite.

No estudo em Belo Horizonte, as imagens de satélite foram empregadas com sucesso – elas precisam de alta resolução para que o computador consiga identificar os padrões.

Embora esse tipo de metodologia pareça custoso, ele gera uma potencial economia ao dispensar a necessidade de visitas presenciais para mapear, casa por casa, áreas suscetíveis à dengue. Em vez disso, os agentes de saúde aproveitariam as informações obtidas remotamente – e processadas com a Inteligência Artificial – para se dirigir aos locais prioritários com mais assertividade.

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 14, 2024 12:00 PM

Derick Rethams spotted this beer mat in his local pub in Maida Vale. For those of you who are interested

Southern Hemisphere IPA Northern Monk Collaboration. Motueka and Wai-Iti hops bring notes of lime, stone fruits and citrus, alongside a lofty dry hopping of Eclipse and Nectaron for a vibrant mix of orange, tropical fruit and a touch of pine to finish.

https://www.eebriatrade.com/products/beer/double-barrelled-brewery/63341-aerial-ping-pong#more-items

MapsintheWild Aerial Ping Pong

by Steven at May 14, 2024 09:00 AM

May 13, 2024

At OPENGIS.CH, we’ve been working lately on improving the DXF Export QGIS functionality for the upcoming release 3.38. In the meantime, we’ve also added nice UX enhancements for making it easier and much more powerful to use!

Let’s see a short review.

DXF Export app dialog and processing algorithm harmonized

You can use either the app dialog or the processing algorithm, both of them offer you equivalent functionality. They are now completely harmonized!

Export settings can now be exported to an XML file

You can now have multiple settings per project available in XML, making it possible to reuse them in your workflows or share them with colleagues.

Load DXF settings from XML.

All settings are now well remembered between dialog sessions

QGIS users told us there were some dialog options that were not remembered between QGIS sessions and had to be reconfigured each time. That’s no longer the case, making it easier to reuse previous choices.

“Output layer attribute” column is now always visible in the DXF Export layer tree

We’ve made sure that you won’t miss it anymore.

DXF Export, output layer attribute

Possibility to export only the current map selection

Filter features to be exported via layer selection, and even combine this filter with the existing map extent one.

DXF Export algorithm, use only selected features

Empty layers are no longer exported to DXF

When applying spatial filters like feature selection and map extent, you might end up with empty layers to be exported. Well, those won’t be exported anymore, producing cleaner DXF output files for you.

Possibility to override the export name of individual layers

It’s often the case where your layer names are not clean and tidy to be displayed. From now on, you can easily specify how your output DXF layers should be named, without altering your original project layers.

Override output layer names for DXF export.

We’ve also fixed some minor UX bugs and annoyances that were present when exporting layers to DXF format, so that we can enjoy using it. Happy DXF exporting!

We would like to thank the Swiss QGIS user group for giving us the possibility to improve the important DXF part of QGIS 🚀🚀🚀

by Germán Carrillo at May 13, 2024 04:48 AM

Scanxiety.

This is where I am right now. Scanxiety.

Each stage of the cancer experience is marked by a particular set of tests, of scans.

I actually managed to get through my first set of scans surprisingly calmly. After getting diagnosed (“there’s some cancer in you”), they send you for “staging”, which is an MRI and CT scan.

These scans both involve large, Star Trek seeming machines, which make amazing noises, and in the case of the CT machine I was put through was decorated with colorful LED lights by the manufacturer (because it didn’t look whizzy enough to start with?).

MRI

I kind of internalized the initial “broad-brush” staging my GI gave me, which was that it was a tumor caught early so I would be early stage, so I didn’t worry. And it turned out, that was a good thing, since the scans didn’t contradict that story, and I didn’t worry.

The CT scan, though, did turn up a spot on my hip bone. “Oh, that might be a bone cancer, but it’s probably not.” Might be a bone cancer?!?!?

How do you figure out if you have “a bone cancer, but it’s probably not”? Another cool scan, a nuclear scan, involving being injected with radioactive dye (frankly, the coolest scan I have had so far) and run through another futuristic machine.

Bone Scan

This time, I really sweated out the week between the scan being done and the radiology coming back. And… not bone cancer, as predicted. But a really tense week.

And now I’m in another of those periods. The result of my major surgery is twofold: the piece of me that hosted my original tumor is now no longer inside of me; and, the lymph nodes surrounding that piece are also outside of me.

They are both in the hands of a pathologist, who is going to tell me if there is cancer in the lymph nodes, and thus if I need even more super unpleasant attention from the medical system in the form of several courses of chemotherapy.

The potential long term side effects of the chemotherapy drugs used for colorectal cancers include permanent “peripheral neuropathy”, AKA numbness in the fingers and toes. Which could put a real crimp in my climbing and piano hobbies.

Climb

So as we get closer to getting that report, I am experiencing more and more scanxiety.

If I escape chemo, I will instead join the cohort of “no evidence of disease” (NED) patients. Not quite cured, but on a regular diet of blood work, scans, and colonoscopy, each one of which will involve another trip to scanxiety town. Because “it has come back” starts as a pretty decent probability, and takes several years to diminish to something safely unlikely.

Yet another way that cancer is a psychological experience as well as a physical one.

Talk to you again soon, inshalla.

May 13, 2024 12:00 AM

May 12, 2024

May 11, 2024

As a self-taught PyQGIS developer, one of my main hurdles has always been to prepare the development environment for PyQGIS. An environment that allow me to run PyQGIS scripts, helps me code faster by providing PyQGIS highlighting and autocompletion, enables me to debug my plugins and scripts as they run, etc… I have been a user (and even a… “cof cof”… maintainer) of the QGIS packages for conda provided by the conda-forge community.

May 11, 2024 01:32 AM

May 10, 2024

GeoSpatial Techno is a startup focused on geospatial information that is providing e-learning courses to enhance the knowledge of geospatial information users, students, and other startups. The main approach of this startup is providing quality, valid specialized training in the field of geospatial information.

( YouTube | LinkedIn | Facebook | X )


How to create Tile Layers with GeoServer

In this session, we will explore “How to create Tile Layers with GeoServer”. If you want to access the complete tutorial, simply click on the link.

Introduction

GeoWebCache is a tiling server that acts as a proxy between a map client and map server, caching tiles to save processing time. It is integrated with GeoServer and can significantly improve the responsiveness and reliability of the server. These settings can be accessed from the left side of the screen under the Tile Caching heading. They include:

  • Tile Layers: This section lists all cached layers for review and parameter modification.
  • Caching Defaults: Caching Defaults is the entry point for these.
  • Gridsets: This option allows you to create new tiling schemes or modify the existing ones.
  • Disk Quota: The Disk Quota and BlobStores options allow you to set predefined amounts of disk space for each layer.

Tile Layers

The tile layers menu shows a listing of all of the layers known to the integrated GeoWebCache. It is similar to the Layer Preview for GeoWebCache, with many of the same options that you can review the status and the main parameters for each layer. For each layer cached by GeoWebCache, the following information is available:

  • Type and Layer Name: These columns display the type of layer and the name value of each.
  • Disk Quota: The maximum amount of disk space that can be used for this layer. This feature cannot be configured in GeoServer as in the GeoWebCache standalone version, so you can only see the N/A value here.
  • Disk Used: The current disk space being used by tiles for this particular layer. This counter will only be updated if disk quotas are enabled.
  • Enabled: Indicates whether tile caching is enabled for this layer. It is possible to have a layer definition here but not to have tile caching enabled (set in the layer properties).
  • Preview: Similar to Layer Preview, this will generate a simple OpenLayers application populated with tiles from one of the available gridset/image format combinations.
  • Actions: Actions have two options, Seed/Truncate and Empty link.
  • Seed/Truncate: Opens the GeoWebCache page for automatically seeding and truncating the tile cache. Use this if you want to pre-populate some of your cache.
  • Empty: This removes all saved tiles from the cache and is identical to a full truncate operation for the layer.
  • Empty All: This link lets you clear the entire cache for all layers, gridsets and filter parameter combinations. Note that, This will truncate all layers in GeoWebCache.

Adding or removing cached layers

We have used the GeoWebCache to store tiles generated by user requests, but you can also precompute tiles for a layer to avoid delays for some users. The process of pre-computing tiles is called seeding. This section will help you understand how it works:

  • Navigate to the Tile Caching > Tile Layers page and click on the Seed/Truncate link for your desired layer.
  • Scroll to the Create a new task section and set the parameters for seeding. For the number of parallel processes (threads) that will request maps from GeoServer, it’s best to select a number of threads equal to “50” percent of the available cores on your server.
  • From the Operation Type you have a few options: Seed, Reseed and Truncate. Select Seed to generate missing tiles or Reseed to regenerate all tiles. Truncate allows you to select specific zoom levels to remove, unlike the Clear operation, which removes all tiles.
  • You must select a gridset and an image format to seed. If you want to precompute cache for more than one gridset and/or image format, you can start another operation immediately after this one.
  • To start seeding, select a subset of the gridset, specify a level range and an area. If no specific area is needed, leave the Bounding Box empty. Select Zoom Start and Zoom Stop, then press Submit to start seeding.
  • The web interface displays a list of running tasks, which can be filtered by layer and stopped if necessary. Refreshing the list shows the number of tiles completed, time elapsed, and time remaining. Seeding multiple layers can be resource intensive and time consuming.
  • When tasks are completed, check for an empty list. Return to the Tile Layers page to see a significant increase in the amount of disk space allocated for your layer’s tiles. Seeding your layers can greatly improve performance, as map requests will now hit the cache in the precomputed layer range. Expect a 10- to 90-fold increase in performance.

Tiles can be generated by GeoWebCache in two ways: on-demand during map viewing, or by seeding in advance. On-demand caching speeds up subsequent views, but reduces the quality of the user experience. Seeding improves the user experience but is time and disk consuming. Typically, a combination of both methods is used, with popular areas or zoom levels seeded and less frequently viewed tiles left uncached.

Disk Quota

The Disk Quotas page lets you manage disk usage for cached tiles and set a global disk quota. By default, disk usage for cached tiles is unlimited, but setting a quota can prevent disk capacity issues, especially with Direct WMS integration. This is important for serving large areas with terabytes of tile cache on disk. The following information is displayed on the Disk Quota page:

  • Enable disk quota: When enabled, the disk quota will be set according to the options listed below. The setting is disabled by default.
  • Disk quota check frequency: This setting determines how often the cache is polled for any overage. Smaller values (more frequent polling) will slightly increase disk activity, but larger values (less frequent polling) may cause the disk quota to be temporarily exceeded. The default is 10 seconds.
  • Maximum tile cache size: The maximum size for the cache. When this value is exceeded and the cache is polled, tiles will be removed according to the policy. Note that the unit options are MegaBytes, GigaBytes and TeraBytes. The default is 500 MiB.
  • Tile removal policy: When disk quota is exceeded, choose between Least Frequently Used or Least Recently Used policy to delete tiles based on access frequency or date last accessed. Optimal configuration varies based on data and server usage.
  • Disk quota store type: GeoServer uses the H2 database for disk quotas by default, with the option to configure an external database. You can choose an in-process database with a local H2 database or an external database with H2, PostgreSQL or Oracle. Don’t forget to press Submit after making changes.

In this session, we explored “How to create Tile Layers with GeoServer”. If you want to access the complete tutorial, simply click on the link.

by Nima Ghasemloo at May 10, 2024 12:00 AM

May 09, 2024

Hoy vamos a hablar acerca de un proyecto que aborda una de las temáticas que actualmente, requiere de los esfuerzos conjuntos entre autoridades, la academia y la sociedad: la violencia en contra de las mujeres y niñas.

Al ser considerada como unas de las violaciones a los derechos humanos más extendidas en el mundo ha llevado a la generación de múltiples propuestas, cuyo propósito es generar alternativas (desde diversos enfoques) para brindar soluciones en este tema. Tal es el caso del proyecto Criteria Taronja, desarrollado por la Universidad Politécnica de Valencia (UPV) por el Grupo de Investigación en Tecnologías Geoespaciales (GITG) de dicha Universidad, que a través de un modelo de datos específico realiza el geoprocesamiento de diversas variables socioeconómicas sobre una unidad territorial homogénea, identificando y visualizando aquellos espacios públicos o zonas potencialmente inseguras para las mujeres. Este modelo se ha complementado con un módulo de cartografía participativa, que con un enfoque de ciencia ciudadana geográfica integra fuentes no tradicionales de información para la generación de conjuntos de datos con una escala de desagregación y análisis con mayor detalle. En este punto es importante resaltar la importancia de la participación de las ciudadanas usuarias del espacio público, ya que sus aportaciones permitirán enriquecer el modelo a través de cada sitio que reporten.

El módulo de cartografía de cartografía participativa se presenta como un aplicativo web el cual puede usarse desde el ordenador o computadora, incluso también desde dispositivos móviles. Al momento de ingresar, un formulario permite especificar nuestro perfil como participantes, proporcionando la garantía del uso de los datos a través de la política de privacidad. Una vez hecho lo anterior, a través de un marcador de posición, podemos especificar los puntos sobre la ciudad en donde tenemos percepción de inseguridad, o bien hayamos tenido o presenciado situaciones de acoso callejero.

En el siguiente videotutorial podemos ver la funcionalidad del aplicativo.

Actualmente el proyecto de cartografía participativa se desarrolla en dos ciudades: Toluca (México) y Valencia (España). En el caso de Toluca, la convocatoria a ciudadanos para el levantamiento de datos ha finalizado. En Valencia, se están llevando acciones para la colaboración de los habitantes de esta municipalidad, a través de redes sociales y otros medios de contacto.

Hacemos una cordial invitación para participar y compartir este proyecto, que nos plantea la posibilidad de contribuir en actividades de investigación y la implementación de soluciones a problemas específicos, mejorando las interacciones entre la sociedad y la ciencia.

Enlace a la plataforma de cartografía participativa: https://geoinformatica.mx/criteriaTaronja/valencia

Fdo. Sandra Lucía Hernández Zetina

Facultad de Geografía
Universidad Autónoma del Estado de México

by Alvaro at May 09, 2024 03:24 PM

Você sabe quais são os 3 conceitos e casos de uso essenciais do Geopandas?

Imagine que você está trabalhando em um projeto GIS e deseja obter insights e trabalhar com seus dados GIS e não existem funções específicas que você precisaria em um software GIS conhecido.

Geopandas é extremamente útil para criar suas próprias funcionalidades GIS com sua excelente documentação e funções simples de usar.

👉 Aqui estão os três conceitos/casos de uso mais importantes no GeoPandas:

📍 Operações de dados geoespaciais: Permite realizar operações espaciais como merge, buffer e análise de sobreposição de camadas.
📍 Visualização de dados espaciais: Você pode traçar e personalizar mapas, adicionar camadas, cores e rótulos para tornar seus dados compreensíveis.
📍 Cálculos Geoespaciais: Você pode fazer cálculos para obter cada vez mais insights de seus dados geoespaciais.

👉 Veja um exemplo de cada conceito/caso de uso:

📍 Crie um GeoDataFrame, por exemplo, representando cidades com dados populacionais e suas localizações – use geopandas.read_file(“path/to/your/file.extension”) para permitir que o GeoPandas crie essa estrutura automaticamente para você.
📍 Visualize as cidades num mapa e pinte-as com base na população.
📍 Execute uma operação espacial por buffer – o código irá armazenar em buffer a geometria de cada cidade e armazená-la em ‘Buffered_Geometry’.
📍 Calcule a densidade populacional de cada cidade dividindo a população pela área aproximada da cidade.

Essas operações mostram um pouco do poder do GeoPandas na visualização e análise de dados geoespaciais.

Você sabia desse potencial todo do GeoPandas? Nos conte nos comentários 👇

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 09, 2024 12:00 PM

May 08, 2024

Thanks to our generous donors and sustaining members, we are in the wonderful position to be able to further extend our 2024 Grant Programme and to fund an additional project that came in very close 6th in the voting results:

On behalf of the QGIS.ORG project, I would like to thank the Danish User Group who stepped up and increased their support to Flagship level, joining or recurring Flagship member Felt and all our wonderful sustaining members and donors.

The next proposals on the wait list is:

with € 4,500 missing to be be funded to make it possible to store encrypted credentials in a PostgreSQL database. This improvement targets QGIS Server, particularly multi-instance setups using kubernetes, with potential QGIS Desktop improvements down the road.

If you want to help make this improvement a reality, please get in touch.

Of course, the same applies if you want to help fund any of the remaining proposals:

by underdark at May 08, 2024 06:12 PM

Tem algumas áreas que é difícil, no primeiro momento, visualizar como o geoprocessamento pode ajudar. E hoje falaremos de uma dessas áreas, que é a jurídica.

👉 Mas como o geoprocesssamento pode fazer a diferença e agregar valor nessa área? Separamos alguns pontos interessantes, que iremos comentar abaixo:

📍 Através de uma plataforma WebGIS, por exemplo, onde você pode ter a representação visual de dados complexos trazendo assim uma melhor compreensão em casos legais.

📍 Apoio na idenficação de padrões e tendências geográficas em investigações e processos judiciais.

📍 Apoio na identificação da localização de cenas, propriedades e outros elementos.

📍 Em processos imobiliários, fornecendo uma visão clara da topografia e distribuição geográfica no gerenciamento de terras, e disputa de limites, por exemplo.

Você já sabia disso? Sabe de mais alguma aplicação que não citamos, deixe nos comentários 👇

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 08, 2024 12:00 PM

I have a profoundly embarassing cancer. Say it with me “rectal cancer”. “Rectal cancer”.

Why is it embarassing?

Poop!?! Maybe we are all still six, somewhere deep inside.

When Ryan Reynolds got a colonoscopy on camera, to raise awareness of colorectal cancer screening, part of the frisson of the whole thing was that yes, somehow having this procedure done is really embarassing.

So, watch the video, it’s really nothing but an ordinary medical procedure that could very well save your life. And Ryan Reynolds is charming.

Cows

Meanwhile, colo-rectal cancers remain tough to talk about, because frankly the colonoscopy is the least of it.

Not having control of your bowels is, well, really embarassing in our culture. What do people say about elderly presidential candidates they hate? They call them incontinent. They intimate that they wear adult diapers (gasp!).

Do you know who else gets to wear adult diapers? Colorectal cancer patients. We get our insides man-handled, irradiated and chopped up, and the results are not great for bowel control. It happens if you’re 55, it happens if you’re 35. It’s normal, it’s usually temporary, it’s what happens when you insult a GI tract badly enough.

Another rite of passage in treatment is the ostomy. Stage III rectal cancer treatment usually involves a temporary ostomy, after radio-chemotherapy during the resection of the part of the rectum that holds the tumor. Patients with a low (near the anus) tumor location will sometimes require a permanent ostomy, because the tumor cannot be removed without damaging the anus.

When I was diagnosed, I was initially terrified of the ostomy. “The bag.”

After researching the different treatments, I got a lot less terrified, since the side effects of some non-bag outcomes in terms of quality of life can be pretty terrible. Meanwhile folks with ostomies are out hiking, biking, and swimming.

If this talk is all a little uncomfortable, may I recommend a colonoscopy?

And after that, a big meal and some poooooping! Poop! Poop! Poop!

I’m in a pooping mood because my surgery (2 weeks ago now) has left me, not incontinent, but I guess “disordered” is a better word. You know how it feels to really need to take a dump? Imagine feeling that 12 hours a day, even when you don’t actually have anything to dump.

By most measures I think I am ahead of the median patient in recovery from LAR surgery, but unfortunately the recovery time for things like bowel regularity and “normalcy” (the “new normal” will always be somewhat worse than the “old normal”) is measured in months, not days, so I am a little impatient to improve more, and faster.

Talk to you again soon, inshalla.

Everyone Poops

May 08, 2024 12:00 AM

May 07, 2024

Se você trabalha com geoprocessamento, certamente está acostumado a lidar com a informação em diversos formatos, como KML, Shapefile, GeoPackage, entre outros.

Nesse post vamos falar um pouco sobre um outros formato que é o GeoJSON, ele se baseia no JSON (JavaScript Object Notation), muito usado como formato de troca de dados entre sistemas e APIs.

O GeoJSON é, portanto, um JSON que traz dados que representam recursos geográficos simples, como pontos, multipontos, linhas, multilinhas, polígonos, multipolígonos, geometrycollections e seus atributos não espaciais.

Por convenção, o GeoJSON sempre deve se basear no datum WGS84, com unidades de latitude e longitude em graus decimais. Ou seja, não se deve utilizar coordenadas em graus minutos e segundos.

Para finalizar vamos falar da precisão, que pode ser definida em 8, 10 ou até mesmo 15 casas decimais, mas na prática pode ser um desperdício de bytes, pois em linhas gerais uma precisão de 6 casas decimais equivale a cerca de 10cm, o que é mais do que suficiente para 99% dos casos.

Você já conhecia o GeoJSON? Deixe nos comentários 👇

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 07, 2024 12:00 PM

May 06, 2024

O EPSG é um registro público de datums geodésicos, sistemas de referência espacial, elipsóides terrestres, transformações de coordenadas e unidades de medida relacionadas. A cada entidade é atribuído um código EPSG entre 1024 e 32767, junto com uma representação padrão de texto no formato WKT.

Foi criado por Jean-Patrick Girbig membro do European Petroleum Survey Group (EPSG) em 1985 para padronizar, melhorar e compartilhar dados espaciais entre os membros do grupo, tornando-se público em 1993.

A maioria dos sistemas de informação geográfica (GIS) e bibliotecas GIS usam códigos EPSG como identificadores de sistema de referência espacial (SRIDs) e dados de definição EPSG para identificar sistemas de referência de coordenadas, projeções e realizar transformações entre esses sistemas.

👉 Códigos EPSG comuns:

📍 EPSG:4326 – WGS 84, sistema de coordenadas de latitude/longitude baseado no centro de massa da Terra, utilizado pelo Sistema de Posicionamento Global entre outros;
📍 EPSG:3857 – Projeção Web Mercator, usada para exibição por muitas ferramentas de mapeamento baseadas na web, incluindo Google Maps e OpenStreetMap;
📍 EPSG:4674 – SIRGAS 2000, sistema de coordenadas de latitude/longitude utilizando oficialmente no Brasil.

Você sabia essas informações sobre o EPSG? Conte nos comentários 👇

Fonte: webgis.tech
Instagram: https://instagram.com/webgis.tech
LinkedIn: https://www.linkedin.com/company/webgis-tech

by Fernando Quadro at May 06, 2024 12:00 PM

May 05, 2024

“Anything that’s human is mentionable, and anything that is mentionable can be more manageable. When we can talk about our feelings, they become less overwhelming, less upsetting, and less scary. The people we trust with that important talk can help us know that we are not alone.”
– Fred Rogers

When I found out I had rectal cancer, I hit the internet hard and immediately found the ColonTown community of communities. It has turned out to be simultaneously the most reassuring and the most anxiety producing place on the internet for me.

On the reassuring side, even though colorectal cancer is third most prevalent cancer world-wide, it is not widely talked about, so the community was reassuring: there are other people out there going through this; many are already through it.

There are also a lot of aspects of the rectal cancer treatment process that the medical system seems ill-equiped to support. I was sent home from major surgery with very little guidance in hand about expected recovery progression, or diet. Fortunately the community was a great resource for this information.

On the anxiety producing side, let me start with this meme.

Survivorship bias

The population of a community like ColonTown is necessarily going to bias towards people who are currently in treatment and the population of people currently in treatment will bias toward folks whose treatment is not necessarily getting them cured.

Survivorship bias in this case manifests in the survivors slowly draining out of the community, leaving behind mostly folks still actively in treatment.

There are a lot of people in Stage IV, keeping fucking going, who have harrowing tales. And there are people who were Stage II (like me) who against the odds progressed to Stage IV. It happens. And the thought that I could follow that path too, is frankly terrifying.

Anyways, as a result, reading posts from this wonderful supportive online community can sometimes throw me into a spiral of anxiety, because oh my god I have this terrible deadly thing. Which I do. But also probably I am going to be fine.

Probably?

Talk to you again soon, inshalla.

May 05, 2024 12:00 AM

May 04, 2024

Trajectools continues growing. Lately, we have started expanding towards public transport analysis. The algorithms available through the current Trajectools development version are courtesy of the gtfs_functions library.

There are a couple of existing plugins that deal with GTFS. However, in my experience, they either don’t integrate with Processing and/or don’t provide the functions I was expecting.

So far, we have two GTFS algorithms to cover essential public transport analysis needs:

The “Extract shapes” algorithm gives us the public transport routes:

The “Extract segments” algorithm has one more options. In addition to extracting the segments between public transport stops, it can also enrich the segments with the scheduled vehicle speeds:

Here you can see the scheduled speeds:

To show the stops, we can put marker line markers on the segment start and end locations:

The segments contain route information and stop names, so these can be extracted and used for labeling as well:

If you want to reproduce the above examples, grab the open Vorarlberg public transport schedule GTFS.

These developments are supported by the Emeralds Horizon Europe project.

by underdark at May 04, 2024 05:50 PM

Marc-Tobias has been visiting Hiroshima in Japan and shred these awesome street maps that he spotted. He said:

“As usual in Japan streets are unnamed. To add confusion many
shops are stacked and not next to each other (‘3F’ means 3rd floor)

I saw only those two sheets. Unclear who made them, which purpose,
if they’re updated. I don’t think people use them for navigation.”

Despite the weirdness (to us at least) of the Japanese addressing system people manage to find their way, mail and goods get delivered. Clearly it works.

MapsintheWild Magnificently Confused in Hiroshima

by Steven at May 04, 2024 10:00 AM

With the QGIS Grant Programme 2023, we were able to support four proposals in the first round and an additional two proposals through the second round that are aimed to improve the QGIS project, including software, infrastructure, and documentation. The following reports summarize the work performed in the proposals. 

  1. QGIS Bug Tracker cleanup #266Report
    While checking existing bugreports we have identified and closed ~291 tickets, among them:

    – 162 bugreports and feature requests which were already fixed or implemented
    – 29 bugreports and feature requests which are invalid (data issues, wrong use of functionality, etc)
    – 57 duplicate bugreports and feature requests
    – 5 won’t fix bugreports
    – 5 bugreports were converted to feature requests
    – 33 tickets were closed (does not contain steps to reproduce, test data and no feedback was provided within several month)

    Additionally we ensured that all tickets has correct tags assigned to to make them easier to find.

  2. Porting to C++ and harmonization of Processing algorithms #271Report
    The work has been completed with several pull requests:
    Port various geometry algorithms to C++ 
    Port Align Rasters tool to Processing 
    Port Raster Calculator algorithm to C++ 
    Port XYZ Tiles algorithms to C++ 

    Existing Processing algorithms Voronoi Polygons and Delaunay Triangulation have been ported to C++ and now use GEOS instead of the unmaintained Python module. Both algorithms have been given additional parameters to disable post-processing step (adding attributes of input points for Voronoi polygons and adding feature IDs of the 3 points forming a triangle for Delaunay triangulation) and thus improve algorithm performance. The Concave Hull algorithm has been ported to C++ and uses GEOS for QGIS builds with GEOS >= 3.11, while a port of the existing implementation based on Delaunay triangulation is used for builds with older GEOS versions.

    Two algorithms for generating XYZ tiles (directory and MBTiles variants) have been ported to C++ using a safer and cleaner multi-threading approach.

    The Align Rasters tool (available from the Raster → Align Rasters menu), which was not exposed to Processing, has been removed and a new Processing algorithm with the same functionality has been added. To make this possible, we have added a new QgsProcessingParameterAlignRasterLayers parameter type to Processing, which provides an adapter to QList<QgsAlignRaster::Item>. The “Raster” menu item now opens the Processing algorithm. We have also added a separate, simplified modeler-only version of the algorithm that accepts a single raster to align and returns a single raster layer output.

    The existing Raster Calculator algorithm has been ported to C++. The algorithm now has two variants: a toolbox version that works the same way as before, and a modeler version that uses the same approach to input naming as the GDAL raster calculator (in the formula, you should use letters instead of layer names). The old algorithm implementation has been retained for compatibility with existing models and will be removed in QGIS 4. We have also added a separate raster calculator algorithm for creating virtual raster layers.
  3. Improve test result handling on QGIS CI #268 Report
    While the original proposal was explicitly stated to be a research project with no guarantees of success, the end result is predominantly a success (with some limitations!). You can see the new failure handling in action in this PR

    What we have now is that any tests which fail a rendering comparison will write a descriptive comment to the PR, as shown in the above link. The comment details which render tests failed, where they are in the code, and includes some helpful pointers to downloading the full test report and the QGIS developer documentation.

    Originally, I hoped to link directly to the full test report or include it as an attachment to the comment. Unfortunately this is NOT possible given the current Github API. There’s a bunch of notes I’ve added to the initial comment which link to the limitations / feature requests on Github’s side, so we can monitor the situation and further improve the reports if/when Github add this functionality.

    As well as the above described improvements on the CI side, I’ve also implemented lots of improvements in running the tests locally and how the render test reports are generated and presented to developers!

  4. Unification of all geometric and topological verification methods #236Report
    Our efforts are primarily focused on two main aspects of the proposal:
    1. Unification of Algorithms:
    This component entails both refactoring existing algorithms and introducing new ones. While the refactoring phase did not reach the level of ambition initially envisioned, considerable cleanup has been achieved, particularly in eliminating code duplications, as evidenced by the PR submitted [1] [2] [3]. Other deduplications could be considered, but the main issue is not to break the small optimizations or adaptations to the logic of the existing code.

    2. Integration of Geometry Checker Algorithms:

    Following the discussions within the QEP framework, we have decided to transition algorithms to those already present in the geometry checker. These algorithms encompass topological considerations. Unfortunately, due to timing constraints with the year-end and the feature freeze period, this integration was not yet completed. [4] [5] 
     
    Processing Integration:
    Efforts have been made to incorporate the geometry checker’s algorithms into the processing framework [6] [7] [8]. Regrettably, these efforts have not yet been fully realized. We encountered challenges in aligning processing with the logic of the geometry checker. However, the ongoing PR by Jacky shows promising progress in overcoming these obstacles. We are optimistic about achieving the integration of geometry (and topology) checker algorithms into processing for the upcoming version. This integration aims to streamline ETL processes for verification and correction tasks, thereby enhancing codebase cleanliness.

  5. The remaining project reports will be added as the projects wrap up.

Thank you to everyone who participated and made this round of grants a great success and thank you to all our sustaining members and donors who make this initiative possible!

by underdark at May 04, 2024 09:41 AM

May 03, 2024

A global canopy height map

Meta and the World Resources Institute launched a global map of tree canopy height at a 1-meter resolution, allowing the detection of single trees on a global scale. The good news is that all canopy height data and artificial intelligence models are free and publicly available (Tolan, Couprie, et al. 2024). The maps were created using machine learning models on high-resolution worldwide Maxar satellite imagery. The details of the model used to create the data set are described in Tolan et al. (2024).

The data is available on Google Earth Engine and on AWS as cloud-optimized GeoTIF tiles. In this blog post, I describe how you can download the tiles that intersect with a user-defined area of interest, import the resulting tiles in GRASS GIS, and stitch them together. For some steps, the command line code () and Python code () are provided. For others steps, where results from one command serve as input for the next, only the Python code is provided.

Warning

Before you start, be advised that the global map has some potential issues. If you zoom in, you’ll see conspicuous sudden changes in tree heights. These are more likely to reflect differences in the dates of acquisition of the source imagery than actual differences on the ground. This is easy to see if you zoom in on the Amazon or Congo basins, for example.

Install the software

I assume you already have working knowledge of GRASS GIS. If not, you are advised to start familiarizing yourself with the tool (check out the various resources on the GRASS GIS website).

To download the data from AWS, you first need to download and install the AWS Command Line Interface (AWS CLI) tool. Go to https://aws.amazon.com/cli/, select the installer for your operating system, and follow the instructions.

Set the environment

Start GRASS GIS and create a new project (or locations as it was called prior to version 8.4) with the coordinate reference system (CRS) ESP 3857. Of course, if you already have a project with the required CRS, you can use that.

g.proj epsg=3857 project=PseudoMercator
g.mapset -c mapset=CanopyHeight project=PseudoMercator
import grass.script as gs
gs.run_command("g.proj", epsg=3857, project="PseudoMercator")
gs.run_command("g.mapset", flags="c", mapset="CanopyHeight", 
               project="PseudoMercator")

The global map of tree canopy heights comes in the form of cloud optimized GeoTIFF files. In the next section, you will select and download the tiles for your area of interest to your computer. To make things easier, create a new (temporary) directory to download the files to and make that the working directory.

Note that I am working on Linux. However, apart from the file paths, the commands should be the same on the Windows command line.

mkdir /home/paulo/data/canopyheight
cd /home/paulo/data/canopyheight
import os
os.mkdir('/home/paulo/data/canopyheight')
os.chdir('/home/paulo/data/canopyheight')

List files and folders

The data is available on AWS. Besides the Cloud-optimized GeoTIFF files, there is also a Geojson file with the polygons showing the location of the tiles. You can use the asw-cli tool you installed earlier to get a list of the files and directories that are available. Click the tab to see the resulting list.

aws s3 ls --no-sign-request s3://dataforgood-fb-data/forests/v1/
import subprocess

subprocess.run(
    [
        "aws",
        "s3",
        "ls",
        "--no-sign-request",
        "s3://dataforgood-fb-data/forests/v1/"
    ]
)

To list the content of a sub-folder, run the same code as above, but add the name of the sub-folder for which you want to list the content. For example, you can use the following code to list the content of the sub-folder alsgedi_global_v6_float.

aws s3 ls --no-sign-request s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/
import subprocess

subprocess.run(
    [
        "aws",
        "s3",
        "ls",
        "--no-sign-request",
        "s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/"
    ]
)

If you click the tab above, you can see that the folder contains the tiles.geojson file and the three sub-folders, chm, metadata, and msk.

The chm sub-folder contains the Geotiff tiles of the canopy height map. Their names consist of a numerical code (the QuadKey) followed by the .tiff extension. The tiles.geojson file is a vector layer with polygons showing the location of each tile. For each polygon, the attribute table contains the QuadKey of the corresponding Geotiff file.

The metadata folder contains a corresponding Geojson file [QuadKey].geojson with polygons that represent the observation date of the input imagery for area marked by that polygon. And lastly, the sub-folder msk contains cloud-cover masks accompanying GeoTIFs. The names are composed of the QuadKey.tif.msk.

Select the tiles

You now know in which folder to find what data. The next step is to identify the tiles covering the area of interest. For this tutorial, the area of interest is the Terai Arc Landscape (TAL) region. This is a region that straddles the border of Nepal and India Figure 1.

Figure 1: The lowlands of the Terai Arc Landscape TAL (light green) and the protected areas (dark green), on the foot of the Himalayas.
Figure 1: The lowlands of the Terai Arc Landscape TAL (light green) and the protected areas (dark green), on the foot of the Himalayas.

You can download the geopackage with the boundaries of the region here in the working directory. Replace this with the vector layer representing your area of interest. Next, import it into the mapset you created earlier.

Note that the CRS of the file is latlong ESP 4326. Use the r.import function to reproject the file ‘on-the-fly’ during import.

v.import input=TALregion.gpkg layer=TALregion output=TALregion

Note that the CRS of the file is latlong ESP 4326. Use the r.import function to reproject the file ‘on-the-fly’ during import.

gs.run_command("v.import", input="TALregion.gpkg", 
               layer="TALregion", output="TALregion")

Now, download the tiles.geojson file and import it into GRASS GIS. Note that according to the documentation, the polygon layer has the same CRS (EPSG 3857) as the GeoTIF files. However, in reality, the file has the EPSG 4326 CRS. So also here, use the r.import function.

# Download the file to the working directory
aws s3 cp --no-sign-request s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/tiles.geojson tiles.geojson

# Import the layer in GRASS GIS. 
v.import input=tiles.geojson output=tiles
# Download the file to the working directory
subprocess.run(
    [
        "aws",
        "s3",
        "cp",
        "--no-sign-request",
        "s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/tiles.geojson",
        "tiles.geojson",
    ]
)

# Import the layer in GRASS GIS. 
gs.run_command("v.import", input="tiles.geojson", output="tiles")

The next step is to identify the tiles that overlap with the TAL region, and to get the corresponding QuadKeys. Because you’ll need the list with QuadKeys later on, it is arguably easier to do this in Python.

# Select the tiles that overlap with the TAL region.
gs.run_command("v.select", ainput="tiles", binput="TALregion",
               output="tiles_TAL", operator="overlap")

# Get a list with the QuadKeys from the attribute table
qk = gs.read_command(
    "v.db.select",
    flags="c",
    map="tiles_TAL",
    columns="tile",
    format="plain",
    separator="comma",
).split("\n")

# Remove the empty strings and convert string to integer
qk = [_f for _f in qk if _f]

Download and import the tiles

Now you can download the tiles that overlap with the TAL region. The following code will download the tiles, import them in GRASS GIS 1, and delete the Geotiff.

# Download and import
from os import path

# Download and import the files
baseurl = "s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/chm"
for i, quad in enumerate(qk):
    downloadlink = f"{path.join(baseurl, quad)}.tif"
    layername = f"tile_{quad}"
    subprocess.run(
        [
            "aws",
            "s3",
            "cp",
            "--no-sign-request",
            downloadlink,
            f"{layername}.tif",
        ]
    )
    gs.run_command(
        "r.in.gdal", input=f"{layername}.tif", output=layername, memory=40000
    )
    os.remove(f"{layername}.tif")

Last step is to patch the tiles together into one layer. This requires you to set the region first, using g.region, using the following parameters: With vector = TALregion you set the extent to match that of the TAL region. With raster = qkl the resolution of the output layer will match that of the input tiles. And with align = qkl[0], the resulting layer will be cleanly aligned with the input tiles.

# Create the list of layers
qkl = [f"tile_{x}" for x in qk]

# Set the region to match the extent of the com
gs.run_command("g.region", vector="TALregion", raster=qkl, align=qkl[0])

# Set a mask to exclude the areas outside the boundaries of the 

You are almost done, what remains is using r.patch to patch the tiles together. Set the -s flag to speed up the process. It will disable the time consuming reading and creation of support files. It means that the output map will have no category labels and no explicit color table. Which is fine, as tree heights represent a continuous variable, and you can generate the color table afterwards. Make sure to set the nprocs and memory parameters to fit the specifications of your system.

# Patch the layers together
gs.run_command("r.patch", flags="s", input=qkl, output="CHMtmp", nprocs=10, memory=40000)

# Set color
color_rules = {
    0: "247:252:245",
    3: "229:245:224",
    6: "199:233:192",
    9: "161:217:155",
    12: "116:196:118",
    15: "65:171:93",
    18: "35:139:69",
    21: "0:109:44",
    24: "0:68:27",
    95: "0:0:0",
}
rules_file = gs.tempfile()
with open(rules_file, "w") as f:
    for value, color in color_rules.items():
        f.write(f"{value} {color}\n")
gs.run_command('r.colors', map="CHMtmp", rules=rules_file)

An additional step is to clip the raster to the boundaries of your area of interest. There are different ways to do this, but in the example the r.clip is used.

# Install the r.clip addon
gs.run_command("g.extension", extension="r.clip")

# Set the region to match the extent of the area of interest
gs.run_command("g.region", vector="TAL_region", align="CHM")

# Set a mask to mask out all areas outside the areas of interest
gs.run_command("r.mask", vector="TAL_region")

# Clip the map
gs.run_command("r.clip", input="CHMtmp", output="CHM")

And below the resulting map, an approximately 1 meter resolution canopy height map of the Terai Arc Landscape.

Figure 2: The canopy height map for the Terai Arc Landscape region in Nepal and India.
Figure 2: The canopy height map for the Terai Arc Landscape region in Nepal and India.

Acquisition dates

It is important to realize that the dataset is based on available satellite imagery from 2009 through 2020. The reason is that cloud cover and seasonality impose limitations on the analyzed image dates.

Figure 3: Crown height in part of the TAL region. If you look closely, you can see some sharp boundaries between areas with higher and lower crown heights.
Figure 3: Crown height in part of the TAL region. If you look closely, you can see some sharp boundaries between areas with higher and lower crown heights.
Figure 4: These sharp boundaries are the result of the CHM being based on images taken at different times, rather than reflecting real differences.
Figure 4: These sharp boundaries are the result of the CHM being based on images taken at different times, rather than reflecting real differences.

It is therefore good to know that you can download for each tile a corresponding Geojson file with the observation date of the input imagery.

Figure 5: For each tile [QuadKey].tiff, the metadata folder contains a corresponding Geojson file [QuadKey].geojson with polygons that represent the observation date of the input imagery for area marked by that polygon. In this figure, the observation dates for the input imagery used to create tile 313132000 are shown.
Figure 5: For each tile [QuadKey].tiff, the metadata folder contains a corresponding Geojson file [QuadKey].geojson with polygons that represent the observation date of the input imagery for area marked by that polygon. In this figure, the observation dates for the input imagery used to create tile 313132000 are shown.

The example code below will download the Geojson files and patch them together. The assumption is that you already imported the file with tiles (tiles_TAL). Now, first get the list of quadkeys

qk = gs.read_command(
    "v.db.select",
    flags="c",
    map="tiles_TAL",
    columns="tile",
    format="plain",
    separator="comma",
).split("\n")
qk = [_f for _f in qk if _f]

With the list of quadkeys, you can download and import the geojson files that overlap with your area of interest.

baseurl = "s3://dataforgood-fb-data/forests/v1/alsgedi_global_v6_float/metadata"
for i, quad in enumerate(qk):
    downloadlink = f"{path.join(baseurl, quad)}.geojson"
    layername = f"tile_{quad}"
    subprocess.run(
        [
            "aws",
            "s3",
            "cp",
            "--no-sign-request",
            downloadlink,
            f"{layername}.geojson",
        ]
    )
    gs.run_command("v.import", input=f"{layername}.geojson", output=layername)
    os.remove(f"{layername}.geojson")

qkl = [f"tile_{x}" for x in qk]

After downloading and importing the files, create a new layer by combining the imported vector tiles. When done, remove the intermediate layers.

gs.run_command("v.patch", flags="e", input=qkl, output="CHM_tmp1")
gs.run_command("g.remove", type="vector", pattern="tile_*", flags="f")

You can optionally merge all polygons with the same date, using the v.dissolve function. A caveat is that the column used to dissolve features must be of type integer or string. So, the first step is to create a new column in which you copy the dates from the acq_date column to a new column acq_strdate as strings.

# Convert column with dates to string, required for the r.dissolve function
gs.run_command("v.db.addcolumn", map="CHM_tmp1", columns="acq_strdate varchar(24)")
gs.run_command(
    "db.execute",
    sql="UPDATE CHM_tmp1 SET acq_strdate = CAST(acq_date AS VARCHAR(24))",
)
gs.run_command("v.db.dropcolumn", map="CHM_tmp1", columns="acq_date")

Now, the polygons with the same dates can be dissolved, using the newly created column as key column. Optionally, you can create a new column that holds the dates as dates, rather than strings.

# Dissolve areas with same acquire date
gs.run_command("v.dissolve", input="CHM_tmp1", column="acq_strdate", 
               output="CHM_tmp2")

# Create a new column holding the dates as date
gs.run_command("v.db.addcolumn", map="CHM_tmp2", 
               columns="acq_date date")
gs.run_command("db.execute", 
               sql="UPDATE CHM_tmp2 SET acq_date = a_acq_strdate")
gs.run_command("v.db.dropcolumn", map="CHM_tmp2", 
                columns="a_acq_strdate")

One last optional step is to clip the polygon layer to the boundaries of the area of interest. In this example, this is to clip the CHM_tmp3 vector layer to the boundaries of the TAL region.

gs.run_command("v.clip", input="CHM_tmp2", clip="TAL_PAS", 
               output="CHM_acq_dates")
gs.run_command("g.remove", type="vector", pattern="CHM_tmp*", flags="f")

It is advisable to use the CHM map in conjunction with the vector map of the acquisition dates of source imagery. And for now, I would only compare tree heights if based on the same source image.

Figure 6: The canopy height map for the Terai Arc Landscape region in Nepal and India, overlaid with the vector map of the acquisition dates of source imagery.
Figure 6: The canopy height map for the Terai Arc Landscape region in Nepal and India, overlaid with the vector map of the acquisition dates of source imagery.

References

Tolan, Jamie, Camille Couprie, John Brandt, Justine Spore, Tobias Tiecke, Tracy Johns, and Patrick Nease. 2024. “Using Artificial Intelligence to Map the Earth’s Forests.” Meta Sustainability. https://sustainability.fb.com/blog/2024/04/22/using-artificial-intelligence-to-map-the-earths-forests/.
Tolan, Jamie, Hung-I Yang, Benjamin Nosarzewski, Guillaume Couairon, Huy V. Vo, John Brandt, Justine Spore, et al. 2024. “Very High Resolution Canopy Height Maps from RGB Imagery Using Self-Supervised Vision Transformer and Convolutional Decoder Trained on Aerial Lidar.” Remote Sensing of Environment 300 (January): 113888. https://doi.org/10.1016/j.rse.2023.113888.

Footnotes

  1. An alternative way to deal with the tiles is using gdalbuildvrt to create a virtual mosaic, and r.external to import the virtual layer in GRASS GIS. See the GRASS GIS Wiki page on large data handling for the details.↩︎

by Paulo van Breugel at May 03, 2024 10:00 PM

Eve Kahn shared these amazing pictures from the Awe of the Arctic exhibition currently showing at the New York Public Library (until July 13, 2024)

“Nutjuitok (Polar Star) After Matthew Henson 1866 is a series of photographic works that Terry Adkins dedicated to the African American explorer Matthew Henson. Henson accompanied Robert Peary on seven voyages to the Arctic, including the American explorer’s last and most famous expedition in 1908-09, on which Peary claimed to have reached the geographic North Pole. Henson was only recognized much later in his life for his many contributions to these expeditions. In the photograph Qikiqtaaluk (the Inuktitut name for Baffin Island), a historical map of the Arctic is shown on Adkins’s skin, which has been whitened for the projection. He wears polar bear fur pants and brown fur mittens- two crucial items for surviving in the Arctic and in his hands, he holds large gemstones, a reference to the natural resources explorers sought to extract from the landscape. By embodying the spirit of Henson and the region in this self-portrait, Adkins directly addresses historical relationships to race, capitalism, and imperialism.”Nutjuitok (Polar Star) After Matthew Henson 1866 is a series of photographic works that Terry Adkins dedicated to the African American explorer Matthew Henson. Henson accompanied Robert Peary on seven voyages to the Arctic, including the American explorer’s last and most famous expedition in 1908-09, on which Peary claimed to have reached the geographic North Pole. Henson was only recognized much later in his life for his many contributions to these expeditions. In the photograph Qikiqtaaluk (the Inuktitut name for Baffin Island), a historical map of the Arctic is shown on Adkins’s skin, which has been whitened for the projection. He wears polar bear fur pants and brown fur mittens- two crucial items for surviving in the Arctic and in his hands, he holds large gemstones, a reference to the natural resources explorers sought to extract from the landscape. By embodying the spirit of Henson and the region in this self-portrait, Adkins directly addresses historical relationships to race, capitalism, and imperialism.”

MapsintheWild Nutjuitok

by Steven at May 03, 2024 10:00 AM

May 02, 2024

QGIS Cloud Support is repeatedly asked: ‘Why is my QGIS Cloud map so slow and what can I do to make it faster?’ Basically, the bottlenecks can be traced back to poor data management and/or a slow internet connection. In this short blog article, I would like to give you a few tips on what you can do to speed up a slow loading map. Optimise the data Make sure that your database is optimised and that you only use the necessary data in your web map.

May 02, 2024 12:00 AM