Mostrando postagens com marcador GDAL. Mostrar todas as postagens
Mostrando postagens com marcador GDAL. Mostrar todas as postagens

sábado, março 01, 2014

Mapas de relevo com a biblioteca GDAL

Elaborei um script utilizando as bibliotecas gdal e imagemagick para criar mapas de relevo com o efeito Hillshade. Para rodá-lo, basta colocá-lo em uma pasta (como mostra a Figura 1) com os arquivos SRTM baixados da site da Embrapa, junto com o shapefile contendo o limite da sua área de interesse e esta palheta de cores. Em seguida, é só digitar ./mosaic.sh

Figura 1 - Pasta contendo os arquivos necessários

Em função das cores definidas para cada cota no arquivo my_classes.ctp, o script cria um novo raster através do comando gdaldem (Figura 2).


Figura 2 - Imagem SRTM colorida

Em seguida é gerado o arquivo Hillshade (Figura 3).

Figura 3 - Relevo sombreado (Hillshade)

E os arquivos anteriores são unidos, como mostra a Figura 4.

Figura 4 - União do relevo sombreado com a imagem SRTM colorida

Por último, em função do limite do shapefile, o arquivo é "clipado". A Figura 5 mostra o resultado da utilização do script para o estado do Mato Grosso do Sul.

Figura 5 - Resultado da utilização do script para o MS.

Lembrando que para utilizar o script, é necessário lhe dar permissões de execução após baixá-lo, através do comando abaixo:
$ sudo chmod +x mosaic.sh

Referências:

domingo, agosto 04, 2013

Renomeando e Ordenando Campos de um Shapefile via Terminal

Quando trabalhamos com dados vetoriais, frequentemente temos que renomear e ordenar os campos que compõem um arquivo .dbf de um shapefile. O que muita gente não sabe é que isso pode ser feito sem a necessidade de abrir um software de SIG, utilizando a biblioteca OGR através do seguinte comando:
ogr2ogr -a_srs <EPSG:num> <novo_shapefile> <shapefile> -sql "SELECT campo1 AS novo_campo1, campo_2 AS novo_campo2, ..., campoN AS novo_campoN FROM shapefile"
Isto pode ser útil principalmente se o arquivo tiver muitos megabytes (MB), pois este processo pode ser bem demorado, independente do software que você esteja utilizando. Vejamos o exemplo a seguir:

Primeiro utilizando o comando ogrinfo, para obter as informações necessárias sobre o shapefile munipb.shp:
ogrinfo -summary munipb.shp munipb
INFO: Open of `munipb.shp'
      using driver `ESRI Shapefile' successful.

Layer name: munipb
Geometry: Polygon
Feature Count: 224
Extent: (-38.764870, -8.302541) - (-34.793335, -6.026173)
Layer SRS WKT:
GEOGCS["GCS_South_American_1969",
    DATUM["South_American_Datum_1969",
        SPHEROID["GRS_1967_Modified",6378160.0,298.25]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]]
NOMEMUNIC: String (50.0)
ARp: Integer (10.0)
UF: String (2.0)
ID_UF: Real (11.0)
REGIAO: String (20.0)
MESOREGIAO: String (50.0)
MICROREGIA: String (50.0)
geocode: String (12.0)
Area_km2: Real (6.2)
per_km: Real (6.2)
id_: Integer (10.0)
O novo shapefile (municipios.shp) será reestruturado com os seguintes campos fid (campo do tipo autoincremento criado automaticamente), geocodigo (geocode), nome (NOMEMUNIC),  micro (MICROREGIA), meso (MESOREGIAO), e area_km2 (Area_km2):
ogr2ogr -a_srs EPSG:4291 municipios.shp munipb.shp -sql "SELECT fid, geocode AS geocodigo, nomemunic AS nome, microregia AS micro, mesoregiao AS meso, area_km2 AS area_km2 FROM munipb"
O resultado pode ser visto com o comando ogrinfo ou em qualquer software de SIG (Figura 1):
ogrinfo -summary municipios.shp municipios
INFO: Open of `municipios.shp'
      using driver `ESRI Shapefile' successful.

Layer name: municipios
Geometry: Polygon
Feature Count: 224
Extent: (-38.764870, -8.302541) - (-34.793335, -6.026173)
Layer SRS WKT:
GEOGCS["SAD69",
    DATUM["South_American_Datum_1969",
        SPHEROID["GRS_1967",6378160,298.247167427]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]]
fid: Integer (10.0)
geocodigo: String (12.0)
nome: String (50.0)
micro: String (50.0)
meso: String (50.0)
area_km2: Real (6.2)
Figura 1 - Tabela reestruturada aberta no QGIS

terça-feira, janeiro 31, 2012

Reprojeção e Clip de Imagens "On the Fly" com o GRASS e GDAL

Lendo o livro "Open Source GIS: A GRASS GIS Approach" (pág. 49) encontrei um procedimento muito útil para resolver o seguinte problema: recortar uma imagem SRTM e reprojetá-la para os mesmos parâmetros definidos no meu projeto. Dessa forma fica bem mais fácil importar a imagem gerada com o comando r.in.gdal. Segue o passo a passo:

1 - Verificar as informações da imagem a ser reprojetada e recortada:
GRASS 7.0.svn (gramame):~ > gdalinfo Geoprocessamento/raster/srtm/SB-25-Y-C.tif
Driver: GTiff/GeoTIFF
Files: Geoprocessamento/raster/srtm/SB-25-Y-C.tif
Size is 1800, 1200
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-36.000415601243731,-7.000415916030761)
Pixel Size = (0.000833333353512,-0.000833333353512)
Metadata:
  TIFFTAG_SOFTWARE=IMAGINE TIFF Support
Copyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved
@(#)$RCSfile: etif.c $ $Revision: 1.10.1.9 $ $Date: 2003/11/13 19:25:49EST $
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -36.0004156,  -7.0004159) ( 36d 0' 1.50"W,  7d 0' 1.50"S)
Lower Left  ( -36.0004156,  -8.0004159) ( 36d 0' 1.50"W,  8d 0' 1.50"S)
Upper Right ( -34.5004156,  -7.0004159) ( 34d30' 1.50"W,  7d 0' 1.50"S)
Lower Right ( -34.5004156,  -8.0004159) ( 34d30' 1.50"W,  8d 0' 1.50"S)
Center      ( -35.2504156,  -7.5004159) ( 35d15' 1.50"W,  7d30' 1.50"S)
Band 1 Block=1800x2 Type=UInt16, ColorInterp=Gray

2 - Imprimir a região definida como padrão na location:
GRASS 7.0.svn (gramame):~ > g.region -p
projection: 1 (UTM)
zone:       -25
datum:      sam69
ellipsoid:  sam69
north:      9206640
south:      9180540
west:       259380
east:       300780
nsres:      90
ewres:      90
rows:       290
cols:       460
cells:      133400

3 - Usar o comando "eval", para setar como variáveis, a saída do comando anterior (veja mais sobre esse comando em: http://goo.gl/70sAi):
GRASS 7.0.svn (gramame):~ > eval `g.region -g`


4 - Verificar o resultado do comando anterior com o comando "echo", observe que é exibido o retângulo envolvente da região definida como default:
GRASS 7.0.svn (gramame):~ > echo "$w $s $e $n"
259380 9180540 300780 9206640


5 - Verificar a projeção definida na location:
GRASS 7.0.svn (gramame):~ > g.proj -wf
PROJCS["UTM Zone 25, Southern Hemisphere",GEOGCS["sam69",DATUM["South_American_Datum_1969",SPHEROID["South_American_Datum_1969",6378160,298.25]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-33],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]


6 - Usar o comando "gdalwarp" para recortar e reprojetar a imagem a partir da projeção definida no grass e do retângulo envolvente da region:
GRASS 7.0.svn (gramame):~ > gdalwarp -t_srs "`g.proj -wf`" -te $w $s $e $n \
Geoprocessamento/raster/srtm/SB-25-Y-C.tif  \ 
Geoprocessamento/raster/srtm/srtm_clip_utm25.tif
Creating output file that is 450P x 283L.
Processing input file Geoprocessamento/raster/srtm/SB-25-Y-C.tif.
0...10...20...30...40...50...60...70...80...90...100 - done.


7 - Verificar a imagem reprojetada e cortada:
GRASS 7.0.svn (gramame):~ > gdalinfo Geoprocessamento/raster/srtm/srtm_clip_utm25.tif 
Driver: GTiff/GeoTIFF
Files: Geoprocessamento/raster/srtm/srtm_clip_utm25.tif
Size is 450, 283
Coordinate System is:
PROJCS["UTM Zone 25, Southern Hemisphere",
    GEOGCS["sam69",
        DATUM["South_American_Datum_1969",
            SPHEROID["South_American_Datum_1969",6378160,298.25,
                AUTHORITY["EPSG","7036"]],
            AUTHORITY["EPSG","6291"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-33],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (259380.000000000000000,9206640.000000000000000)
Pixel Size = (92.000000000000000,-92.226148409893995)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  259380.000, 9206640.000) ( 35d10'43.72"W,  7d10'19.95"S)
Lower Left  (  259380.000, 9180540.000) ( 35d10'47.83"W,  7d24'29.27"S)
Upper Right (  300780.000, 9206640.000) ( 34d48'14.67"W,  7d10'25.80"S)
Lower Right (  300780.000, 9180540.000) ( 34d48'18.07"W,  7d24'35.32"S)
Center      (  280080.000, 9193590.000) ( 34d59'31.09"W,  7d17'27.72"S)
Band 1 Block=450x9 Type=UInt16, ColorInterp=Gray


8 - Importar a imagem reprojetada e cortada:
GRASS 7.0.svn (gramame):~ > r.in.gdal input=Geoprocessamento/raster/srtm/srtm_clip_utm25.tif output=mde
Projection of input dataset and current location appear to match
 100%
r.in.gdal complete. Raster map created.

Um abraço e até o próximo post :)

sexta-feira, dezembro 17, 2010

Alternativa Open Source para o ArcCatalog e o ArcToolbox

O wxGIS possui as mesmas características dos módulos ArcCatalog e ArcToolBox do ArcGIS, e o que é melhor, é open source :)





Para instalar o wxGIS no Ubuntu é necessário: