[gdal-dev] mrsid : coordinate system problem

Didrik Pinte lists at dipole-consulting.com
Wed Sep 17 09:33:10 EDT 2008


Hi,

I have succesfully build the DSDK 7.0 mrsid extension against gdal-svn
and then mapserver 5.2.0 with gdal.

The problem is that the coordinate system of the mrsid file does not
seem to be recognized.

The output of gdalinfo on the file is :
=========================================================================================
[root at ssc data]# gdalinfo ortho/5227y04.sid 
Driver: MrSID/Multi-resolution Seamless Image Database (MrSID)
Files: ortho/5227y04.sid
Size is 10000, 10000
Coordinate System is:
PROJCS["Belge_Lambert_1972",
    GEOGCS["GCS_Belge_1972",
        DATUM["Reseau_National_Belge_1972",
            SPHEROID["International_1924",6378388.0,297.0]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["False_Easting",150000.01256],
    PARAMETER["False_Northing",5400088.4378],
    PARAMETER["Central_Meridian",4.367486666666666],
    PARAMETER["Standard_Parallel_1",49.8333339],
    PARAMETER["Standard_Parallel_2",51.16666733333333],
    PARAMETER["Latitude_Of_Origin",90.0],
    UNIT["Meter",1.0]]
Origin = (-0.001049999962561,21.001049251179211)
Pixel Size = (0.002099999925122,-0.002099999925122)
Metadata:
  IMAGE__INPUT_NAME=N:\DLT\2\P10000\5227Y04.tif
  IMAGE__INPUT_FILE_SIZE=300080326.000000
  IMAGE__COMPRESSION_VERSION=1,4,1
  IMAGE__TARGET_COMPRESSION_RATIO=29.999962
  IMAGE__COMPRESSION_NLEV=5
  IMAGE__COMPRESSION_WEIGHT=4.000000
  IMAGE__COMPRESSION_GAMMA=2.000000
  IMAGE__COMPRESSION_BLOCK_SIZE=512
  IMAGE__CREATION_DATE=Wed Jul 14 07:53:51 2004

Corner Coordinates:
Upper Left  (  -0.0010500,  21.0010493) (  2d18'19.78"E, 49d17'38.76"N)
Lower Left  (  -0.0010500,   0.0010500) (  2d18'19.81"E, 49d17'38.08"N)
Upper Right (  20.9989493,  21.0010493) (  2d18'20.82"E, 49d17'38.78"N)
Lower Right (  20.9989493,   0.0010500) (  2d18'20.85"E, 49d17'38.10"N)
Center      (  10.4989496,  10.5010496) (  2d18'20.31"E, 49d17'38.43"N)
Band 1 Block=1024x128 Type=Byte, ColorInterp=Red
  Minimum=40.000, Maximum=247.000, Mean=91.441, StdDev=25.039
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
  Metadata:
    LAYER_TYPE=athematic
Band 2 Block=1024x128 Type=Byte, ColorInterp=Green
  Minimum=48.000, Maximum=249.000, Mean=94.453, StdDev=23.292
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
  Metadata:
    LAYER_TYPE=athematic
Band 3 Block=1024x128 Type=Byte, ColorInterp=Blue
  Minimum=54.000, Maximum=249.000, Mean=94.428, StdDev=20.477
  Overviews: 5000x5000, 2500x2500, 1250x1250, 625x625, 313x313
  Metadata:
    LAYER_TYPE=athematic
=========================================================================================

Thus, the definition of the raster in mapserver is : 

=========================================================================================
MAP
    IMAGECOLOR 255 255 255 #0 0 0
    SIZE 380 380
    STATUS ON
    SHAPEPATH "/var/data/rasters/"
    RESOLUTION 81.28
    UNITS meters
    FONTSET "/var/data/fontset.txt"

    OUTPUTFORMAT
        NAME png24
        DRIVER "GD/PNG"
        MIMETYPE "image/png"
        IMAGEMODE RGBA
        EXTENSION "png"
        TRANSPARENT OFF
    END

    EXTENT 141193 161468 157995 178169 # extent of URBMAP_PB

    PROJECTION
        "+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90
+lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl
+units=m +no_defs"
    END

  LAYER
        NAME "ortho"
        DATA "ortho/5227y04.sid"
        TYPE RASTER
        PROCESSING "BANDS=1,2,3"
        DEBUG 5
        PROJECTION
                "+proj=lcc +lat_2=51.16666723333333 +lat_1=49.8333339
+lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438
+ellps=intl +units=m +no_defs"
        END
        STATUS DEFAULT
  END
=========================================================================================

When I ask MapServer to render the raster layer on it's extent :

=========================================================================================
[root at aliwen data]# shp2img  -m mrsid.map  -map_debug 9999 -all_debug
9999 -e 153614 164568 154431 165383 o > img.png
[Wed Sep 17 15:25:57 2008].23287 msLoadMap(): 0.003s
[Wed Sep 17 15:25:57 2008].24736 msDrawRasterLayerLow(ortho): entering.
[Wed Sep 17 15:25:57 2008].35788 msResampleGDALToMap(): no overlap ...
no result.
[Wed Sep 17 15:25:57 2008].37246 msDrawMap(): Layer 0 (ortho), 0.013s
[Wed Sep 17 15:25:57 2008].37420 msDrawMap(): Drawing Label Cache,
0.000s
[Wed Sep 17 15:25:57 2008].37532 msDrawMap() total time: 0.014s
[Wed Sep 17 15:25:57 2008].61124 msSaveImage() total time: 0.023s
[Wed Sep 17 15:25:57 2008].61480 msFreeMap(): freeing map at 0x9591c40.
[Wed Sep 17 15:25:57 2008].60993 shp2img total time: 0.042s
=========================================================================================

It says there is no overlap ...

When trying with the following command, I got the output ... :

=========================================================================================
[root at aliwen data]# shp2img  -m mrsid.map  -map_debug 9999 -all_debug
9999 -e 153614 164568 154431 165383 -layer otrho > img.png
[Wed Sep 17 15:25:57 2008].23287 msLoadMap(): 0.003s
[Wed Sep 17 15:25:57 2008].24736 msDrawRasterLayerLow(ortho): entering.
[Wed Sep 17 15:25:57 2008].35788 msResampleGDALToMap(): no overlap ...
no result.
[Wed Sep 17 15:25:57 2008].37246 msDrawMap(): Layer 0 (ortho), 0.013s
[Wed Sep 17 15:25:57 2008].37420 msDrawMap(): Drawing Label Cache,
0.000s
[Wed Sep 17 15:25:57 2008].37532 msDrawMap() total time: 0.014s
[Wed Sep 17 15:25:57 2008].61124 msSaveImage() total time: 0.023s
[Wed Sep 17 15:25:57 2008].61480 msFreeMap(): freeing map at 0x9591c40.
[Wed Sep 17 15:25:57 2008].60993 shp2img total time: 0.042s
=========================================================================================

I have quickly tried to load the data in qgis and the raster appears at
the right place ... 

Do someone have an idea on the source of the problem ? Is this related
to the "Corner Coordinates" that does not seem to be in the projected
coordinate system ?

Thanks in advance for any help.

Didrik

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 197 bytes
Desc: This is a digitally signed message part
Url : http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080917/e8074b36/attachment.bin


More information about the gdal-dev mailing list