Suporte Postgresql para rasters

From Geomaster, Lda
Revision as of 06:52, 20 February 2018 by Jgrocha (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Suporte Postgis para Raster, workshop realizada no SASIG 2017, no Porto.

http://suite.opengeo.org/docs/latest/dataadmin/pgGettingStarted/raster2pgsql.html

https://duncanjg.wordpress.com/2012/11/20/the-basics-of-postgis-raster/

Ortos 2012

jgr@dusseldorf:/media/bonn/cme/ortos2012/Finais$ gdalinfo 16004200.tif
Driver: GTiff/GeoTIFF
Files: 16004200.tif
       16004200.tfw
Size is 17000, 17000
Coordinate System is `'
Origin = (-40050.000000000000000,125049.899999999994179)
Pixel Size = (0.300000000000000,-0.300000000000000)
Metadata:
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_XRESOLUTION=2116
  TIFFTAG_YRESOLUTION=2116
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -40050.000,  125049.900) 
Lower Left  (  -40050.000,  119949.900) 
Upper Right (  -34950.000,  125049.900) 
Lower Right (  -34950.000,  119949.900) 
Center      (  -37500.000,  122499.900) 
Band 1 Block=128x128 Type=Byte, ColorInterp=Red
  Overviews: 8500x8500, 4250x4250, 2125x2125, 1062x1062, 531x531, 265x265, 132x132, 66x66, 33x33, 16x16
Band 2 Block=128x128 Type=Byte, ColorInterp=Green
  Overviews: 8500x8500, 4250x4250, 2125x2125, 1062x1062, 531x531, 265x265, 132x132, 66x66, 33x33, 16x16
Band 3 Block=128x128 Type=Byte, ColorInterp=Blue
  Overviews: 8500x8500, 4250x4250, 2125x2125, 1062x1062, 531x531, 265x265, 132x132, 66x66, 33x33, 16x16
Band 4 Block=128x128 Type=Byte, ColorInterp=Undefined
  Overviews: 8500x8500, 4250x4250, 2125x2125, 1062x1062, 531x531, 265x265, 132x132, 66x66, 33x33, 16x16

Preparação

Atribuir uma projeção

#!/bin/bash
for f in *.tif
do
	echo "Processing $f"
  	# assign projection to $f
	gdal_edit.py -a_srs "EPSG:3763" $f
done

Notas do Prof. José Alberto Gonçalves

Notas da MapBox, sobre processamento de imagens de satélite.

Exemplo "Ortos do Porto"

gdalbuildvrt estarreja.vrt *.tif

gdal_translate -b 1 -scale 31 255 0 255 estarreja.vrt R.tif
gdal_translate -b 2 -scale 36 255 0 255 estarreja.vrt G.tif
gdal_translate -b 3 -scale 40 255 0 255 estarreja.vrt B.tif

gdalbuildvrt -separate estarreja_all.vrt R.tif G.tif B.tif

# gdalwarp -t_srs "EPSG:3763" -dstalpha -co TILED=YES -co COMPRESS=JPEG estarreja_all.vrt estarreja.tif
# gdalwarp -t_srs "EPSG:3763" -dstalpha -co TILED=YES -co COMPRESS=JPEG -multi estarreja_all.vrt estarreja.tif
gdal_translate estarreja_all.vrt estarreja.tif -b 1 -b 2 -b 3 -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co TILED=YES --config GDAL_TIFF_INTERNAL_MASK YES

gdaladdo -r average estarreja.tif 4 8 16 32 64 128 256

Alternativa (a funcionar com o QGIS e com o Geoserver)

gdalbuildvrt estarreja.vrt *.tif
gdal_translate estarreja.vrt translate/estarreja.tif -b 1 -b 2 -b 3 -mask 4 -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co TILED=YES --config GDAL_TIFF_INTERNAL_MASK YES
gdaladdo -r average translate/estarreja.tif 4 8 16 32 64 128 256

Recriar os tiffs (o geotools estoirava nalgumas imagens)

http://gis.stackexchange.com/questions/89444/file-size-inflation-normal-with-gdalwarp

#!/bin/bash
mkdir warp
for f in *.tif
do
	echo "Processing $f"
  	# gdal_translate -a_srs "EPSG:3763" $f warp/$f
	gdalwarp -t_srs "EPSG:3763" $f warp/$f
done

Raster support PostgreSQL

Nota: Para adicionar ao QGIS, usar o DBManager, em vez do habitual Adicionar camada PostGIS.

Variáveis de ambiente. http://postgis.net/docs/postgis_installation.html#install_short_version

sudo vi /etc/postgresql/9.5/main/environment
POSTGIS_ENABLE_OUTDB_RASTERS=1
POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL

Ortos registados na base de dados (mas guardados fora da base de dados):

raster2pgsql -s 3763 -d -I -C -M -R -l 4 /media/bonn/cme/ortos2012/Finais/15504150.tif -F -t 1000x1000 ortos2012x1000x1000 | psql -h localhost -U geobox -d ortos

-a append -C enforce constraints

raster2pgsql -s 3763 -d -I -M -R -c -l 4 /media/bonn/cme/ortos2012/Finais/16004200.tif -F -t 1000x1000 ortos2012x1000x1000 | psql -h localhost -U geobox -d ortos
raster2pgsql -s 3763 -d -I -M -R -a -l 4 /media/bonn/cme/ortos2012/Finais/15504150.tif -F -t 1000x1000 ortos2012x1000x1000 | psql -h localhost -U geobox -d ortos
raster2pgsql -s 3763 -d -I -C -R -M /media/bonn/cme/ortos2012/Finais/translate/estarreja.tif -F -t 100x100 public.ortos2012 | psql -h localhost -U geobox -d ortos

Dá... Mas fica muito lento!

raster2pgsql -s 3763 -d -I -C -R -M -l 2,4,8,16 /media/bonn/cme/ortos2012/Finais/translate/estarreja.tif -F -t 100x100 public.ortos2012 | psql -h localhost -U geobox -d ortos

Muito mais rápido, mas ainda não está perfeito.

File:Piramides com -l 2 4 8 16.png

raster2pgsql -s 3763 -d -I -C -R -M -l 2,4,8,16,32,64,128,256 /media/bonn/cme/ortos2012/Finais/translate/estarreja.tif -F -t 100x100 public.ortos2012 | psql -h localhost -U geobox -d ortos

Muito mais rápido! Quase perfeito.

File:Piramides com -l 2,4,8,16,32,64,128,256.png

raster2pgsql -s 3763 -d -I -C -R -M -l 2,4,8,16,32,64,128,256,512 /media/bonn/cme/ortos2012/Finais/translate/estarreja.tif -F -t 100x100 public.ortos2012 | psql -h localhost -U sig ortos
No servidor 192.168.100.45
root@cme-sig01:~# adduser geoserver
20170223
root@cme-sig01:~# mkdir ~geoserver/data
root@cme-sig01:~# chmod a+rwx ~geoserver/
root@cme-sig01:~# chmod a+rwx ~geoserver/data
root@cme-sig01:~# mkdir ~geoserver/data/ortos
root@cme-sig01:~# chmod a+rwx ~geoserver/data/ortos
scp estarreja.tif admsig@192.168.100.45:/home/geoserver/data/ortos
apt install postgis
sudo vi /etc/postgresql/9.6/main/environment
POSTGIS_ENABLE_OUTDB_RASTERS=1
POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL
sudo service postgresql restart
raster2pgsql -s 3763 -d -I -C -R -M -l 2,4,8,16,32,64,128,256,512 /home/geoserver/data/ortos/estarreja.tif -F -t 100x100 public.ortos2012 | psql -h localhost -U sig ortos

Virtual layer

Funciona muito bem (no QGIS).

cd /media/bonn/cme/ortos2012/Finais
gdalbuildvrt all_tif.vrt *.tif -a_srs "EPSG:3763"

Gerar uma shapefile com os rectangulos onde encaixam os ortos:

cd /media/bonn/cme/ortos2012/Finais
gdaltindex -t_srs "EPSG:3763" mosaico.shp *.tif

Geoserver: ImageMosaic

Ao criar um store ImageMosaic no Geoserver, é gerada uma shapefile Finais.* (de índice, como a criada pelo gdaltindex) e um documento de configuração Finais.properties.