[postgis-devel] GeomRound and GeomTrunc functions
Gino Lucrezi
gino-postgisdevel at lucrezi.net
Fri Jul 2 02:10:12 PDT 2004
They are working, now.
In the end, I decided that these functions are only meant for simplifying display of points or structures, and not for further processing, so I haven't yet implemented any check on identical consecutive points and the such.
That would have involved more pointer handling than I liked, maybe I'll do it in the future.
I left a warning to that effect in the code.
I aree with strk's thoughts on the subject, though.
The number of digits might be negative:
GeomRound( GeomFromText('POINT(358808.488 701948.639)',23033), -5 )
will result in POINT(400000 4700000), even though this is a rather extreme example.
Since the number of digits can't be too high, I used an int2 as parameter type, so I had to do some SQL type casting. I don't know if it's a good idea.
Also, I left this parameter optional (default value being 0), and I did the overloading in SQL; it might have been better to do it in the C code, I don't know.
Since this was my first attempt at hacking into PostGIS I didn't know what the standard practices are.
The code is shamelessly based on the translate function.
I am enclosing diffs with existing files, it's just a matter of adding a few functions.
If there is a better way to submit code, tell me. The server I develop on has no internet access, so I can't connect to a CVS server.
Gino Lucrezi
Penta Consulting Services
diff ./postgis_fn.c ../postgis-cvs-originale/postgis_fn.c
985,1203d984
< //round geometry
< // rounds off all points in a geometry up to ndigits
< // WARNING: might generate invalid geometries.
< // It is only meant for simplifying display of points or structures.
< // Different points in a line or polygon might "collapse" into the same
< // coordinates. A possible improvement could be removing consecutive
< // points with the same coordinates, but lines and polygons might
< // degenerate into single points.
< // Obviously, this function is not a OGC standard.
< // It was inspired by the SQL ROUND() function (which only works on scalar data)
< // ndigits might be negative:
< // GeomRound( GeomFromText('POINT(358808.488 4701948.639)',23033), -5 )
< // will result in POINT(400000 4700000)
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
<
< PG_FUNCTION_INFO_V1(geom_round);
< Datum geom_round(PG_FUNCTION_ARGS)
< {
< GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
< GEOMETRY *geom1;
< int2 ndigits = PG_GETARG_INT16(1);
< int32 *offsets1;
< char *o1;
< int32 type,j,i;
< POLYGON3D *poly;
< POINT3D *point,*pts;
< LINE3D *line;
< int numb_points;
<
< //make a copy of geom so we can return a new version
< geom1 = (GEOMETRY *) palloc (geom->size);
< memcpy(geom1, geom, geom->size); //Will handle SRID and grid
<
< offsets1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs ) ;
<
< //now have to do a scan of each object
<
< for (j=0; j< geom1->nobjs; j++) //for each object in geom1
< {
< o1 = (char *) geom1 +offsets1[j] ;
< type= geom1->objType[j];
<
< if (type == POINTTYPE) //point
< {
< point = (POINT3D *) o1;
< round_points(point, 1,ndigits);
<
< }
<
< if (type == LINETYPE) //line
< {
< line = (LINE3D *) o1;
< round_points(&(line->points[0]), line->npoints,ndigits);
< }
<
< if (type == POLYGONTYPE) //POLYGON
< {
< poly = (POLYGON3D *) o1;
< //find where the points are and where to put them
< numb_points =0;
< for (i=0; i<poly->nrings;i++)
< {
< numb_points += poly->npoints[i];
< }
< pts = (POINT3D *) ( (char *)&(poly->npoints[poly->nrings] ) );
< pts = (POINT3D *) MAXALIGN(pts);
< round_points(pts, numb_points,ndigits);
< }
< }
< //round the bounding box as well
< geom1->bvol.LLB.x = my_round( geom1->bvol.LLB.x, ndigits );
< geom1->bvol.LLB.y = my_round( geom1->bvol.LLB.y, ndigits );
< geom1->bvol.LLB.z = my_round( geom1->bvol.LLB.z, ndigits );
< geom1->bvol.URT.x = my_round( geom1->bvol.URT.x, ndigits );
< geom1->bvol.URT.y = my_round( geom1->bvol.URT.y, ndigits );
< geom1->bvol.URT.z = my_round( geom1->bvol.URT.z, ndigits );
<
< PG_RETURN_POINTER(geom1);
< }
<
<
< // Round a single point to ndigits decimal places
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
<
< void round_points(POINT3D *pt, int npoints, int2 ndigits)
< {
< int i;
<
< if (npoints <1)
< return; //nothing to do
<
< for (i=0;i<npoints;i++)
< {
< //elog(LOG,"before: [%g,%g,%g]\n", pt[i].x, pt[i].y,pt[i].z);
<
< pt[i].x = my_round( pt[i].x, ndigits );
< pt[i].y = my_round( pt[i].y, ndigits );
< pt[i].z = my_round( pt[i].z, ndigits );
<
< //elog(LOG,"after: [%g,%g,%g]\n", pt[i].x, pt[i].y,pt[i].z);
< }
< }
<
< // round a double to a specified number of digits
< // ndigits can be negative
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
< double my_round( double X, int2 ndigits )
< {
< return round( X * pow(10, ndigits) ) * pow(10, -ndigits );
< }
<
<
< //trunc geometry
< // truncates extra digits off all points in a geometry up to ndigits
< // WARNING: might generate invalid geometries.
< // It is only meant for simplifying display of points or structures.
< // See geom_round for details
< // Obviously, this function is not a OGC standard.
< // It was inspired by the SQL TRUNC() function (which only works on scalar data)
< // ndigits might be negative:
< // GeomTrunc( GeomFromText('POINT(358808.488 4701948.639)',23033), -5 )
< // will result in POINT(300000 4700000)
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
<
< PG_FUNCTION_INFO_V1(geom_trunc);
< Datum geom_trunc(PG_FUNCTION_ARGS)
< {
< GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
< GEOMETRY *geom1;
< int2 ndigits = PG_GETARG_INT16(1);
< int32 *offsets1;
< char *o1;
< int32 type,j,i;
< POLYGON3D *poly;
< POINT3D *point,*pts;
< LINE3D *line;
< int numb_points;
<
< //make a copy of geom so we can return a new version
< geom1 = (GEOMETRY *) palloc (geom->size);
< memcpy(geom1, geom, geom->size); //Will handle SRID and grid
<
< offsets1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs ) ;
<
< //now have to do a scan of each object
<
< for (j=0; j< geom1->nobjs; j++) //for each object in geom1
< {
< o1 = (char *) geom1 +offsets1[j] ;
< type= geom1->objType[j];
<
< if (type == POINTTYPE) //point
< {
< point = (POINT3D *) o1;
< trunc_points(point, 1,ndigits);
<
< }
<
< if (type == LINETYPE) //line
< {
< line = (LINE3D *) o1;
< trunc_points(&(line->points[0]), line->npoints,ndigits);
< }
<
< if (type == POLYGONTYPE) //POLYGON
< {
< poly = (POLYGON3D *) o1;
< //find where the points are and where to put them
< numb_points =0;
< for (i=0; i<poly->nrings;i++)
< {
< numb_points += poly->npoints[i];
< }
< pts = (POINT3D *) ( (char *)&(poly->npoints[poly->nrings] ) );
< pts = (POINT3D *) MAXALIGN(pts);
< trunc_points(pts, numb_points,ndigits);
< }
< }
< //truncates the bounding box as well
< geom1->bvol.LLB.x = my_trunc( geom1->bvol.LLB.x, ndigits );
< geom1->bvol.LLB.y = my_trunc( geom1->bvol.LLB.y, ndigits );
< geom1->bvol.LLB.z = my_trunc( geom1->bvol.LLB.z, ndigits );
< geom1->bvol.URT.x = my_trunc( geom1->bvol.URT.x, ndigits );
< geom1->bvol.URT.y = my_trunc( geom1->bvol.URT.y, ndigits );
< geom1->bvol.URT.z = my_trunc( geom1->bvol.URT.z, ndigits );
<
< PG_RETURN_POINTER(geom1);
< }
<
<
< // Truncate a single point to ndigits decimal places
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
<
< void trunc_points(POINT3D *pt, int npoints, int2 ndigits)
< {
< int i;
<
< if (npoints <1)
< return; //nothing to do
<
< for (i=0;i<npoints;i++)
< {
< //elog(LOG,"before: [%g,%g,%g]\n", pt[i].x, pt[i].y,pt[i].z);
<
< pt[i].x = my_trunc( pt[i].x, ndigits );
< pt[i].y = my_trunc( pt[i].y, ndigits );
< pt[i].z = my_trunc( pt[i].z, ndigits );
<
< //elog(LOG,"after: [%g,%g,%g]\n", pt[i].x, pt[i].y,pt[i].z);
< }
< }
<
< // truncate a double to a specified number of digits
< // ndigits can be negative
< // 2004/07/01 Gino Lucrezi (gino-pgsc at lucrezi.net)
< double my_trunc( double X, int2 ndigits )
< {
< return trunc( X * pow(10, ndigits) ) * pow(10, -ndigits );
< }
diff ./postgis.h ../postgis-cvs-originale/postgis.h
418,421d417
< void round_points(POINT3D *pt, int npoints, int2 ndigits);
< void trunc_points(POINT3D *pt, int npoints, int2 ndigits);
< double my_round( double X, int2 ndigits );
< double my_trunc( double X, int2 ndigits );
545,546d540
< Datum geom_round(PG_FUNCTION_ARGS);
< Datum geom_trunc(PG_FUNCTION_ARGS);
diff ./postgis_sql_common.sql.in ../postgis-cvs-originale/postgis_sql_common.sql.in
509,538d508
< CREATE FUNCTION GeomRound(geometry,smallint)
< RETURNS geometry
< AS '@MODULE_FILENAME@', 'geom_round'
< LANGUAGE 'C' WITH (isstrict) ;
<
< CREATE FUNCTION GeomRound(geometry, numeric)
< RETURNS geometry AS
< 'SELECT GeomRound( $1, $2::int2)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATE FUNCTION GeomRound(geometry)
< RETURNS geometry AS
< 'SELECT GeomRound( $1, 0)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATE FUNCTION GeomTrunc(geometry,smallint)
< RETURNS geometry
< AS '@MODULE_FILENAME@', 'geom_trunc'
< LANGUAGE 'C' WITH (isstrict) ;
<
< CREATE FUNCTION GeomTrunc(geometry, numeric)
< RETURNS geometry AS
< 'SELECT GeomTrunc( $1, $2::int2)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATE FUNCTION GeomTrunc(geometry)
< RETURNS geometry AS
< 'SELECT GeomTrunc( $1, 0)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
diff ./postgis.sql.in ../postgis-cvs-originale/postgis.sql.in
956,985d955
< CREATEFUNCTION GeomRound(geometry,smallint)
< RETURNS geometry
< AS '@MODULE_FILENAME@', 'geom_round'
< LANGUAGE 'C' WITH (isstrict) ;
<
< CREATEFUNCTION GeomRound(geometry, numeric)
< RETURNS geometry AS
< 'SELECT GeomRound( $1, $2::int2)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATEFUNCTION GeomRound(geometry)
< RETURNS geometry AS
< 'SELECT GeomRound( $1, 0)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATEFUNCTION GeomTrunc(geometry,smallint)
< RETURNS geometry
< AS '@MODULE_FILENAME@', 'geom_trunc'
< LANGUAGE 'C' WITH (isstrict) ;
<
< CREATEFUNCTION GeomTrunc(geometry, numeric)
< RETURNS geometry AS
< 'SELECT GeomTrunc( $1, $2::int2)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
< CREATEFUNCTION GeomTrunc(geometry)
< RETURNS geometry AS
< 'SELECT GeomTrunc( $1, 0)'
< LANGUAGE 'sql' WITH (isstrict) ;
<
More information about the postgis-devel
mailing list