[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["unnamed",
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Geostationary_Satellite"],
PARAMETER["central_meridian",0],
PARAMETER["satellite_height",35785831],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0]]</SRS>
<Metadata/>
<GCPList Projection="PROJCS["unnamed",
GEOGCS["WGS 84",
DATUM["unknown",
SPHEROID["WGS84",6378137,298.257223563]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433]],
PROJECTION["Geostationary_Satellite"],
PARAMETER["central_meridian",0],
PARAMETER["satellite_height",35785831],
PARAMETER["false_easting",0],
PARAMETER["false_northing",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