[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