[gdal-dev] Meteosat 2nd Gen projection (MSG GEO)

Jose Luis Gomez Dans josegomez at gmx.net
Mon Mar 3 13:42:21 EST 2008


Hi,

Forgot to send this message sorry!

> Jose Luis Gomez Dans wrote:
> > You are absolutely right (as usual!). The file contains two bands, one
> with
> > Latitude and one with Longitude values for each pixel's location. In a
> > sense, it is a set of GCPs stuck in a raster. Since my data are
> referenced
> > with respect to these pixels, I can then work out the Long/Lat values
> from
> > here. i guess that thinking about it, what I would like to do is to
> shove my
> > data into a raster like this (unreference) ENVI raster, and be able to
> > project it into something else. So I could still extract the data from
> my
> > reference ENVI raster into GCPs and use that work out the GeoTransform,
> and
> > assign a projection to my rasters.
> > 
> > Since I have a GCP for each pixel, what would be the best way of dealing
> > with this? Extracting gcp lines is not a problem, how can I tell
> > gdal_translate to use a 2000x2000 (ish) list of GCPs?
> 
> Jose,
> 
> There are a few approaches.
> 
> One is to extract a subset of GCPs (a 10x10 grib of GCPs would likely be
> sufficient).

This looked promising, so I tried it. A simple python script reads my raster file, and sets the GCPs in a VRT file. This file looks like this:

<VRTDataset rasterXSize="2319" rasterYSize="2477">
  <SRS>PROJCS[&quot;unnamed&quot;,
    GEOGCS[&quot;WGS 84&quot;,
        DATUM[&quot;unknown&quot;,
            SPHEROID[&quot;WGS84&quot;,6378137,298.257223563]],
        PRIMEM[&quot;Greenwich&quot;,0],
        UNIT[&quot;degree&quot;,0.0174532925199433]],
    PROJECTION[&quot;Geostationary_Satellite&quot;],
    PARAMETER[&quot;central_meridian&quot;,0],
    PARAMETER[&quot;satellite_height&quot;,35785831],
    PARAMETER[&quot;false_easting&quot;,0],
    PARAMETER[&quot;false_northing&quot;,0]]</SRS>
  <Metadata/>
  <GCPList Projection="PROJCS[&quot;unnamed&quot;,
    GEOGCS[&quot;WGS 84&quot;,
        DATUM[&quot;unknown&quot;,
            SPHEROID[&quot;WGS84&quot;,6378137,298.257223563]],
        PRIMEM[&quot;Greenwich&quot;,0],
        UNIT[&quot;degree&quot;,0.0174532925199433]],
    PROJECTION[&quot;Geostationary_Satellite&quot;],
    PARAMETER[&quot;central_meridian&quot;,0],
    PARAMETER[&quot;satellite_height&quot;,35785831],
    PARAMETER[&quot;false_easting&quot;,0],
    PARAMETER[&quot;false_northing&quot;,0]]">
    <GCP Id="1" Pixel="1.0000" Line="1.0000" X="-2.708186721802E+01" Y="3.807548141479E+01"/>
    <GCP Id="2" Pixel="1001.0000" Line="1.0000" X="9.776592254639E+00" Y="3.738645935059E+01"/>
    <GCP Id="3" Pixel="2001.0000" Line="1.0000" X="6.308378601074E+01" Y="4.104884719849E+01"/>
    <GCP Id="4" Pixel="1.0000" Line="1001.0000" X="-2.030329895020E+01" Y="6.382459163666E+00"/>
    <GCP Id="5" Pixel="1001.0000" Line="1001.0000" X="7.539474964142E+00" Y="6.312354087830E+00"/>
    <GCP Id="6" Pixel="2001.0000" Line="1001.0000" X="3.916501617432E+01" Y="6.590889930725E+00"/>
    <GCP Id="7" Pixel="1.0000" Line="2001.0000" X="-2.211191177368E+01" Y="-2.196564483643E+01"/>
    <GCP Id="8" Pixel="1001.0000" Line="2001.0000" X="8.158615112305E+00" Y="-2.168662071228E+01"/>
    <GCP Id="9" Pixel="2001.0000" Line="2001.0000" X="4.384607315063E+01" Y="-2.284282493591E+01"/>
  </GCPList>
  <VRTRasterBand dataType="Int16" band="1">
    <Metadata/>
    <Description>Resize (Band 1:MSG1-SEVI-MSG15-0100-NA-20040206091237.690000000Z-3382.natHRVDNs.img)</Description>
    <SimpleSource>
      <SourceFilename relativeToVRT="0">Input_file.dat</SourceFilename>
      <SourceBand>1</SourceBand>
      <SrcRect xOff="0" yOff="0" xSize="2319" ySize="2477"/>
      <DstRect xOff="0" yOff="0" xSize="2319" ySize="2477"/>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

So, it seems to make sense. If I now translate the VRT into say WGS84, I get a very strange projection issue: 
gdalwarp -rn -of GTiff -t_srs "EPSG:4326" <input.vrt> <output.tif>
gdalinfo <output.tif>
Driver: GTiff/GeoTIFF
Size is 3206, 2429
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235629972,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-0.000237909788845,0.000372871476952)
Pixel Size = (0.000000284910289,-0.000000284910289)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  -0.0002379,   0.0003729) (  0d 0'0.86"W,  0d 0'1.34"N)
Lower Left  (  -0.0002379,  -0.0003192) (  0d 0'0.86"W,  0d 0'1.15"S)
Upper Right (   0.0006755,   0.0003729) (  0d 0'2.43"E,  0d 0'1.34"N)
Lower Right (   0.0006755,  -0.0003192) (  0d 0'2.43"E,  0d 0'1.15"S)
Center      (   0.0002188,   0.0000268) (  0d 0'0.79"E,  0d 0'0.10"N)
Band 1 Block=3206x1 Type=Int16, ColorInterp=Gray

This is clearly wrong, so there is something I'm not doing right. Does anyone have any clues?

Cheers,
Jose
-- 
Der GMX SmartSurfer hilft bis zu 70% Ihrer Onlinekosten zu sparen! 
Ideal für Modem und ISDN: http://www.gmx.net/de/go/smartsurfer


More information about the gdal-dev mailing list