[mapserver-users] How to project a shapefile ?

MOREL Stephane stephane.morel at vanoise.com
Thu Oct 3 04:55:05 EDT 2002


Hi Thank for your help
but it still doesn't work
I have this error message :
Fatal error: MapServer Error in msProcessProjection(): projection not named 
in .../../my php files
I saw you load a dll called php_proj.dll. I don't if we need to do the same 
with Linux with something like dl('php_proj.so')
Here is my code
.....
$oLine->addXY($xmin, $ymin);
$oldproj = ms_newprojectionobj("proj=latlong");
$newproj = ms_newprojectionobj("<27585>proj=lcc a=6378249.2 b=6356515.0 
lat_0=46d48'0.0''N lon_0=2d20'14.025''E lat_1=45d53'56.108''N 
lat_2=47d41'45.652''N x_0=600000.0 y_0=2200000.0 towgs84=-168,-60,+320 
units=m no_defs");
$oLine = $oLine->project($oldproj,$newproj);
$oShp->add( $oLine );
...

if I omit <27585> in the $newproj variable, I get another error message :
Fatal error: MapServer Error in msProcessProjection(): unknown projection 
id in

Stephane MOREL
GIS manager
www.vanoise.com

Hi

   This script creates shape file and reproject in to polyconic. U want to
chane in French projection, change the values in $proj2 variable.

Let me know if u have any quetions.

Good luck,

Anwar
www.atsincorp.com

<?php

dl("php_mapscript.dll");
dl("php_dbase.dll");
dl("php_proj.dll");

function convert($n)
{
   $t = explode(' ',$n);

   $d  = trim($t[0]);
   $m  = trim($t[1]);
   $ds = trim($t[2]);

   $dnum  = (int) $d;
   $mnum  = (int) $m;
   $dsnum = (int) $ds;

   $dd = $dnum + ($mnum / 60) + ($dsnum / 3600);
   if ($dnum > 43)
   { $dd = 0 - $dd; }

   //echo $dnum,"  ,  ",$mnum,"  ,  ",$dsnum,"<BR>";
   return $dd;

}

# End of the function


# Function to generate SHAPE Files

function createShape( $x1, $y1, $a1, $b1, $scale, $programId, $pname )
{
    GLOBAL $shpFile,$dbfFile,$proj1,$proj2;

    $x = convert($x1);
    $y = convert($y1);
    $a = convert($a1);
    $b = convert($b1);

    $oShp = ms_newShapeObj(MS_SHP_POLYGON);
    $oLine = ms_newLineObj();

    $opoint = ms_newPointObj();
    $opoint->setXY($x , $y );

    $opoint1 = ms_newPointObj();
    $opoint1->setXY($a , $b );

    $mmwidth = $scale * 228;
    $mwidth = $mmwidth / 1000;
    $swidth = $mwidth / 23.4;
    $dwidth = $swidth / 3600;

    $endx = $x ;
    $endy = $y ;
    $startx = $a ;
    $starty = $b ;

    $indy = $endy - $starty;
    $indx = $endx - $startx;
    echo $indx,' , ', $indy,'  , ';

    if ($indy == 0)
	{
	    $dy = ($dwidth / 2);
		$p1 = ms_newPointObj();
		$p1->setXY( $endx, $endy+$dy );
		$p1=$p1->project($proj1, $proj2);
		$oLine->add($p1);

		$p2 = ms_newPointObj();
		$p2->setXY( $endx, $endy-$dy );
		$p2=$p2->project($proj1, $proj2);
		$oLine->add($p2);

		$p3 = ms_newPointObj();
		$p3->setXY( $startx, $starty-$dy );
		$p3=$p3->project($proj1, $proj2);
		$oLine->add($p3);

		$p4 = ms_newPointObj();
		$p4->setXY( $startx, $starty+$dy );
		$p4=$p4->project($proj1, $proj2);
		$oLine->add($p4);

		$p1 = ms_newPointObj();
		$p1->setXY( $endx, $endy+$dy );
		$p1=$p1->project($proj1, $proj2);
		$oLine->add($p1);
	}

    elseif ($indx == 0)
    {
 		 $dx = ($dwidth / 2);
		 $p1 = ms_newPointObj();
		 $p1->setXY( $endx-$dx, $endy );
		 $p1=$p1->project($proj1, $proj2);
		 $oLine->add($p1);

		 $p2 = ms_newPointObj();
		 $p2->setXY( $endx+$dx, $endy );
		 $p2=$p2->project($proj1, $proj2);
		 $oLine->add($p2);

		 $p3 = ms_newPointObj();
		 $p3->setXY( $startx+$dx, $starty );
		 $p3=$p3->project($proj1, $proj2);
		 $oLine->add($p3);

		 $p4 = ms_newPointObj();
		 $p4->setXY( $startx-$dx, $starty );
		 $p4=$p4->project($proj1, $proj2);
		 $oLine->add($p4);

		 $p1 = ms_newPointObj();
		 $p1->setXY( $endx-$dx, $endy );
		 $p1=$p1->project($proj1, $proj2);
		 $oLine->add($p1);

     }
    else
	{
		$inz = $opoint->distanceToPoint( $opoint1 );
		$dx = ($dwidth / 2) * ($indy / $inz);
		$dy = ($dwidth / 2) * ($indx / $inz);

		$p1 = ms_newPointObj();
		$p1->setXY( $endx-$dx, $endy+$dy );
		$p1=$p1->project($proj1, $proj2);
		$oLine->add($p1);

		$p2 = ms_newPointObj();
		$p2->setXY( $endx+$dx, $endy-$dy );
		$p2=$p2->project($proj1, $proj2);
		$oLine->add($p2);

		$p3 = ms_newPointObj();
		$p3->setXY( $startx+$dx, $starty-$dy );
		$p3=$p3->project($proj1, $proj2);
		$oLine->add($p3);

		$p4 = ms_newPointObj();
		$p4->setXY( $startx-$dx, $starty+$dy );
		$p4=$p4->project($proj1, $proj2);
		$oLine->add($p4);

		$p1 = ms_newPointObj();
		$p1->setXY( $endx-$dx, $endy+$dy );
		$p1=$p1->project($proj1, $proj2);
		$oLine->add($p1);
	}

    $oShp->add($oLine);

    $shpFile->addShape($oShp);
    dbase_add_record($dbfFile, array($programId, $pname ));

}


# End of Function generate SHAPE Files


$shpFname = 'E:\Program Files\Apache
Group\Apache\htdocs\php406\pams\flightline';
$shpFile = ms_newShapeFileObj( $shpFname, MS_SHP_POLYGON);
$dbfFile = dbase_create($shpFname.".dbf",
array(array("PROG_ID","N",5,0),array("PROJ_NAME","C",10)));

$proj1  = ms_newprojectionobj("proj=latlong,ellps=GRS80");
$proj2  =
ms_newprojectionobj("proj=poly,ellps=GRS80,lat_0=40.925n,lon_0=77.75w");

createShape( "-076 17 34","40 03 59","-076 18 34","40 04 59",3000, 1,
"PAMS");

echo "Shapes Created.<BR>";

$shpFile->free();
echo "Shape File ($shpFname) closed.<BR>";

dbase_close($dbfFile);
echo "Dbase file closed. <br>";

?>


-----Original Message-----
From: MOREL Stephane [mailto:stephane.morel at vanoise.com]
Sent: Wednesday, October 02, 2002 12:52 PM
To: 'mapserver-users at lists.gis.umn.edu'
Subject: [mapserver-users] How to project a shapefile ?


Hi list
I'm a newbie on Mapserver
I use Mapserver 3.6.1 / Suse 7.2
MapServer version 3.6.1 OUTPUT=PNG OUTPUT=JPEG OUTPUT=WBMP SUPPORTS=PROJ
SUPPORTS=TTF SUPPORTS=WMS_SERVER INPUT=TIFF INPUT=EPPL7 INPUT=JPEG
INPUT=OGR INPUT=GDAL INPUT=SHAPEFILE
I wrote a small php script to create a shapefile filled with datas
retrieved from a relational database.
Differents users can query the database and i would like to allow them to
cartography their results through Mapserver.
The datas in the database is in latlong projection and all the rest of my
geographic database are in Lambert II etendu which is a french projection.
I found in the archive a post from Claude Philipona about this french
projection parameters.
I've updated my epsg file with those parameters.
The script works perfectly but I still have shapefiles in latlon projection
how can I transform the created shapefile in the french projection ?
In fact I don't know how to use projectionobj with Mapscript

Any help will be appreciated

TIA

Stephane MOREL
GIS manager
www.vanoise.com



More information about the mapserver-users mailing list