[mapserver-users] Re: Projection issue with vrt file

Marc-Andre Morin Marc-Andre.Morin at dfo-mpo.gc.ca
Fri Apr 27 11:53:42 EDT 2012 

Hi Frank,

Thanks for your reply.  I've made some tests and got some more details about
this issue:

1- If I translate my BSB file (KAP format) to Geotiff format, this geotiff
file is well shown and georeferenced at the right place if I use it
directly:

	DATA "BSB/KAP/131701.kap.warped.tiff"	
	PROJECTION
		'proj=merc'
		'lon_0=-70'
		'k=1'
		'x_0=0'
		'y_0=0'
		'ellps=WGS84'
		'towgs84=0,0,0,0,0,0,0'
		'units=m'
		'no_defs'
        END


Command line to generate the geotiff file:
gdal_translate 131701.kap 131701.kap.ori.tiff
gdalwarp 131701.kap.ori.tiff 131701.kap.warped.tiff


Gdalinfo...

E:\ms4w\apps\marinfo\data\BSB\KAP>gdalinfo 131701.kap.warped.tiff
Driver: GTiff/GeoTIFF
Files: 131701.kap.warped.tiff
Size is 11387, 8093
Coordinate System is:
PROJCS["Global Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",-70],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-107047.991228487620000,5945211.986471482600000)
Pixel Size = (7.342576435115765,-7.342576435115765)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -107047.991, 5945211.986) ( 70d57'41.86"W, 47d12'39.62"N)
Lower Left  ( -107047.991, 5885788.515) ( 70d57'41.86"W, 46d50'45.62"N)
Upper Right (  -23438.073, 5945211.986) ( 70d12'37.97"W, 47d12'39.62"N)
Lower Right (  -23438.073, 5885788.515) ( 70d12'37.97"W, 46d50'45.62"N)
Center      (  -65243.032, 5915500.251) ( 70d35'9.92"W, 47d 1'43.75"N)
Band 1 Block=11387x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
...


2- Now, the VRT file that I try to generate and use comes from the
tilers_tools' python script called map2gdal.py
http://code.google.com/p/tilers-tools/

The big advantage with this VRT file is the fact that this VRT file points
directly to my KAP file, and all metadata related to this file are extracted
and available for future processes...

This VRT file looks like this:

<VRTDataset rasterXSize="11388" rasterYSize="8091">
  <SRS>+proj=merc +lat_ts=47.166667 +datum=NAD83 +nodefs</SRS>
  <GeoTransform>-5.3802490306424974e+006,
5.0001511315648770e+000,-1.8211665160056521e-004,
4.0492551038299361e+006,-2.5595698020313497e-004,-5.0010081232801378e+000</GeoTransform>
  <Metadata>
    <MDI
key="BSB_KNP">SC=50000,GD=NAD83,PR=MERCATOR,PP=47.166667,PI=5.000000,SP=,SK=0.000000,TA=90.000000,UN=METRES,SD=LNT,DX=5.000000,DY=5.000000</MDI>
    <MDI
key="BSB_KNQ">EC=RF,GD=NAR,VC=ZZZ,SC=ZZZ,PC=MC,P1=-70.582029,P2=47.166667,P3=0.0,P4=NOT_APPLICABLE,GC=GE,RM=NONE</MDI>
    <MDI key="CUTLINE">MULTIPOLYGON(((257.76284645559468
7745.7076975815671,2490.025011025306 7746.5301512821934,2492.2163879532559
7181.6425479723066,3276.1412438760321 7210.5951029770404,4429.0270760545409
6817.9091369725729,4890.9737230091405 6736.5650032022368,11250.811878184026
354.19735628215676,11251.785400765162 244.45756921568665,3537.7418346315158
244.33482436588929,3523.5069250460119 2029.8367618071302,1881.7430005409099
2028.8818991008357,1868.6640235103462 2447.8575431625686,258.92730402280057
2434.9666589315798)))</MDI>
    <MDI key="DESCRIPTION">SAULT-AU-COCHON À/TO QUÉBEC</MDI>
  </Metadata>
  <GCPList Projection="+proj=merc +lat_ts=47.166667 +datum=NAD83 +nodefs">
    ...
  </GCPList>
  <VRTRasterBand dataType="Byte" band="1">
    <Metadata />
    <ColorInterp>Palette</ColorInterp>
    <ColorTable>
      <Entry c1="0" c2="0" c3="0" c4="255" />
      <Entry c1="255" c2="255" c3="255" c4="255" />
      <Entry c1="255" c2="214" c3="140" c4="255" />
      <Entry c1="107" c2="165" c3="222" c4="255" />
      <Entry c1="182" c2="213" c3="240" c4="255" />
      <Entry c1="100" c2="160" c3="140" c4="255" />
      <Entry c1="0" c2="86" c3="158" c4="255" />
      <Entry c1="75" c2="75" c3="75" c4="255" />
      <Entry c1="123" c2="8" c3="90" c4="255" />
      <Entry c1="100" c2="100" c3="100" c4="255" />
      <Entry c1="243" c2="137" c3="71" c4="255" />
      <Entry c1="44" c2="99" c3="207" c4="255" />
      <Entry c1="245" c2="255" c3="80" c4="255" />
      <Entry c1="255" c2="0" c3="0" c4="255" />
    </ColorTable>
    <SimpleSource>
      <SourceFilename
relativeToVRT="0">E:\ms4w\apps\marinfo\data\BSB\KAP\131701.kap</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="11388" RasterYSize="8091"
DataType="Byte" BlockXSize="11388" BlockYSize="1" />
      <SrcRect xOff="0" yOff="0" xSize="11388" ySize="8091" />
      <DstRect xOff="0" yOff="0" xSize="11388" ySize="8091" />
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>


If I change the mapfile to now point on my VRT file...

	DATA "BSB/KAP/131701.vrt"	
	PROJECTION
		'proj=merc'
		'lon_0=-70'
		'k=1'
		'x_0=0'
		'y_0=0'
		'ellps=WGS84'
		'towgs84=0,0,0,0,0,0,0'
		'units=m'
		'no_defs'
        END

...  Nothing appears on the map.


3- If I change the GeoTransform element by the one generated with gdal from
the geotiff file


gdalbuildvrt 131701.vrt 131701.kap.warped.tiff


131701.vrt:

<VRTDataset rasterXSize="11387" rasterYSize="8093">
  <SRS>PROJCS[&quot;Global Mercator&quot;,GEOGCS[&quot;WGS
84&quot;,DATUM[&quot;WGS_1984&quot;,SPHEROID[&quot;WGS
84&quot;,6378137,298.257223563,AUTHORITY[&quot;EPSG&quot;,&quot;7030&quot;]],AUTHORITY[&quot;EPSG&quot;,&quot;6326&quot;]],PRIMEM[&quot;Greenwich&quot;,0],UNIT[&quot;degree&quot;,0.0174532925199433],AUTHORITY[&quot;EPSG&quot;,&quot;4326&quot;]],PROJECTION[&quot;Mercator_1SP&quot;],PARAMETER[&quot;central_meridian&quot;,-70],PARAMETER[&quot;scale_factor&quot;,1],PARAMETER[&quot;false_easting&quot;,0],PARAMETER[&quot;false_northing&quot;,0],UNIT[&quot;metre&quot;,1,AUTHORITY[&quot;EPSG&quot;,&quot;9001&quot;]]]</SRS>
  <GeoTransform>-1.0704799122848762e+005, 7.3425764351157650e+000,
0.0000000000000000e+000, 5.9452119864714826e+006,
0.0000000000000000e+000,-7.3425764351157650e+000</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1">
    <ColorInterp>Palette</ColorInterp>
    ...
    </ColorTable>
    <SimpleSource>
      <SourceFilename
relativeToVRT="1">131701.kap.warped.tiff</SourceFilename>
      <SourceBand>1</SourceBand>
      <SourceProperties RasterXSize="11387" RasterYSize="8093"
DataType="Byte" BlockXSize="11387" BlockYSize="1"/>
      <SrcRect xOff="0" yOff="0" xSize="11387" ySize="8093"/>
      <DstRect xOff="0" yOff="0" xSize="11387" ySize="8093"/>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>


So the new VRT file will look like this:

<VRTDataset rasterXSize="11388" rasterYSize="8091">
  <SRS>+proj=merc +lat_ts=47.166667 +datum=NAD83 +nodefs</SRS>
  <GeoTransform>-1.0704799122848762e+005, 7.3425764351157650e+000,
0.0000000000000000e+000, 5.9452119864714826e+006,
0.0000000000000000e+000,-7.3425764351157650e+000</GeoTransform>
  <Metadata>
    <MDI
key="BSB_KNP">SC=50000,GD=NAD83,PR=MERCATOR,PP=47.166667,PI=5.000000,SP=,SK=0.000000,TA=90.000000,UN=METRES,SD=LNT,DX=5.000000,DY=5.000000</MDI>
    <MDI
key="BSB_KNQ">EC=RF,GD=NAR,VC=ZZZ,SC=ZZZ,PC=MC,P1=-70.582029,P2=47.166667,P3=0.0,P4=NOT_APPLICABLE,GC=GE,RM=NONE</MDI>
    <MDI key="CUTLINE">MULTIPOLYGON(((257.76284645559468
7745.7076975815671,2490.025011025306 7746.5301512821934,2492.2163879532559
7181.6425479723066,3276.1412438760321 7210.5951029770404,4429.0270760545409
6817.9091369725729,4890.9737230091405 6736.5650032022368,11250.811878184026
354.19735628215676,11251.785400765162 244.45756921568665,3537.7418346315158
244.33482436588929,3523.5069250460119 2029.8367618071302,1881.7430005409099
2028.8818991008357,1868.6640235103462 2447.8575431625686,258.92730402280057
2434.9666589315798)))</MDI>
    <MDI key="DESCRIPTION">SAULT-AU-COCHON À/TO QUÉBEC</MDI>
  </Metadata>
...
<SourceFilename
relativeToVRT="0">D:\ms4w\apps\marinfo\data\BSB\KAP\131701.kap</SourceFilename>
...


And now, the BSB file is shown properly...

Why?  Is it because the BSB geotransform element has originally a rotation
in the matrix, that is eliminated when the warped geotiff file is generated?

Effectively, if I try to generate a VRT file from my original BSB file, I
get the following warning:

E:\ms4w\apps\marinfo\data\BSB\KAP>gdalbuildvrt 131701.vrt 131701.kap
0...10...20...30...40...50...60...70...80...90...100 - done.
Warning 6: gdalbuildvrt does not support rotated geo transforms. Skipping
131701
.kap

Any help on this will be appreciated.

Regards,

Marc-André



--
View this message in context: http://osgeo-org.1560.n6.nabble.com/Projection-issue-with-vrt-file-tp4903846p4934316.html
Sent from the Mapserver - User mailing list archive at Nabble.com.

More information about the mapserver-users mailing list