0

I'm trying to re-project Meteor-M 2 satellite image with gdalwrap. I have GCP points which I put into .vrt file:

<VRTDataset rasterXSize="1568" rasterYSize="480">
  <GCPList Projection="EPSG:4326">
    <GCP Id="0" Pixel="1568.5" Line="0.5" X="62.20476647354344" Y="42.09684997440528" Z="0.0" />
    ... lots of them
    </GCPList>
    <VRTRasterBand dataType="Byte" band="1">
      <SimpleSource>
        <SourceFilename relativeToVRT="1">output.png</SourceFilename>
        <SourceBand>1</SourceBand>
      </SimpleSource>
    </VRTRasterBand>
</VRTDataset>

Then I executed gdalwrap with the following arguments:

gdalwarp -tps -overwrite -t_srs epsg:3857 -of GTIFF output.vrt output10.tif

Then I cut this geotiff into tiles:

gdal2tiles.py --profile=mercator -z 3-8 output10.tif tiles

So far so good. But when I draw both GCP and tiles onto open street map I'm getting:

meteor-m2

Re-projected image and GCPs differ greatly.

gdalinfo:

Driver: GTiff/GeoTIFF
Files: output10.tif
Size is 1656, 583
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (2926801.391389018390328,6009567.154649123549461)
Pixel Size = (2411.811113223014218,-2411.811113223014218)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( 2926801.391, 6009567.155) ( 26d17'30.86"E, 47d24'42.83"N)
Lower Left  ( 2926801.391, 4603481.276) ( 26d17'30.86"E, 38d10'12.31"N)
Upper Right ( 6920760.595, 6009567.155) ( 62d10'12.90"E, 47d24'42.83"N)
Lower Right ( 6920760.595, 4603481.276) ( 62d10'12.90"E, 38d10'12.31"N)
Center      ( 4923780.993, 5306524.215) ( 44d13'51.88"E, 42d57'51.12"N)
Band 1 Block=1656x4 Type=Byte, ColorInterp=Gray

So the question is:

How to re-project wide-angle camera? Looks like gdal incorrectly process it. I.e. pixel size in the middle of the picture is ~1km. While gdalinfo makes it ~2km. Which might be ok on the edges.

Another question:

Why only "-tps" can handle non-linear distortions on the edges?

0

Issue was on my side. When I was generating GCPs, I took (originalImageHeight / someRate) points. Effectively missing the last line. So the correct number of lines should be (originalImageHeight / someRate) + 1.

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.