[Gfoss] Riproiettare 2 (Ovvero, .bil e grid)

Mon Sep 10 16:10:33 CEST 2007

(cc lista per uso di GDAL tools)

On Mon, Sep 10, 2007 at 03:24:15AM -0700, Luca Moiana wrote:
> Ciao Markus,
>
> Sono riuscito con gdalwarp a riproiettare un gtiff e usarlo per le mie 
> carte.
>
> Ti allego comunque i files zippati, mi ineteressa sapere come riproiettarli 
> nella loro forma originale.

allora:

## file originale:
gdalinfo mydata.bil
Driver: EHdr/ESRI .hdr Labelled
Files: mydata.bil
       mydata.hdr
       mydata.clr
Size is 360, 360
Coordinate System is `'
Origin = (14.999999995000000,9.999999995000000)
Pixel Size = (0.008333330000000,-0.008333330000000)
Corner Coordinates:
Upper Left  (  15.0000000,  10.0000000)
Lower Left  (  15.0000000,   7.0000012)
Upper Right (  17.9999988,  10.0000000)
Lower Right (  17.9999988,   7.0000012)
Center      (  16.4999994,   8.5000006)
Band 1 Block=360x1 Type=UInt32, ColorInterp=Palette
  NoData Value=-500
  Color Table (RGB with 9253 entries)
    0: 0,0,0,255
    1: 0,150,150,255
    2: 0,150,150,255
...

## ... il file originale NON contiene la proiezone.
## Sembra LatLong (ma anche WGS84? o qualcos'altro di ellisoide?).
##
## Assegno LatLong/WGS84 (spero che sia giusto), scrivo un GeoTIFF:
gdal_translate -a_srs epsg:4326 mydata.bil mydata.tif
Input file size is 360, 360
0Warning 1: Unable to export color table to GeoTIFF file.  Color tables
can only be written to 1 band Byte or UInt16 GeoTIFF files.
...10...20...30...40...50...60...70...80...90...100 - done.

## .... si vede che GeoTIFF non va bene per i colori.

## Allora, faccio un ERDAS/IMG file:
gdal_translate -a_srs epsg:4326 -of HFA mydata.bil mydata_latlong.hfa
Input file size is 360, 360
0...11...22..30...41..50...61...72..80...91..100 - done.

## Controllo:
gdalinfo mydata_latlong.hfa
Driver: HFA/Erdas Imagine Images (.img)
Files: mydata_latlong.hfa
Size is 360, 360
Coordinate System is:
GEOGCS["WGS_1984",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016],
        TOWGS84[0,0,0,0,0,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (14.999999995000000,9.999999995000000)
Pixel Size = (0.008333330000000,-0.008333330000000)
Corner Coordinates:
Upper Left  (  15.0000000,  10.0000000) ( 15d 0'0.00"E, 10d 0'0.00"N)
Lower Left  (  15.0000000,   7.0000012) ( 15d 0'0.00"E,  7d 0'0.00"N)
Upper Right (  17.9999988,  10.0000000) ( 18d 0'0.00"E, 10d 0'0.00"N)
Lower Right (  17.9999988,   7.0000012) ( 18d 0'0.00"E,  7d 0'0.00"N)
Center      (  16.4999994,   8.5000006) ( 16d30'0.00"E,  8d30'0.00"N)
Band 1 Block=64x64 Type=UInt32, ColorInterp=Palette
  Description = Layer_1
  Metadata:
    LAYER_TYPE=thematic
  Color Table (RGB with 9253 entries)
    0: 0,0,0,255
    1: 0,150,150,255
    2: 0,150,150,255
    3: 0,150,150,255
...

## ... bello. Cosi' il file e' pronto per 'gdalwarp':
##
## Riproiezione in Utm33 N:
## Il codice EPSG prendo dal file /usr/local/share/proj/epsg:

# WGS 84 / UTM zone 33N
<32633> +proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs  <>

## Attenzione: bisogna (meglio) indicare la resoluzione metrica
## La risoluzione 0.0083333333 (vedi sopra) corresponde a Trento a
## nsres=926.27658097 meters, ewres=644.38333656 meters
##
## Se voi tenere raster celle quadratiche, devi indicare bene la
## risoluzione in uscita con -tr xres yres (p.e. 750m x 750m).
## Ci sono alcuno metodi da sceglere per il resampling:
##   "Available resampling methods:
##      near (default), bilinear, cubic, cubicspline, lanczos.
##   "
## Vedi anche questa pagina dove c'e' un paragone tra i metodi:
##   http://169.237.35.250/~dylan/gdalwarp/image_matrix.html
##
gdalwarp -of HFA -t_srs epsg:32633 -tr 750 750 mydata_latlong.hfa mydata_utm33n_wga84.hfa
Copying color table from mydata_latlong.hfa to new file.
Creating output file that is 360P x 362L.
Processing input file mydata_latlong.hfa.
:0...10...20...30...40...50...60...70...80...90...100 - done.

## controllo:
gdalinfo mydata_utm33n_wga84.hfa
Driver: HFA/Erdas Imagine Images (.img)
Files: mydata_utm33n_wga84.hfa
Size is 441, 444
Coordinate System is:
PROJCS["UTM Zone 33, Northern Hemisphere",
    GEOGCS["WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.2572235630016],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1],
    AUTHORITY["EPSG","32633"]]
Origin = (499999.999447746900842,1106908.852485561044887)
Pixel Size = (750.000000000000000,-750.000000000000000)
Corner Coordinates:
Upper Left  (  499999.999, 1106908.852) ( 15d 0'0.00"E, 10d 0'48.72"N)
Lower Left  (  499999.999,  773908.852) ( 15d 0'0.00"E,  7d 0'5.21"N)
Upper Right (  830749.999, 1106908.852) ( 18d 0'59.74"E,  9d59'59.46"N)
Lower Right (  830749.999,  773908.852) ( 17d59'35.37"E,  6d59'30.94"N)
Center      (  665374.999,  940408.852) ( 16d30'8.75"E,  8d30'16.73"N)
Band 1 Block=64x64 Type=UInt32, ColorInterp=Palette
  Description = Layer_1
  Metadata:
    LAYER_TYPE=thematic
  Color Table (RGB with 9253 entries)
    0: 0,0,0,255
    1: 0,150,150,255
    2: 0,150,150,255
...

## viz:
qgis mydata_utm33n_wga84.hfa

Joy and happiness...
tutto a posto, funziona.

ciao
Markus

-- 
Markus Neteler  <neteler itc it>  http://mpa.itc.it/markus/
FBK-irst -  Centro per la Ricerca Scientifica e Tecnologica
MPBA - Predictive Models for Biol. & Environ. Data Analysis
Via Sommarive, 18        -       38050 Povo (Trento), Italy

------------------
ITC -> dall'1 marzo 2007 Fondazione Bruno Kessler
ITC -> since 1 March 2007 Fondazione Bruno Kessler
------------------