From fc0aa78676ce7ddb4634f9d78c3b444b5087f5c5 Mon Sep 17 00:00:00 2001 From: No Body Date: Sat, 21 Jul 2001 00:43:33 +0000 Subject: [PATCH] This commit was manufactured by cvs2svn to create tag 'postgis_0_5'. git-svn-id: http://svn.osgeo.org/postgis/tags/postgis_0_5@33 b70326c6-7e19-0410-871a-916f4a2858ee --- CHANGES | 36 + CREDITS | 8 +- Makefile | 20 +- README.postgis | 4 +- TODO | 11 +- doc/postgis.xml | 227 +++-- examples/wkb_reader/Makefile | 3 + examples/wkb_reader/comp | 2 +- examples/wkb_reader/readwkb | Bin 21420 -> 35744 bytes examples/wkb_reader/readwkb.c | 12 +- loader/Makefile | 11 + loader/README | 63 ++ loader/dbfopen.c | 1074 ++++++++++++++++++++++ loader/shapefil.h | 351 +++++++ loader/shp2pgsql.c | 858 +++++++++++++++++ loader/shpopen.c | 1623 +++++++++++++++++++++++++++++++++ postgis.h | 69 +- postgis.sql.in | 183 +++- postgis_debug.c | 9 + postgis_fn.c | 534 ++++++++++- postgis_inout.c | 302 +++++- postgis_proj.c | 338 +++++++ 22 files changed, 5572 insertions(+), 166 deletions(-) create mode 100644 loader/Makefile create mode 100644 loader/README create mode 100644 loader/dbfopen.c create mode 100644 loader/shapefil.h create mode 100644 loader/shp2pgsql.c create mode 100644 loader/shpopen.c create mode 100644 postgis_proj.c diff --git a/CHANGES b/CHANGES index 04a7552aa..cf2ae25b7 100644 --- a/CHANGES +++ b/CHANGES @@ -1,3 +1,39 @@ +PostGIS 0.5 + +- New functions + - Dimension() + - GeometryType() + - Envelope() + - X(), Y(), Z() + - NumPoints() + - PointN() + - ExteriorRing() + - NumInteriorRings() + - InteriorRingN() + - NumGeometries() + - GeometryN() + - Length_Spheroid() + - Length3D_Spheroid() + - AsBinary() + XDR and NDR variants + - force_collection() +- Removed functions + - wkb_ndr() + - wkb_xdr() +- New Objects + - SPHEROID(,,) + To be used with the length_spheroid functions for accurate + length calculations on lat/lon data. +- Minor bug fixes +- Internal Functions + - Extra constructors to make geometry manipulation easier +- Structural Reorganization + - Broke postgis.c up into four new files + postgis_debug.c -- debugging functions + postgis_fn.c -- generic functions (like length()) + postgis_ops.c -- operators and indexing functions + postgis_inout.c -- type support functions and data conversion functions + + PostGIS 0.2 2001/06/19 diff --git a/CREDITS b/CREDITS index f0b535b06..df1c609b3 100644 --- a/CREDITS +++ b/CREDITS @@ -8,5 +8,11 @@ participated in the project: Dave Blasby - Core server objects and indexing (the hard stuff!) Paul Ramsey - Extensions to the PostgreSQL JDBC driver. Jeff Lounsbury - Shape file loader/dumper. - + +Version 0.2 of PostGIS includes contributions from: + Norman Vine on CygWin compilation issues + +Version 0.5 of PostGIS includes contributions from: + Geographic Data BC (David Skea and Mark Sondheim) with funding and + direction on adding calculations on a spheroid to PostGIS. diff --git a/Makefile b/Makefile index 1af4bd05d..8bf4f27ee 100644 --- a/Makefile +++ b/Makefile @@ -4,16 +4,18 @@ subdir = contrib/postgis # Root of the pgsql source tree -top_builddir = ../.. -#top_builddir = /data3/postgresql-7.1.2/src +ifeq (${PGSQL_SRC},) + top_builddir = ../.. + installlibdir = $(libdir)/contrib +else + top_builddir = ${PGSQL_SRC} + installlibdir = ${PWD} +endif include $(top_builddir)/src/Makefile.global test_db = geom_regress -# override libdir to install shlib in contrib not main directory -libdir := $(libdir)/contrib - # shared library parameters NAME=postgis SO_MAJOR_VERSION=0 @@ -24,7 +26,7 @@ SO_MINOR_VERSION=3 override CPPFLAGS := -g -I$(srcdir) $(CPPFLAGS) -DFRONTEND -DSYSCONFDIR='"$(sysconfdir)"' override DLLLIBS := $(BE_DLLLIBS) $(DLLLIBS) -OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o +OBJS=postgis_debug.o postgis_ops.o postgis_fn.o postgis_inout.o postgis_proj.o # Add libraries that libpq depends (or might depend) on into the # shared library link. (The order in which you list them here doesn't @@ -37,17 +39,15 @@ all: all-lib $(NAME).sql include $(top_srcdir)/src/Makefile.shlib $(NAME).sql: $(NAME).sql.in - sed -e 's:@MODULE_FILENAME@:$(libdir)/$(shlib):g' < $< > $@ + sed -e 's:@MODULE_FILENAME@:$(installlibdir)/$(shlib):g' < $< > $@ -#$(NAME).sql: $(NAME).sql.in -# sed -e 's:@MODULE_FILENAME@:/data1/Refractions/Projects/PostGIS/geom/$(shlib):g' < $< > $@ install: all installdirs install-lib $(INSTALL_DATA) $(srcdir)/README.$(NAME) $(docdir)/contrib $(INSTALL_DATA) $(NAME).sql $(datadir)/contrib installdirs: - $(mkinstalldirs) $(docdir)/contrib $(datadir)/contrib $(libdir) + $(mkinstalldirs) $(docdir)/contrib $(datadir)/contrib $(installlibdir) uninstall: uninstall-lib @rm -f $(docdir)/contrib/README.$(NAME) $(datadir)/contrib/$(NAME).sql diff --git a/README.postgis b/README.postgis index 7b0746b35..4c5b6508c 100644 --- a/README.postgis +++ b/README.postgis @@ -1,7 +1,7 @@ PostGIS - Geographic Information Systems Extensions to PostgreSQL ~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -VERSION: 0.2 (2001/06/19) +VERSION: 0.5 (2001/07/20) MORE INFORMATION: http://postgis.refractions.net @@ -17,6 +17,8 @@ Directory structure: the GIS objects. ./doc Documentation on the code, objects and functions provided. + ./loader A program to convert ESRI Shape files into SQL text + suitable for uploading into a PostGIS/PostgreSQL database. ./examples Small programs which demonstrate ways of accessing GIS data. diff --git a/TODO b/TODO index 56c3d5906..0dc303da4 100644 --- a/TODO +++ b/TODO @@ -1,8 +1,7 @@ -2001/06/19 +2001/07/20 -- Change function names to be OGC compliant (following the "Simple Features - Specification for SQL") -- Add some more of the OGC spec functions -- Tie in to simple GIS clients - GML import/export routines -- Document Java code better. Usage example better and perhaps client app. +- Possible re-write to become more OGC compliant (?) +- Inclusion of OGC simple feature predicate functions for overlaps, + intersections, etc. +- Connectivity to commercial GIS software. diff --git a/doc/postgis.xml b/doc/postgis.xml index 08dd80a57..806c11ecb 100644 --- a/doc/postgis.xml +++ b/doc/postgis.xml @@ -40,10 +40,6 @@ functions. -To Do List -This is the To Do list. - - More Information The latest software, documentation and news items are available at the PostGIS web site, http://postgis.refractions.net. @@ -61,18 +57,18 @@ functions. The PostGIS module is a extension to the PostgreSQL backend server. As such, PostGIS requires a full copy of the PostgreSQL source tree in order to compile. The PostgreSQL source code is available athttp://www.postgresql.org. -PostGIS 0.2 has been built and tested against PostgreSQL 7.1.2. It +PostGIS 0.5 has been built and tested against PostgreSQL 7.1.2. It will probably work with any of the 7.1.x versions. It will not work with any version prior to 7.1. Before you can compile the postgis server modules, you must compile and install the PostgreSQL package. -Retrieve the PostGIS source archive from http://postgis.refractions.net/postgis-0.2.tar.gz. Uncompress and untar the +Retrieve the PostGIS source archive from http://postgis.refractions.net/postgis-0.5.tar.gz. Uncompress and untar the archive in the "contrib" directory of the PostgreSQL source tree.cd [postgresql source tree]/contrib -zcat postgis-0.2.tar.gz | tar xvf - +zcat postgis-0.5.tar.gz | tar xvf - Once your PostgreSQL installation is up-to-date, enter the - "postgis" directory, and run the compile and install commands. cd ./postgis-0.2 + "postgis" directory, and run the compile and install commands. cd ./postgis-0.5 make make install @@ -102,7 +98,8 @@ make install Loader/Dumper - No documentation yet. + The data loader should compile very simple.cd postgis-0.5/loader +makeThe loader is called "shp2pgsql" and converts ESRI Shape files into SQL suitable for loading in PostGIS/PostgreSQL. @@ -300,8 +297,8 @@ INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES (6,'LINESTRING(198231 263418,19821 Using the Loader - This section to be done. - + The data loader converts ESRI Shape files into SQL suitable for insertion into a PostGIS/PostgreSQL database. The loader has several operating modes distinguished by command line flags:-dDrops the database table before creating a new table with the data in the Shape file.-aAppends data from the Shape file into the database table. Note that to use this option to load multiple files, the files must have the same attributes and same data types.-cCreates a new table and populates it from the Shape file. This is the default mode.-dumpCreates a new table and populates it from the Shape file. This uses the PostgreSQL "dump" format for the output data and is much faster to load than the default "insert" SQL format. Use this for very large data sets.An example session using the loader to create an input file and uploading it might look like this:shp2pgsql shaperoads roads.sql +psql -d roadsdb -f roads.sqlA conversion and upload can be done all in one step using UNIX pipes:shp2pgsql shaperoads | psql -d roadsdb Retrieving GIS Data @@ -536,44 +533,81 @@ if( geom.getType() = Geometry.POLYGON ) { PostGIS Reference The functions given below are the ones which a user of PostGIS is likely to need. There are other functions which are required support functions to the PostGIS objects which are not of use to a general user. - Functions + OpenGIS Functions - length3d(geometry) - Returns the 3-dimensional length of the geometry if it is a - linestring or multi-linestring. + AsBinary(geometry) + Returns the geometry in the OGC "well-known-binary" format, + using the endian encoding of the server on which the database is running. This is useful in binary cursors to pull data out + of the database without converting it to a string representation. + + Dimension(geometry) + Returns '2' if the geometry is two dimensional and '3' if the geometry is three dimensional. + - length2d(geometry) - Returns the 23-dimensional length of the geometry if it is - a linestring or multi-linestring. - + Envelope(geometry) + Returns a POLYGON representing the bounding box of the geometry. - area2d(geometry) - Returns the area of the geometry if it is a polygon or - multi-polygon. - + GeometryType(geometry) + Returns the type of the geometry. Eg: LINESTRING, POLYGON, MULTIPOINT, etc. + + X(geometry) + Find and return the X coordinate of the first point in the geometry. Return NULL if there is no point in the geometry. - perimeter3d(geometry) - Returns the 3-dimensional perimeter of the geometry, if it - is a polygon or multi-polygon. - + Y(geometry) + Find and return the Y coordinate of the first point in the geometry. Return NULL if there is no point in the geometry. - perimeter2d(geometry) - Returns the 2-dimensional perimeter of the geometry, if it - is a polygon or multi-polygon. - + Z(geometry) + Find and return the Z coordinate of the first point in the geometry. Return NULL if there is no point in the geometry. - truly_inside(geometryA,geometryB) - Returns true if any part of A is within the bounding box - of B. - + NumPoints(geometry) + Find and return the number of points in the first linestring in the geometry. Return NULL if there is no linestring in the geometry. + + + PointN(geometry,integer) + Return the N'th point in the first linestring in the geometry. Return NULL if there is no linestring in the geometry. + + + ExteriorRing(geometry) + Return the exterior ring of the first polygon in the geometry. Return NULL if there is no polygon in the geometry. + + NumInteriorRings(geometry) + Return the number of interior rings of the first polygon in the geometry. Return NULL if there is no polygon in the geometry. + + + InteriorRingN(geometry,integer) + Return the N'th interior ring of the first polygon in the geometry. Return NULL if there is no polygon in the geometry. + + + NumGeometries(geometry) + If geometry is a GEOMETRYCOLLECTION return the number of geometries, otherwise return NULL. + + + GeometryN(geometry) + Return the N'th geometry if the geometry is a GEOMETRYCOLLECTION. Otherwise, return NULL. + + + + + + + + + + + + + Other Functions + + + A <& B The "<&" operator returns true if A's bounding box @@ -618,34 +652,28 @@ if( geom.getType() = Geometry.POLYGON ) { - npoints(geometry) - Returns the number of points in the geometry. - - - - nrings(geometry) - If the geometry is a polygon or multi-polygon returns the - number of rings. + area2d(geometry) + Returns the area of the geometry if it is a polygon or + multi-polygon. - mem_size(geometry) - Returns the amount of space (in bytes) the geometry takes. - + asbinary(geometry,'NDR') + Returns the geometry in the OGC "well-known-binary" format, + using little-endian encoding. This is useful in binary cursors to pull data out + of the database without converting it to a string representation. - - - numb_sub_objects(geometry) - Returns the number of objects stored in the geometry. This - is useful for MULTI-geometries and GEOMETRYCOLLECTIONs. + + asbinary(geometry,'XDR') + Returns the geometry in the OGC "well-known-binary" format, + using big-endian encoding. This is useful in binary cursors to pull data out + of the database without converting it to a string representation. - - - summary(geometry) - Returns a text summary of the contents of the geometry. + + box3d(geometry) + Returns a BOX3D representing the maximum extents of the geometry. - - + extent(geometry) The extent() function is an "aggregate" function in the terminology of PostgreSQL. That means that it operators on lists of data, in @@ -654,14 +682,11 @@ if( geom.getType() = Geometry.POLYGON ) { all features in the table. Similarly, "SELECT EXTENT(GEOM) FROM GEOMTABLE GROUP BY CATEGORY" will return one extent result for each category. - - - translate(geometry,float8,float8,float8) - Translates the geometry to a new location using the numeric - parameters as offsets. Ie: translate(geom,X,Y,Z). + + force_collection(geometry) + Converts the geometry into a GEOMETRYCOLLECTION. This is useful for simplifying the WKB representation. - - + force_2d(geometry) Forces the geometries into a "2-dimensional mode" so that all output representations will only have the X and Y coordinates. This is @@ -675,23 +700,73 @@ if( geom.getType() = Geometry.POLYGON ) { all output representations will have the X, Y and Z coordinates. + - wkb_ndr(geometry) - Returns the geometry in the OGC "well-known-binary" format, - using little-endian encoding. This is useful in binary cursors to pull data out - of the database without converting it to a string representation. + length2d(geometry) + Returns the 23-dimensional length of the geometry if it is + a linestring or multi-linestring. + + + length3d(geometry) + Returns the 3-dimensional length of the geometry if it is a + linestring or multi-linestring. + + + length_spheroid(geometry,spheroid) + Calculates the length of of a geometry on an elipsoid. This is useful if the coordinates of the geometry are in latitude/longitude and a length is desired without reprojection. The elipsoid is a separate database type and can be constructed as follows:SPHEROID[<NAME>,<SEMI-MAJOR AXIS>,<INVERSE FLATTENING>] +Eg: SPHEROID["GRS_1980",6378137,298.257222101]An example calculation might look like this:SELECT length_spheroid(geometry_column,'SPHEROID["GRS_1980",6378137,298.257222101]') from geometry_table; + + + length3d_spheroid(geometry,spheroid) + Calculates the length of of a geometry on an elipsoid, taking the elevation into account. This is just like length_spheroid except vertical coordinates (expressed in the same units as the spheroid axes) are used to calculate the extra distance vertical displacement adds. + + + mem_size(geometry) + Returns the amount of space (in bytes) the geometry takes. + + + npoints(geometry) + Returns the number of points in the geometry. - wkb_xdr(geometry) - Returns the geometry in the OGC "well-known-binary" - format, using big-endian encoding. This is useful in binary cursors to pull - data out of the database without converting it to a string representation. + nrings(geometry) + If the geometry is a polygon or multi-polygon returns the + number of rings. - - - - + + numb_sub_objects(geometry) + Returns the number of objects stored in the geometry. This + is useful for MULTI-geometries and GEOMETRYCOLLECTIONs. + + + + perimeter2d(geometry) + Returns the 2-dimensional perimeter of the geometry, if it + is a polygon or multi-polygon. + + perimeter3d(geometry) + Returns the 3-dimensional perimeter of the geometry, if it + is a polygon or multi-polygon. + + point_inside_circle(geometry,float,float,float) + The syntax for this functions is point_inside_circle(<geometry>,<circle_center_x>,<circle_center_y>,<radius>). Returns the true if the geometry is a point and is inside the circle. Returns false otherwise. + + + summary(geometry) + Returns a text summary of the contents of the geometry. + + + translate(geometry,float8,float8,float8) + Translates the geometry to a new location using the numeric + parameters as offsets. Ie: translate(geom,X,Y,Z). + + + truly_inside(geometryA,geometryB) + Returns true if any part of A is within the bounding box + of B. + + Programming Information The following sections provide information about the internal structure of the PostGIS spatial objects, and some information about certain design decisions. diff --git a/examples/wkb_reader/Makefile b/examples/wkb_reader/Makefile index eeef03f2c..717386db9 100644 --- a/examples/wkb_reader/Makefile +++ b/examples/wkb_reader/Makefile @@ -13,3 +13,6 @@ CFLAGS=-I$(PG_INC) -L$(PG_LIB) -lpq all: readwkb readwkb: readwkb.c + +clean: + @rm -f readwkb diff --git a/examples/wkb_reader/comp b/examples/wkb_reader/comp index b915e13c9..bf67b16d2 100644 --- a/examples/wkb_reader/comp +++ b/examples/wkb_reader/comp @@ -1,2 +1,2 @@ -gcc -g -I/usr/include/pgsql -I/usr/local/include/pgsql -L/data3/postgresql-7.1.2/lib readwkb.c -o readwkb -lpq +gcc -g -I/data3/postgresql-7.1.2/include -I/usr/include/pgsql -I/usr/local/include/pgsql -L/data3/postgresql-7.1.2/lib readwkb.c -o readwkb -lpq diff --git a/examples/wkb_reader/readwkb b/examples/wkb_reader/readwkb index 6d5e5313c5f6eae63120ae9ce34298c415ff4bf4..659cec97aa7940d3be368f7a9f6f837ab6c0dd7f 100755 GIT binary patch literal 35744 zcmd^odw7)9wfEjh#DF{!5R@v~5xGS}?g@cNjV_RzjQt?82u*XxOMNKW*#GukbIccTM`K`4t^G*n@=lP!J z`{Vmwc=P_&-fQo@_S$Q&eR(hQWl!xQDWxz+mdF-F&D)xjTaIhaCKXdITw;WnAubXZ z;3@>gmA8WjXE!jNXW-1inT_)}ort6J9AG*rj`FspAx_}TpCJUDlwRbi_yW)c&*kLO zS&Rg=IH_EU&%$LV;&&pB&K;oXq4-jb6`rwcTlCa$P8+*s7wG_y4vNv$tT zL<>t*9I|HF7B5=?q{?A{KVxxTjB^am&)_^CCmob09iwr67Uy`JBXM4ca}-X9AtvEG z7pDv71e`+;(ia_6FLSUiKu{jT9VBZm&T%+Nw{)DTF5)~~Nf&3~oQU%RoFxC*I4{CE z+(A4&x98+8#Fg@#j`Jw!rJ63nm2`~@^7_Nz5>(9nc$zMyir$fvdowPRaaQ5X1AQaH zDZc|aU7!nbCHe}S#h^)FMAzVa2<71d`Yh3d>pa9ik1ORzioEuRIk}|ovCw86u3p4n zhbyIDspA1zJ{RcKh@Yg(zY5_o(A0Ls-;A>zbgiae#dRC#GELurtEt~OXl@Mz!b#A9 z)}TL*Ydn}pwI-J*{mE1UWJ@r)*58^!4$WXh018wwk>+r)wTXC|BmOopEUv9vFu&GU zw`kFFPlK;v{(@SMk7AK75dsLt~KUExHZq$}iJ(ed8qBjcA-ml4<8^5WwqLTvtV6H&=C zdK0Kf*ILng7An%a2)SJI%8Q5@tnecB;M4~ECCpeLitF3YG-)A3 zD(vfKnzRvubfWKRrb#Oy7t!09Chde~5WS6Q(o%@pv#*_L(pIR7=or(awa`4GS2Io8 z3oRkq%QR^*dRxAv>b{NEtn>4ht?B)^lu25 zN?H%K6Wzl!X+N};=>1GnD}=TY-OV($M(9DJpJtj`CA6LB?Mzeagmw_U4RpzeU+Vbn z=EDQthWZfs%Mdy7D_4eYdj!625Ddq+WgOh5xtr&fBNvhST>rSoih;VkIUN%Sbc(_K z9J0Aeqz*qf4<2$$cXHI`J;y^;i0N(o)MJqPlr$OQZR5Nqc6xH}%-gzR5ZBHXqwgHGHMi3<`py-DTk|NnuEU-I zM0;^=?mpAo`yzC`c~2e+M94f9Ji^=iq%(MPcb1Ac!`pj*W<<7%$no~xs3J1tRJkga zY;W&X&a}F$EN?H_4@97psva};dQL)3Vy}m$W8x*x?(*XMA(&hU15j(B_T_$z1JzUdQ@9D}ytz(ym)@lZP&zITCQLXkB8 z+3V@0BDZyVdaFH0HscRfo1?VEltlA|Dr=pRqzQOkanl4q;rL_a|Ev6DuhU7n*< z^Uf6~<$iB3O&n0kPR~*4InIAi4tjd~AL!aSjq*o3b*?zRdH2W;G{y1IKT%dYR95}j zuoc^W4a>DE&Uwi*DErUS=|QkQ@K@_&(CT9zv@ziA-Jr5--!y2nfxrfs7Q}kV^NzS> z(?Bn#bZgH7A#d>vh;~G68jxeQL}9cf)k7rJIFrpgV$5ySpGfVnGnf;p{bg@&o>74b zABFI4Z|}e-S~Vto0>Yp3_V%c77HlTU;Y{kRUODYN2ZyYOnBhODwr63;a^0Xy~DfjzS?gxLN+PX+@JBj`^ z(|PTi-Vw>MXzsl{*X|#q?27|vlHE5h4-Vo46rVL_%OVUCl){moS8;3$F!Jggv%Tx zl1wt>-!+Qeq|2q8yQ}Yh)McCM&uo;v9lJaB@9N3v@{a7-pA9$3^}E6;wI|Ip*6^Ta z7dYg#Msuf8jFR^zOAcc+uaR8O46xC_GxXJR&(yBqVGdWub3XZsrN=v7mA{m)c8&fo zZjUSBIQqvaFVgXP$E&;E&FT95$d25;7cmin|A0MZkJ-|U{>_#~u8s`rbt-h8%KQ|= zT4D-y8rCJM4qOt8@mkmDhe?Fa^W|$D6Zd)!P~OV!4pX?X@59xejd3i+W&4=si-o)rphD~{gUS`7=8~4!2F$_ zo@&ogl>d^Nf30Ej=vpY^8UnXUBTxZfC zzt=hr?)q(xGHWugLu6j#RZRvxhx_t*4ALF+`G@k4rk>ngh{$rbWbMCG_S~%#Qj`3?(PnJ}Fy1Zu_ zJBH0yIt8cU%gK1IVMGl_IfS+bztQMO3EgJN(;55|eTZrGyVC*PnE*y&)+}Fy2-Td6{fkAnJjXU zu|Duw`cK3{l(1A>UUg(h9n)kz|`bzG*EqOZge)(f- z-sI1%kA0fM)no_$9Sft;o!1>7#)q>9A9>1LXp)2l2j(`pj{fzAnmrxLn6)24|9TBwyM5CkkwpJGqu*-2xC7V3@6ZWd8I z4h>~Y;xfBmiIEtYjO_3naLnzv_aBFJd${+b3!ceyw6WVrIo$n!_hF`Z|7mEoR(H<- zRE(o=n*S*s@`8D%V^a62=W_v*>yZ9G{B+B>aQOPqrv?3#^SOziDn&;BKf=lx(*JFP z^4PY>K~n!W3O67BIe@N<8bjbh!$(*k=`H6h4{=dLP-hodfEGm%1O-!?S9e3ZF)1U({r-gGkW8Q{&MnX zNX{EmV_wLs&%`IaA=kIz)h<=p4SG)Y@98*^&d-#6zKZYY-aH^vQz3Wg3M6}sZU%XK zH}#ob&vGs@dv4g(1QUmSUa@p^cX&>6o+rJ%Uo&Bw2eOi!=k|1-`I&kCk%tVLXP{!qxA%`G7Wz^lF&mV zRf9>U26-yI((Q=;Dz!f}5OO_XwfnG2Y*M`7AR+%%c5{$5vfad@iy#Xwy>@8nR)6=I z81#Wm{jE-yf1$}z>HKpJ66L?CHp282-RzgJEXUu|wVv%cym3^Y2jjMm$Nv9eZ)iKp z=q*y&IrmRYP2EPi*`6Mx|1+&6wSF)p!R74@a*5bX@h?RgfB{vH%+mzgv98)Y(Zn01Rxellz(FRS>D zEzfZ#^!Ke(GWmhcGqg^#hvonaK+>}*=}Z*|9d~)&@OE+$-rk#3$Dj?EBj}D$cSxywF)CAa=_TpU1`3_hOcTEz9z$m>|xmxlw2j<$^)nGf$_bq@ymaFIc~9S zo=blRlcPJcH@|ptaR27&!PN0Q%{+L~pa zU%ND1jRV)(Ku3C&-@T^k4Sh}LNp+;xJ*qKmAYWCB}H}W#qf(^-_>yoKqSID1m z@xu!)5HqH(hlr7=OFtYDjz(O`sH>?F8&{4a5p2bi6Rx@%x2riGZF9w1eDnZ?Yh4JB zBe($RSqIl#*YtI38mCt~5^5NTrdpdKNf$k3k%AmIq=NAcR1hA6XmUlvIGbH9!Dw4B zN#XQtM#k~?`J2L&cRZL(#Uri~S3J7zVv!Ujvr7w0Dys@BiVMp`1rC8nAmVdl(L}NZ zk8_A*sR#kjotkjre=Hhj*wjc>L10Gt3)Q0|Mmu;YqS2oSx~5M}Om{W=!>z%l!chWg z7I+rdEK`4>(Y7`}YNz5e^_h`52)Acp?fiO=t8GId6;DLtt_3yA=GU)sEnHE*ysqB0 z2&ZlWdM;x|B-NH*;PNLL!x4Xc!;EAo=p*G7Opi3hr{^QWg?3Oar7B|b*ZB<Ixku9E{hjcr~)lk2} zqtK%HwaYb}iibrga%zjdCe&3dtgTy~R;XB5w{&SugQ?RJcm#vtxdsWF>rO}9`?u%h zp1?Wcj-1?caZbZoiE|Ool{lMmuElvf&adJ84$fzBzJ&8vIN!&40_TXWNRM+G&PtqA z9>z+3=7p}^IPqu-KRZUsqh}D{@5Kx}_lFs!$h{4D;5q=B;_2BKdM*w(2ko}dLw8+J}69x zitHCJRs2n$notJ{H@E=!a-4Kf_irbF<9eL!1aRDjvz-8r=Ww915ZF@4S%%*uXf-Y9r&vb{Hz0i?7-t_YemaY>Ab#NY{}9Onv`P@T16|WLPi}Hw!iX^%at3rtkwfLVUzBe3y#h!<3s1 zOKle}4!eNy#f-}sPi9=h_!7pe8Bb-rf$=oP+Zf~NEF3!+monbPxQy{Z#^sEUF|J@d z8vAtZf}TP}7jF&L}GUFiQcE-0c z-pu$h#<5adl@eQCVTULL5htUZ^bi$k24Rm z@E>D)2OQxSgQC#wwcw~0=S9hV@Sa9&4FZp48L5cdDmII0Qng4MPUdvcUij4!gHAT5n~lzHJQRy z-|;h6;kjtScA@Nf8e?U{iy5zGrZ&cke;4C6CVs>C2F7EiP(DewoZ^gCa9^9Y+o`TgQjnM~z zx8r>!OnJq-8e=-e{I}ssJh%uquQmx0p2+WSDSW;IU%>BhsqjhsZkNIr3Ksz!ZnTRF z*M5}M4i2+6??%`NoOJZybm639iaw8le+u4kD~CL1PbZmOEHm-A`OP;%geTj4+~N{` zw@rml5z1Q=&s0&%+=QozpvF{I80UG&e+qWzwxg^PE^CIPtbEN)X)kq@b(y2A0$Ud4 zISc0;&a+q?K%Q4~o<*EzKx5K59d4m~i^3(2{LAdSc5YF@djN{3%Ax0Kt2eipiycOF zxW(tNWohsNQE0K7--9F4Ef(APnOiK;mQLX{;)|C53#=c-Un@4)@TF)=b-3~7+k+Ng zDIT%7-rB2Mtl;+{DatLb<~;+2SBaM`{})A%#n*`UEnX^)SbUv&uMD(Xd`V zagN2S#e^aF2I#z9lnC6t8N8Wsslc6-!Fw2&2|O7!xLsrHnG8O~xLn{i?cfU<;~v}K zON=W7cFPB;I#B zt|zhv&rJ^7+bCQVh{KIq!iDkSM(kz|GPM1@ejobi+?0WWSM+^EH1M6c`?P}7etxGKNYhq-X&@+ zeo@q0yhp6Gc(1tD;#Wkg#jlCA7Vi_=E&jQn$F)_xej$!qykGpq;@8E;x7Faw|7F&FlthIQQthe}V zxzggX@>+|}k*yX_kZUbISAN^#^W>WrPn09IvvZ5{7xDe|ipPnEkYo-X%VJVXA<;!EXw7LSzlZYUjYQ7Fe)Tr4M8Tp~*>E|p6yu8?&W zSIVHpRq`f_Pl!&7tL6O`&z0Y@c%FRH;sx?qix*1T`&6>I<=Yl}_gj38 z{Hw*+$?TCvhA+u+7W?GI7W-wf#f|bBi-U5t#m%zA;uiUD7Ki0?7GE#-THGpsWpS+> zuo&K4O?RQ6n9Q>{E+<)>l(dJc^p=uE7O$1GpQ^&w$*U}0FR!!sM%ildP4ZTYH_8r+ zH_0xGzbqfK_!jwu#kb1uS-eGl-{L#uPb}Ule{OMy{J`Qnp@*|5MlYh1N zTXOVh!}GYDW$_d8a*KCJ+V@rY{Ik5*;_t}^EPhgc$KrpH&sf|g-?aE&hcR{I@@pS1V``Mkyba-YRV$@(hcQ$w?OfNlv!-gq&va zpJjo?pUBHC{;OPM@!w?1;z9X@#UksV#aUVJSe%{pk;OS#f3O4nsLR{@v{S@#G4N1M?o+_qasMtC zadgk+He4y3Sl!?;ha0PUJzYZh3j2J4TdcI}L@KKQ=N`^~m6-Cd*(17|^M8}ey4tZ; zRr_F6h8sHt>HM#C;4eAWr`HM(G2?KH)sFS1!fJi}MP8c|Pd84JXAf{KPC7_VFV5ew zoFS2foPS~2T3Alvaod$D;jqoeja943(=8$n*`j)PpZMcqK4pT#js3dx9>`h|Q^5#d zCvLa+M)6gPH;V6Byh-e`xLxeG_*Se`)Zs>p{kO$i1busl4!2PF%zS0fv|sUA#@DmW z6HhlzQp0fB|@v6m#9qV_jH|<*Z zeKFSJUg5I%15spgzgS}NhmJM)N8(l+{;{~%;y*gpDW|qvdZb9wi4X9wX1oQH|sl zW2MXDadNK3=*Xgu1+_+CNd(L;>y6ExQBV0=H1H9uhd6&`C|VEk3F zA3XoT_-kSxa1Y}Lcntc8@q>ciJ(zGQ`ID~;di;7O<8KIhWO^>+hXg(Kypr+5f*xLu zF@8kQW5-`+{HUOB?{qQVF6gC*M;JdQ=*fJtf^M863E=3#Y4kP)_*SjMCg2fx9!}6F zxe5cKD3yp8CE|gidhgPr1uJT5-9;_*RjZkJ+_SK(D3C}MiIza1)Q3kEWBh0)QYPml zXA~EdkhU_Nu*6Yn8devlOLrh$Y0m2uX=bG} zO=m<^ap*H6uS5vX*;-9hspM7p)icXVXH`~J(_d6wah+1yY{fpCnxnF!qO1a(O0~Ad zN{Y2a#mKXyv`RIb%_+@iOreA_E3H6Q&8<e{jGD@5(|n8%1W!lZT=Qbmz7dLI-x@>(}a~;E|a&j zeA6Pc@~c@u!*hzDTK6STzLIn;9>Z33gq}&Z+4-WduqYe}w5FPZ^tH*r%*2Ma#%ODz zF!YalVPh(cPg91=sw&m@DU0xYd2_f02{W@qUurdh$yO@F8KO#V4&vK-cw*b|snEoR zL{XqA*sO!w!ihkUGmPB?@+*QbfDDaA;R*;tPh`jaYl2#;(!$vl@EN9#_);En%q*uh z8t}IURhYs>HmDV~T05s!80^Jw#;EFwY6^nOxsIuJM$Q2K zg$!UoQa0Ws*2zRFzppD5SF7CQ-F!aZ{Nws^Y4C{-&n5O+pnX^|FxwaN?%mTFbzVPYNQP zg{{IclyJf_3U6_QQ%#W%hle9klTRtt6X7g%TeL~HH`TmXfqHC(D|tt`SSqtsdH z=0v*aGF|I3>itzdA6~^mZ-aI+dVWbY84V2@mF4g=`0lH-iON)US;|001lN~b7m_5y z_&lsJepO+XKe4W{)oQd%RhT0J=;P@KM~$1Cd~N>38r_eSMRV&TA`vu2=o*)!f(Z;k z%^2jUOLOhYQQNd%D?|0D%2iaslHyWT#WleVM%;2`&#KeJ0zOC6m8-gQ??62Q7*VtU z<;tSDDk>)3O6A(9xs{p|LlTk>AP*a&iFBlO+Iw@NIo#?jR5!m)t0Em~wNhAHu>EvQ zh12?BvW^&S{-vB}(wD#l#BaN2g;EXIGKsJHW=1OI=+ct0WICBr4m3)`6S?JN2(!?p z1UYK<+A$SS#BNv<#v3DzVH@h)nLQB*;T5Qe#Pzb{`i5XQNF*9bH5#|X2TEHO}_7}7C z5IRc@p_zPUL>5L2f6{1e7LNhQ1@^%MbyG0u4}?f0o^5ddI{@AOg~YIT-)%iYrF!w*6N$Hrf|HUO7nGgxu)gXe9?6go(egtK+J3K&gzf`Ys=6A zH5t3@46!om%kBuGk!6u&B)22qA$5*M0@{-!iSoHdF7g$rgNWu0bzWSvyl7oCz6P_2wZS5+ zbCUR!d{J<{pY|vcMR;kLTRIpo64c_L&3>C%A-DOHA^3S-T$?b*`aB~XBZgg}XGF9_ zQbmbms?qKvn4EYl38h+sG%Fff3`bJ=YKfkaqm9F8$lDVx8bnC~IAAaC9BnyYMrC^cw|H6ZpTZB+^2QTuWG1 z0Jei)IY{V9A$3+v?NA3j?qy2DhGS+q&gG?9b*hD%V|iJrJ1$Q^Q`FbVJh1b@otCF; zDb}(IRy~#{WX|Pf4hZBh`+p(h>)=a~m78sT<%iEthkbhuq!qM3MtU8!dE>E$IN*d}xW$MO;} zR8MY86k(T<2v^nQlGDpegi|F=73>-X;cB(TnoOjZmxxf!%Mt12CBoI9&ILG^mx$4| zRu*BEp_Z4r4(a730(7l2mzN07=*53;c?n)sU6#_ZyhH@oS7|i8x{wkQLK!qSetLC@aNWvgb?Ngt)|Y5xwCRxbCFeEl`V#EQpt*)&=lT*+ zs`96-FA<|mnk$pJzC@h1YHlX8gj90_M3NR@ChQzcWV#cTI}JW_eTg)>{dL;R^(CUU zT3BC>^(DHG)BGLlOT=jVA8vh#Sfv`S<&gCyB9(Gu%uK#cnNH2tUYS}nb2>Hqw9~2CXb|#* zdOD@&&~WCwXhj6w0U^2D{P8uI7Gt=qqGD*>6TxR#FgRjUo4pu?!_;Ukm}pEja||-) zy;0ua<|s_K^dv#exiQr+T0!2ds64Nml0+LiEOBZEk@`ysEeUALh+ET})ljx--<(9x4 zG$W0=ML`Z!xHCYGWdB8fm9hOfxGOO9kI^FiSVAB(;S>s0!S6 zAo(2vp)(g$7V~=3mx|DjJYhPC+QO%#vo67n1g=%(tWpx6b`+sW(^GMegB%nYAjZ9L zRM}c{gk~^RJdwbKrVY*H#W02`E@cjejb{(UhaQs2mk0*e_$(WF5Mo33Xl9{QW*8k3 z!88}UQ*fd3nyVSTM>$l!Bf^rvNoC9}c{ayWNz4G-Fe9{MINv=PHVu9+D`ID!#shIi z!^WsP8pB5U6KzGqM#NeYH}Jj9*bOtA=~iNf&?=S}W7wYX!Lfj&j_<`{A3@E@;G$?n zQ;Ow3mNX0NmMv3X8S~YB0qYqDyOIT7b<37}=9FhduV`38WmPx>@q1>TWh>`YV#%}+-^W_cKY`XjE1@JwO5+F% zY(z)P2fd3?09IYK_=MSF-=g^#E8M6-W^lb{z8le~NM`s|^(-Umk{J%c%7z5Z$4A&! zQYET|mzKHEtT)<<#%mcWCXE3Q6Q-m!4AawMdX{<{R{5@6;i+F`Wu)&&&0ppwWf<`q zR(NYsDTElo7S?%J`RZ0QSg|;uW*OvBlB1c@ny%1+%j%Xbn%{umiu2Uh*VRL>iq{!7 z3-$5f8)`!r(wY!Pl#SxV!3&%?b z$_PU?es>Qsfgm<30{A7oC~e+Slx`bZGF-|5#vNc_j@{50jkan=HJc>&<-L-B3a>6k z_gEQ(u!y&V1sIIE^xAMJf;z0=A=fy@Swk6{bTnED2RaJ+QZLMWC2=~9M zs)Ore;G}&k&4|$$w~z4)ei)v0P&wbzrXK@b&UgB0{83{wj}9atZx9sdTls#qr=Y_0 zc+9^gxQzD_C|5KVW)zf>9+Y9+>FTn|c+o?=G(%E9c}PDHh|<)C1)_Lkth5UIsMrcJ zvj!gQf|0dr@#K>#LQOEtGgLZIRY2XcA~}@ot^UR!-Na(a)UyXo;Y5t@)Mm&_p^hMEz8xAr|Oa7cA<4y zs+{LGl$J+U9v4-rNHWMRm=3DeO;#cJ1f9Y|37lRxhiDGR6G>k*!rNt(Pg*!Dk=ae7hL7-m86!+=pcM*3cT|XT zr69Dg3R!DFI8nwPh{~{qsw|u#Dl3Gj+N#~Te6gOgfr}FAE(E}hCPY=7pYUO2pmh4Z zOUNySS9*c56i^51xbXN%x--YP zsZlu{CL)XY7freP`2yc+)QFvv#w?+#D|U`qEwsA4J{7^OAzd|iLSK4c6m3EQrm3Ms zsy99`>d299&GAhn)D%58J=~O5l&kSZP;XPW(7~hoGqpH+f(lmsTs_=p88KY&ZXg&{ zx3V!9EjR`Z&RhQ0)d}i7X-duT9o6A4l#;+CKkd$>rJX%{ww6{bUbt|MYX*LqIp3uo zJ#&>HwxqbM*fj%>IAMSdV$FkxaQ|OoTLyuK3Miqr5N2H#-^A3g)d{KQ9JrZRl={naDjds@ej zcJSQ{Om(97bLeMh9|9nqn#0uRHUt>H{ovaVzL-%4;)svV4hP>nd~0tWgvr4-pUgoz zrsE!*CO`V#tb2<3jdza3AIggkdS8k7NH43ww;Ft}==V7gBwX|d$HO>@kD?0j9om8! z>YKJ8`7rr?3zXrbU&E~epNS(zgTL?Kqc3=04n9=H9)@qHgU<~<_i6Zk?BJvC;M)OA_Cg2o1*fA}W+_VqA)V-e;o zcfv+tznw$urYyRW4VYu+Ekd+IQ5Pbh8Yk6*Y-=P=XMW{d(Ptb1CO@#9a*#iwz-(PA z`5i^8lF7d)T4hLnMA1|_+6DeX(WSMTMaLui}+vsGlpEFM|5BqNx_pgQZ6? zP5qNVHhVJoIZIOu3PJsiuG#uy1Jth+tzJf^{-fwIbm4xYXw@sJKPZ}71>a7zw0g0V z>|4bvk4<*0X!20XUKKroX|hX2lV?gF*^cjYQE0aPI+Wtc7Wn8M4r5o@_U%3sPoqfs z$hL1ss<3QWg|nS=#2*GpL-g%W#W6}abD;RiL*mg*>|xT!9O=suOFyBH0~Puqh)c`qW`&{Y0G(DQ8l#tAB~O7m0v(DJCio9~)a z{f`m9O+R+n_AvHIE zde4pG>Dx`Kk)Ns$<}ILOps6QP{sGYKpldZ90Zrp9Na2G0Zvedw@nt%`6SUFyIDB!B z-XY!Y$nWc*p9Xz{=6@1&H|X;uc|6#q_8Dmn zq2y3h`VlX+{{}gl8i(0)!sc5ut<#0iy=vL~r8NtMPbmgZ_66|(8CFti5rbz#nxpgP zEynXx3w@=9#Yjj$_#u4FF<)q%Xy%>SW+y{>^T)9>Me1>GgQU02kYPmq)P!@(MW9?P zwdf&|e;^NT4dDijQvlQ)OYXR!=ZxdYjLP$78gj+rOUKg23IcG?I-o*u(=i;$%ua32 zsQr0}+rrOsrGkQ5PiJ$g8gN505(YPfNTw2^E!Y-_ZOF_Qx5Ketl5p%GfxnqHAXpVn zW`8mocGA9}GM&Y>bqjEJp>EM4+_~{J%wJIJfqh_mamdqzjx9Nn7)DsOL|qri0MW5( aMpMV*(fCqq-(wpbT56%!%UG+h!v6ul4e$H_ delta 5956 zcmZu#4OG-s7JvU?7={mV_yz@L1R-?6G|@C&m&)nrhMvjP)cinU^sgZx43_yZ(1=1r zLD$_)7Yo-utmA22sSc&vV$y1@m9=ZfT5HM08WIf+lagY8_x=BWpw15Ga^LUWk9XgF z_r00-n*Lxq_<~X&A7ZOETeF1_HfKZ1s-)lyA(#TuBWyw_U}3;vEoN(2!@juST>!fP zb~PM(z-fmx925tN2l2IX_9=0bhNI64 ziU9E^NMj-f+yV*##exz*2vFigw)9p49o%8IHe(1L^?A@y;9PLV91VvK2Jct;XF#uK z{~*o)J_Y@4kf(s=g3bY-1ZVsAL0!P(bX)-52h6j9+$%mtBW$PHdK8=;dA;*?#zml9U}Hk9z)OKAf=6O};Y6bl z9h@6BI04!Myg|o};Ew@wE?NFEs1ewu;~&8f0{eArLLf(h$LQGki-6P7pJ`-tHd49> z3X6u!CYZkDfcdp52tEO(=sU+7+0%DlL!$er&g zScSH`w6wT%j=QXE$}BF}0% zYkp<8Nth6zBL53vqGtk=P@W3@UZZ1yQQl>kE`REg6AOx)7dH?EI|Z0iF`YSFe>S0Cbz?w`<; z#)(%e7`IEz3!qjrJ}ogXgj&zIRbpNc^(n@OB<6)tn;17q%nPI*U|bIjE+$-BUh{>o zqt`WmUX;qmnX2MUx@Vy}|7mRG8?@Xy!W#`MQ9k2h18Pt5YpoPy%t`#IVgF} zIuzELP4ji@I7RVMdWhYUgVVr2gYFKQ8oPj-{eAuYzBU_=c^~6%XitdKvY4CC(uI(@ zN;zeR-fmgS9pSV&bd(ZCdqd+AdT!`;ee0~XPOGmq##QHRarw{Dd!dUg`K-R2eifEr zD!4&UhsK(6Z_vJIeAk3!#;#}Kmp`z*gU9tS{+OD>rUh-}=I?%>xx-_Xy%aWVmPH16 zGc6f5Q5i{(4;y2;@;yB+lYDE&L>g7PBbs3{`c z^usk7-0_G^iwtfTU5nV3+{oT-*Q9E{v!ip&k3-_@rH06)pr={oOpi>_XLKsEdLWLd z$x1vmL}e-$=}r89lP*T-A>`2P=q%+HdN4Xe38tp#94#UHud;UsCuGLZgw)bc(PJzx z@jy!{D<)1UAZJXTmWt1>$ly46G8OOBk(fl~LpmFis|B)x#>P%D-F1ami7H~FP0>cv zN^iu5Dd*`#tR7YYrNw2L7Isr{cr2}pi#CnzmPw&FeZ|_o83-^-4^U?4Abl3MO)Jt2 zS{a`nbcAE=yG)PA$0}xeE?y6I+vT4JTSOxhQk8Y|n*{fGSy=zM^mBy|4BWTQjH<3( za8@~uy6$%QQ4I;MAZfJX65VgR)pXY-W6e#`+#JtS8?h)z)?q;xEsV$vTG%E153Ksw zE^UG%^;K`9w@2!$o^w&E@~YpguezKoCYYG7&Gm2eZsPQClctCjSf}{krs$;EceeBN z)4$do^j&PL`HpwKJSyU$h^e^zv6{T5)z_<(r}0_Q53_yz9qICKqNbz-uj4r!O00WV zr}S0RTI1}K_I)n@lSZGf*W{7*c~IBw(=?$m)>XM$e5%wsVH!7Rs%WN7gQj2mO8?~G)<=Oyd; zO0sS~q)#}zt^r}vwazXVyR3EgfcjkO4fNMlorDvYKUofgLn*n(=<51vAQR{0bk4c_ zPmyAe$#y(HF#GG8lCMom7E~^1@z=DvROg`H2lQKeR>&-L)HcN381SvggNPAT~Tu8sF^!qaX+mWy?XVC9dLsg4=$q)#?*R?N< z`V-C$UXeOax{&$A3~icbP3u%=$Ci-JU*mU8y(Y(Ubzr6#!+$p|A2niH-k@qbRIlNS zaZt4#D++go7qlgFo|5aW_fDqfQKOVLIy1_ic#IwXSyQ_hEwxU&)PHdAZnC8$DutAl zl4SXqI~J2G#jcE@k`#N;=iCu}fqtKojC20g6uVdUe9$2WtlWAk3RjKt+fWC!P9C;H z3+Ne5j@?xCirg4jLvBEdVZUUjk7+8k&QsHzC(CF5(!;V}dO%`edO477@1eZZ81H_q z1p1xrojo`|IDdSW9MH}ZnV7Tx)fc{3{T7CDEw^iZywlGb>v3pM&vZc7Tl#EYcd`A#p~C*i#%nSjt3uDwW6q?G{d91 zSFCXtuOH)BS;n1?)uqMj9cg8xY_(Hsy0tmuCaZ#jWgfBhY!! zRgh_?*%}E-0o?++12hk`3RDVu1oR~6PoTq~w?H3(&V#OkOt{{&fBpq$(>dpcuj2Tj zm~Q3!H2-YmOKL_X?knK!KHOixyMbB0jILzFc{33z-+vB)`1Z~>biVnwRGF>ZJ^;?Q zN*)Gx4|$76l@(W`AHRi*IfQQTS-@z@5BCYNL!`ohJNb^0jgYzDXlJ60i(>Q6F=>`i z7K*h+Ydz!gX~Vsa=1^B(B&BTLOs5~-+`M=5w~BI-8n)U)IP`9K=~n2Ct;xYWqxJwH z*3cJQ?ZF(uEKGC;1#hzla(va4wJkO909prV0a}O9Iz(I1;@P#(b7=AJYL(tYs}Ze6 zy0$HK0>7Tv5O@baZ|EPMKYyM8@jN(_oZ+5ST>t6AZBew;8y4)qm!D%ae#p}XuRS;n zxs~rMQhv}I7RbNAB=QuP9`QH0`%7%&ZdUBT5x|wjx@}3#dUH3Yp*w)S4}RN;OuG`mt1@} z#~S8plZB~cenrB)mJ`mc09m&vF}au=|?qv`T_;D)^k_A0OuvVwiL8F-iOpXEDPAFCb@uHd1mUWi*n z(+~xI_>QzW_Yw&?a65PNd!2?JHYRr~vqch~xNgU>@cZV8_5TA_aT(>-sGJrPnVZ z&lu2;65j&HQ$m&9)U%&%s~X|uL#Ta7og*0s=!!~kP~kD?TECB`Iy zjRTZl^rit1ebKuEc<75>KJYo%pVaK>pH(s5p)Y#dAwPievw@Qa0A|06R0p$zgEx>&5&%?gtfNurPgg@g&aX#wKH>d?RG8~D`)&N`xx$)w-2Y3tc23`Lk@MFMu zdy)gb#`egNOX3i5%!`gjC}n<{{iR- BXM6wv diff --git a/examples/wkb_reader/readwkb.c b/examples/wkb_reader/readwkb.c index be40e62c8..33c07cdbe 100644 --- a/examples/wkb_reader/readwkb.c +++ b/examples/wkb_reader/readwkb.c @@ -157,7 +157,7 @@ int WKB_OID; double *double_val; char *char_val; char *wkb_val; - char *table_name = "tt"; + char *table_name = "t"; char query_str[1000]; /* @@ -172,7 +172,7 @@ int WKB_OID; pgport = "5555"; pgoptions = "user=postgres"; pgtty = NULL; - dbName = "test1"; + dbName = "t2"; @@ -217,7 +217,9 @@ int WKB_OID; * fetch rows from the pg_database, the system catalog of * databases */ - sprintf(query_str, "DECLARE mycursor BINARY CURSOR FOR select interesting, comments, wkb_ndr(the_geom) as wkb from %s",table_name); + sprintf(query_str, "DECLARE mycursor BINARY CURSOR FOR select text(num), asbinary(the_geom,'ndr') as wkb from %s",table_name); + + printf(query_str); printf("\n"); res = PQexec(conn, query_str); @@ -271,7 +273,7 @@ int WKB_OID; double_val = (double *) PQgetvalue(res, row, field); printf("%s: %g\n",field_name,*double_val); } - if (field_type ==1043 )//varchar + if ( (field_type ==1043 ) || (field_type==25) )//varchar { char_val = (char *) PQgetvalue(res, row, field); printf("%s: %s\n",field_name,char_val); @@ -300,4 +302,4 @@ int WKB_OID; PQfinish(conn); return 0; -} \ No newline at end of file +} diff --git a/loader/Makefile b/loader/Makefile new file mode 100644 index 000000000..4d2ee46c8 --- /dev/null +++ b/loader/Makefile @@ -0,0 +1,11 @@ +OBJS = shp2pgsql.o shpopen.o dbfopen.o +PROG = shp2pgsql + +all: $(PROG) + +$(PROG): $(OBJS) + +clean: + @rm -f $(OBJS) $(PROG) + + diff --git a/loader/README b/loader/README new file mode 100644 index 000000000..a1b5a9082 --- /dev/null +++ b/loader/README @@ -0,0 +1,63 @@ +shp2pgsql - Convert Shape file to PostGIS +~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +VERSION: 0.5 (2001/07/20) + +MORE INFORMATION: http://postgis.refractions.net + +INTRODUCTION: +This program takes in ESRI shape files and output formatted text suitable +for uploading to the PostGIS/PostgreSQL spatial database using the psql +terminal monitor. + +This application uses functionality from shapelib 1.2.8 +by Frank Warmerdam to read from ESRI +Shape files. + + +INSTALLATION: + +To build shp2pgsql just run 'make'. +Copy the binary wherever you like. :) + + +USAGE: + +shp2pgsql -d -c -a -dump + +The is the name of the shape file, without any extension +information. For example, 'roads' would be the name of the shapefile +comprising the 'roads.shp', 'roads.shx', and 'roads.dbf' files. + +The is the name of the database table you want the data stored +in in the database. Within that table, the geometry will be placed in +the 'geo_value' column by default. + +The options are as follows: + + -d Delete mode. Delete the database table before uploading the + data into a new empty database table in 'insert' format. + -c Create mode. This is the default mode. Create a new table and + upload the data into that table in 'insert' format. + -a Append mode. Do not delete the target table or try to create + a new table, simple insert the data into the existing table. + (A table will have to exist for this to work, it is usually + used after a create mode as been run once.) + -dump Dump mode. Create a new table and upload the data into that + table in 'dump' format. Dump format is used by PostgreSQL for + large data dumps and uploads. Use this mode if your upload + dataset is very large. + + +EXAMPLES: + +Loading directly: + + shp2pgsql roads1 -c | psql -d roadsdatabase + shp2pgsql roads2 -a | psql -d roadsdatabase + +Saving to an intermiate file: + + shp2pgsql roads1 > roads.sql + psql -d roadsdatabase -f roads.sql + diff --git a/loader/dbfopen.c b/loader/dbfopen.c new file mode 100644 index 000000000..790350650 --- /dev/null +++ b/loader/dbfopen.c @@ -0,0 +1,1074 @@ +/****************************************************************************** + * $Id$ + * + * Project: Shapelib + * Purpose: Implementation of .dbf access API documented in dbf_api.html. + * Author: Frank Warmerdam, warmerda@home.com + * + ****************************************************************************** + * Copyright (c) 1999, Frank Warmerdam + * + * This software is available under the following "MIT Style" license, + * or at the option of the licensee under the LGPL (see LICENSE.LGPL). This + * option is discussed in more detail in shapelib.html. + * + * -- + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************** + * + * $Log$ + * Revision 1.1 2001/07/18 21:42:12 pramsey + * Initial add of the data loader code. + * + * Revision 1.22 1999/12/15 13:47:24 warmerda + * Added stdlib.h to ensure that atof() is prototyped. + * + * Revision 1.21 1999/12/13 17:25:46 warmerda + * Added support for upper case .DBF extention. + * + * Revision 1.20 1999/11/30 16:32:11 warmerda + * Use atof() instead of sscanf(). + * + * Revision 1.19 1999/11/05 14:12:04 warmerda + * updated license terms + * + * Revision 1.18 1999/07/27 00:53:28 warmerda + * ensure that whole old field value clear on write of string + * + * Revision 1.1 1999/07/05 18:58:07 warmerda + * New + * + * Revision 1.17 1999/06/11 19:14:12 warmerda + * Fixed some memory leaks. + * + * Revision 1.16 1999/06/11 19:04:11 warmerda + * Remoted some unused variables. + * + * Revision 1.15 1999/05/11 03:19:28 warmerda + * added new Tuple api, and improved extension handling - add from candrsn + * + * Revision 1.14 1999/05/04 15:01:48 warmerda + * Added 'F' support. + * + * Revision 1.13 1999/03/23 17:38:59 warmerda + * DBFAddField() now actually does return the new field number, or -1 if + * it fails. + * + * Revision 1.12 1999/03/06 02:54:46 warmerda + * Added logic to convert shapefile name to dbf filename in DBFOpen() + * for convenience. + * + * Revision 1.11 1998/12/31 15:30:34 warmerda + * Improved the interchangability of numeric and string attributes. Add + * white space trimming option for attributes. + * + * Revision 1.10 1998/12/03 16:36:44 warmerda + * Use r+b instead of rb+ for binary access. + * + * Revision 1.9 1998/12/03 15:34:23 warmerda + * Updated copyright message. + * + * Revision 1.8 1997/12/04 15:40:15 warmerda + * Added newline character after field definitions. + * + * Revision 1.7 1997/03/06 14:02:10 warmerda + * Ensure bUpdated is initialized. + * + * Revision 1.6 1996/02/12 04:54:41 warmerda + * Ensure that DBFWriteAttribute() returns TRUE if it succeeds. + * + * Revision 1.5 1995/10/21 03:15:12 warmerda + * Changed to use binary file access, and ensure that the + * field name field is zero filled, and limited to 10 chars. + * + * Revision 1.4 1995/08/24 18:10:42 warmerda + * Added use of SfRealloc() to avoid pre-ANSI realloc() functions such + * as on the Sun. + * + * Revision 1.3 1995/08/04 03:15:16 warmerda + * Fixed up header. + * + * Revision 1.2 1995/08/04 03:14:43 warmerda + * Added header. + */ + +static char rcsid[] = + "$Id$"; + +#include "shapefil.h" + +#include +#include + +typedef unsigned char uchar; + +#ifndef FALSE +# define FALSE 0 +# define TRUE 1 +#endif + +static int nStringFieldLen = 0; +static char * pszStringField = NULL; + +/************************************************************************/ +/* SfRealloc() */ +/* */ +/* A realloc cover function that will access a NULL pointer as */ +/* a valid input. */ +/************************************************************************/ + +static void * SfRealloc( void * pMem, int nNewSize ) + +{ + if( pMem == NULL ) + return( (void *) malloc(nNewSize) ); + else + return( (void *) realloc(pMem,nNewSize) ); +} + +/************************************************************************/ +/* DBFWriteHeader() */ +/* */ +/* This is called to write out the file header, and field */ +/* descriptions before writing any actual data records. This */ +/* also computes all the DBFDataSet field offset/size/decimals */ +/* and so forth values. */ +/************************************************************************/ + +static void DBFWriteHeader(DBFHandle psDBF) + +{ + uchar abyHeader[XBASE_FLDHDR_SZ]; + int i; + + if( !psDBF->bNoHeader ) + return; + + psDBF->bNoHeader = FALSE; + +/* -------------------------------------------------------------------- */ +/* Initialize the file header information. */ +/* -------------------------------------------------------------------- */ + for( i = 0; i < XBASE_FLDHDR_SZ; i++ ) + abyHeader[i] = 0; + + abyHeader[0] = 0x03; /* memo field? - just copying */ + + /* date updated on close, record count preset at zero */ + + abyHeader[8] = psDBF->nHeaderLength % 256; + abyHeader[9] = psDBF->nHeaderLength / 256; + + abyHeader[10] = psDBF->nRecordLength % 256; + abyHeader[11] = psDBF->nRecordLength / 256; + +/* -------------------------------------------------------------------- */ +/* Write the initial 32 byte file header, and all the field */ +/* descriptions. */ +/* -------------------------------------------------------------------- */ + fseek( psDBF->fp, 0, 0 ); + fwrite( abyHeader, XBASE_FLDHDR_SZ, 1, psDBF->fp ); + fwrite( psDBF->pszHeader, XBASE_FLDHDR_SZ, psDBF->nFields, psDBF->fp ); + +/* -------------------------------------------------------------------- */ +/* Write out the newline character if there is room for it. */ +/* -------------------------------------------------------------------- */ + if( psDBF->nHeaderLength > 32*psDBF->nFields + 32 ) + { + char cNewline; + + cNewline = 0x0d; + fwrite( &cNewline, 1, 1, psDBF->fp ); + } +} + +/************************************************************************/ +/* DBFFlushRecord() */ +/* */ +/* Write out the current record if there is one. */ +/************************************************************************/ + +static void DBFFlushRecord( DBFHandle psDBF ) + +{ + int nRecordOffset; + + if( psDBF->bCurrentRecordModified && psDBF->nCurrentRecord > -1 ) + { + psDBF->bCurrentRecordModified = FALSE; + + nRecordOffset = psDBF->nRecordLength * psDBF->nCurrentRecord + + psDBF->nHeaderLength; + + fseek( psDBF->fp, nRecordOffset, 0 ); + fwrite( psDBF->pszCurrentRecord, psDBF->nRecordLength, 1, psDBF->fp ); + } +} + +/************************************************************************/ +/* DBFOpen() */ +/* */ +/* Open a .dbf file. */ +/************************************************************************/ + +DBFHandle DBFOpen( const char * pszFilename, const char * pszAccess ) + +{ + DBFHandle psDBF; + uchar *pabyBuf; + int nFields, nRecords, nHeadLen, nRecLen, iField, i; + char *pszBasename, *pszFullname; + +/* -------------------------------------------------------------------- */ +/* We only allow the access strings "rb" and "r+". */ +/* -------------------------------------------------------------------- */ + if( strcmp(pszAccess,"r") != 0 && strcmp(pszAccess,"r+") != 0 + && strcmp(pszAccess,"rb") != 0 && strcmp(pszAccess,"rb+") != 0 + && strcmp(pszAccess,"r+b") != 0 ) + return( NULL ); + +/* -------------------------------------------------------------------- */ +/* Compute the base (layer) name. If there is any extension */ +/* on the passed in filename we will strip it off. */ +/* -------------------------------------------------------------------- */ + pszBasename = (char *) malloc(strlen(pszFilename)+5); + strcpy( pszBasename, pszFilename ); + for( i = strlen(pszBasename)-1; + i > 0 && pszBasename[i] != '.' && pszBasename[i] != '/' + && pszBasename[i] != '\\'; + i-- ) {} + + if( pszBasename[i] == '.' ) + pszBasename[i] = '\0'; + + pszFullname = (char *) malloc(strlen(pszBasename) + 5); + sprintf( pszFullname, "%s.dbf", pszBasename ); + + psDBF = (DBFHandle) calloc( 1, sizeof(DBFInfo) ); + psDBF->fp = fopen( pszFullname, pszAccess ); + + if( psDBF->fp == NULL ) + { + sprintf( pszFullname, "%s.DBF", pszBasename ); + psDBF->fp = fopen(pszFullname, pszAccess ); + } + + free( pszBasename ); + free( pszFullname ); + + if( psDBF->fp == NULL ) + { + free( psDBF ); + return( NULL ); + } + + psDBF->bNoHeader = FALSE; + psDBF->nCurrentRecord = -1; + psDBF->bCurrentRecordModified = FALSE; + +/* -------------------------------------------------------------------- */ +/* Read Table Header info */ +/* -------------------------------------------------------------------- */ + pabyBuf = (uchar *) malloc(500); + fread( pabyBuf, 32, 1, psDBF->fp ); + + psDBF->nRecords = nRecords = + pabyBuf[4] + pabyBuf[5]*256 + pabyBuf[6]*256*256 + pabyBuf[7]*256*256*256; + + psDBF->nHeaderLength = nHeadLen = pabyBuf[8] + pabyBuf[9]*256; + psDBF->nRecordLength = nRecLen = pabyBuf[10] + pabyBuf[11]*256; + + psDBF->nFields = nFields = (nHeadLen - 32) / 32; + + psDBF->pszCurrentRecord = (char *) malloc(nRecLen); + +/* -------------------------------------------------------------------- */ +/* Read in Field Definitions */ +/* -------------------------------------------------------------------- */ + + pabyBuf = (uchar *) SfRealloc(pabyBuf,nHeadLen); + psDBF->pszHeader = (char *) pabyBuf; + + fseek( psDBF->fp, 32, 0 ); + fread( pabyBuf, nHeadLen, 1, psDBF->fp ); + + psDBF->panFieldOffset = (int *) malloc(sizeof(int) * nFields); + psDBF->panFieldSize = (int *) malloc(sizeof(int) * nFields); + psDBF->panFieldDecimals = (int *) malloc(sizeof(int) * nFields); + psDBF->pachFieldType = (char *) malloc(sizeof(char) * nFields); + + for( iField = 0; iField < nFields; iField++ ) + { + uchar *pabyFInfo; + + pabyFInfo = pabyBuf+iField*32; + + if( pabyFInfo[11] == 'N' || pabyFInfo[11] == 'F' ) + { + psDBF->panFieldSize[iField] = pabyFInfo[16]; + psDBF->panFieldDecimals[iField] = pabyFInfo[17]; + } + else + { + psDBF->panFieldSize[iField] = pabyFInfo[16] + pabyFInfo[17]*256; + psDBF->panFieldDecimals[iField] = 0; + } + + psDBF->pachFieldType[iField] = (char) pabyFInfo[11]; + if( iField == 0 ) + psDBF->panFieldOffset[iField] = 1; + else + psDBF->panFieldOffset[iField] = + psDBF->panFieldOffset[iField-1] + psDBF->panFieldSize[iField-1]; + } + + return( psDBF ); +} + +/************************************************************************/ +/* DBFClose() */ +/************************************************************************/ + +void DBFClose(DBFHandle psDBF) +{ +/* -------------------------------------------------------------------- */ +/* Write out header if not already written. */ +/* -------------------------------------------------------------------- */ + if( psDBF->bNoHeader ) + DBFWriteHeader( psDBF ); + + DBFFlushRecord( psDBF ); + +/* -------------------------------------------------------------------- */ +/* Update last access date, and number of records if we have */ +/* write access. */ +/* -------------------------------------------------------------------- */ + if( psDBF->bUpdated ) + { + uchar abyFileHeader[32]; + + fseek( psDBF->fp, 0, 0 ); + fread( abyFileHeader, 32, 1, psDBF->fp ); + + abyFileHeader[1] = 95; /* YY */ + abyFileHeader[2] = 7; /* MM */ + abyFileHeader[3] = 26; /* DD */ + + abyFileHeader[4] = psDBF->nRecords % 256; + abyFileHeader[5] = (psDBF->nRecords/256) % 256; + abyFileHeader[6] = (psDBF->nRecords/(256*256)) % 256; + abyFileHeader[7] = (psDBF->nRecords/(256*256*256)) % 256; + + fseek( psDBF->fp, 0, 0 ); + fwrite( abyFileHeader, 32, 1, psDBF->fp ); + } + +/* -------------------------------------------------------------------- */ +/* Close, and free resources. */ +/* -------------------------------------------------------------------- */ + fclose( psDBF->fp ); + + if( psDBF->panFieldOffset != NULL ) + { + free( psDBF->panFieldOffset ); + free( psDBF->panFieldSize ); + free( psDBF->panFieldDecimals ); + free( psDBF->pachFieldType ); + } + + free( psDBF->pszHeader ); + free( psDBF->pszCurrentRecord ); + + free( psDBF ); + + if( pszStringField != NULL ) + { + free( pszStringField ); + pszStringField = NULL; + nStringFieldLen = 0; + } +} + +/************************************************************************/ +/* DBFCreate() */ +/* */ +/* Create a new .dbf file. */ +/************************************************************************/ + +DBFHandle DBFCreate( const char * pszFilename ) + +{ + DBFHandle psDBF; + FILE *fp; + char *pszFullname, *pszBasename; + int i; + +/* -------------------------------------------------------------------- */ +/* Compute the base (layer) name. If there is any extension */ +/* on the passed in filename we will strip it off. */ +/* -------------------------------------------------------------------- */ + pszBasename = (char *) malloc(strlen(pszFilename)+5); + strcpy( pszBasename, pszFilename ); + for( i = strlen(pszBasename)-1; + i > 0 && pszBasename[i] != '.' && pszBasename[i] != '/' + && pszBasename[i] != '\\'; + i-- ) {} + + if( pszBasename[i] == '.' ) + pszBasename[i] = '\0'; + + pszFullname = (char *) malloc(strlen(pszBasename) + 5); + sprintf( pszFullname, "%s.dbf", pszBasename ); + free( pszBasename ); + +/* -------------------------------------------------------------------- */ +/* Create the file. */ +/* -------------------------------------------------------------------- */ + fp = fopen( pszFullname, "wb" ); + if( fp == NULL ) + return( NULL ); + + fputc( 0, fp ); + fclose( fp ); + + fp = fopen( pszFullname, "rb+" ); + if( fp == NULL ) + return( NULL ); + + free( pszFullname ); + +/* -------------------------------------------------------------------- */ +/* Create the info structure. */ +/* -------------------------------------------------------------------- */ + psDBF = (DBFHandle) malloc(sizeof(DBFInfo)); + + psDBF->fp = fp; + psDBF->nRecords = 0; + psDBF->nFields = 0; + psDBF->nRecordLength = 1; + psDBF->nHeaderLength = 33; + + psDBF->panFieldOffset = NULL; + psDBF->panFieldSize = NULL; + psDBF->panFieldDecimals = NULL; + psDBF->pachFieldType = NULL; + psDBF->pszHeader = NULL; + + psDBF->nCurrentRecord = -1; + psDBF->bCurrentRecordModified = FALSE; + psDBF->pszCurrentRecord = NULL; + + psDBF->bNoHeader = TRUE; + + return( psDBF ); +} + +/************************************************************************/ +/* DBFAddField() */ +/* */ +/* Add a field to a newly created .dbf file before any records */ +/* are written. */ +/************************************************************************/ + +int DBFAddField(DBFHandle psDBF, const char * pszFieldName, + DBFFieldType eType, int nWidth, int nDecimals ) + +{ + char *pszFInfo; + int i; + +/* -------------------------------------------------------------------- */ +/* Do some checking to ensure we can add records to this file. */ +/* -------------------------------------------------------------------- */ + if( psDBF->nRecords > 0 ) + return( -1 ); + + if( !psDBF->bNoHeader ) + return( -1 ); + + if( eType != FTDouble && nDecimals != 0 ) + return( -1 ); + +/* -------------------------------------------------------------------- */ +/* SfRealloc all the arrays larger to hold the additional field */ +/* information. */ +/* -------------------------------------------------------------------- */ + psDBF->nFields++; + + psDBF->panFieldOffset = (int *) + SfRealloc( psDBF->panFieldOffset, sizeof(int) * psDBF->nFields ); + + psDBF->panFieldSize = (int *) + SfRealloc( psDBF->panFieldSize, sizeof(int) * psDBF->nFields ); + + psDBF->panFieldDecimals = (int *) + SfRealloc( psDBF->panFieldDecimals, sizeof(int) * psDBF->nFields ); + + psDBF->pachFieldType = (char *) + SfRealloc( psDBF->pachFieldType, sizeof(char) * psDBF->nFields ); + +/* -------------------------------------------------------------------- */ +/* Assign the new field information fields. */ +/* -------------------------------------------------------------------- */ + psDBF->panFieldOffset[psDBF->nFields-1] = psDBF->nRecordLength; + psDBF->nRecordLength += nWidth; + psDBF->panFieldSize[psDBF->nFields-1] = nWidth; + psDBF->panFieldDecimals[psDBF->nFields-1] = nDecimals; + + if( eType == FTString ) + psDBF->pachFieldType[psDBF->nFields-1] = 'C'; + else + psDBF->pachFieldType[psDBF->nFields-1] = 'N'; + +/* -------------------------------------------------------------------- */ +/* Extend the required header information. */ +/* -------------------------------------------------------------------- */ + psDBF->nHeaderLength += 32; + psDBF->bUpdated = FALSE; + + psDBF->pszHeader = (char *) SfRealloc(psDBF->pszHeader,psDBF->nFields*32); + + pszFInfo = psDBF->pszHeader + 32 * (psDBF->nFields-1); + + for( i = 0; i < 32; i++ ) + pszFInfo[i] = '\0'; + + if( strlen(pszFieldName) < 10 ) + strncpy( pszFInfo, pszFieldName, strlen(pszFieldName)); + else + strncpy( pszFInfo, pszFieldName, 10); + + pszFInfo[11] = psDBF->pachFieldType[psDBF->nFields-1]; + + if( eType == FTString ) + { + pszFInfo[16] = nWidth % 256; + pszFInfo[17] = nWidth / 256; + } + else + { + pszFInfo[16] = nWidth; + pszFInfo[17] = nDecimals; + } + +/* -------------------------------------------------------------------- */ +/* Make the current record buffer appropriately larger. */ +/* -------------------------------------------------------------------- */ + psDBF->pszCurrentRecord = (char *) SfRealloc(psDBF->pszCurrentRecord, + psDBF->nRecordLength); + + return( psDBF->nFields-1 ); +} + +/************************************************************************/ +/* DBFReadAttribute() */ +/* */ +/* Read one of the attribute fields of a record. */ +/************************************************************************/ + +static void *DBFReadAttribute(DBFHandle psDBF, int hEntity, int iField, + char chReqType ) + +{ + int nRecordOffset; + uchar *pabyRec; + void *pReturnField = NULL; + + static double dDoubleField; + +/* -------------------------------------------------------------------- */ +/* Have we read the record? */ +/* -------------------------------------------------------------------- */ + if( hEntity < 0 || hEntity >= psDBF->nRecords ) + return( NULL ); + + if( psDBF->nCurrentRecord != hEntity ) + { + DBFFlushRecord( psDBF ); + + nRecordOffset = psDBF->nRecordLength * hEntity + psDBF->nHeaderLength; + + fseek( psDBF->fp, nRecordOffset, 0 ); + fread( psDBF->pszCurrentRecord, psDBF->nRecordLength, 1, psDBF->fp ); + + psDBF->nCurrentRecord = hEntity; + } + + pabyRec = (uchar *) psDBF->pszCurrentRecord; + +/* -------------------------------------------------------------------- */ +/* Ensure our field buffer is large enough to hold this buffer. */ +/* -------------------------------------------------------------------- */ + if( psDBF->panFieldSize[iField]+1 > nStringFieldLen ) + { + nStringFieldLen = psDBF->panFieldSize[iField]*2 + 10; + pszStringField = (char *) SfRealloc(pszStringField,nStringFieldLen); + } + +/* -------------------------------------------------------------------- */ +/* Extract the requested field. */ +/* -------------------------------------------------------------------- */ + strncpy( pszStringField, pabyRec+psDBF->panFieldOffset[iField], + psDBF->panFieldSize[iField] ); + pszStringField[psDBF->panFieldSize[iField]] = '\0'; + + pReturnField = pszStringField; + +/* -------------------------------------------------------------------- */ +/* Decode the field. */ +/* -------------------------------------------------------------------- */ + if( chReqType == 'N' ) + { + dDoubleField = atof(pszStringField); + + pReturnField = &dDoubleField; + } + +/* -------------------------------------------------------------------- */ +/* Should we trim white space off the string attribute value? */ +/* -------------------------------------------------------------------- */ +#ifdef TRIM_DBF_WHITESPACE + else + { + char *pchSrc, *pchDst; + + pchDst = pchSrc = pszStringField; + while( *pchSrc == ' ' ) + pchSrc++; + + while( *pchSrc != '\0' ) + *(pchDst++) = *(pchSrc++); + *pchDst = '\0'; + + while( *(--pchDst) == ' ' && pchDst != pszStringField ) + *pchDst = '\0'; + + } +#endif + + return( pReturnField ); +} + +/************************************************************************/ +/* DBFReadIntAttribute() */ +/* */ +/* Read an integer attribute. */ +/************************************************************************/ + +int DBFReadIntegerAttribute( DBFHandle psDBF, int iRecord, int iField ) + +{ + double *pdValue; + + pdValue = (double *) DBFReadAttribute( psDBF, iRecord, iField, 'N' ); + + return( (int) *pdValue ); +} + +/************************************************************************/ +/* DBFReadDoubleAttribute() */ +/* */ +/* Read a double attribute. */ +/************************************************************************/ + +double DBFReadDoubleAttribute( DBFHandle psDBF, int iRecord, int iField ) + +{ + double *pdValue; + + pdValue = (double *) DBFReadAttribute( psDBF, iRecord, iField, 'N' ); + + return( *pdValue ); +} + +/************************************************************************/ +/* DBFReadStringAttribute() */ +/* */ +/* Read a string attribute. */ +/************************************************************************/ + +const char *DBFReadStringAttribute( DBFHandle psDBF, int iRecord, int iField ) + +{ + return( (const char *) DBFReadAttribute( psDBF, iRecord, iField, 'C' ) ); +} + +/************************************************************************/ +/* DBFGetFieldCount() */ +/* */ +/* Return the number of fields in this table. */ +/************************************************************************/ + +int DBFGetFieldCount( DBFHandle psDBF ) + +{ + return( psDBF->nFields ); +} + +/************************************************************************/ +/* DBFGetRecordCount() */ +/* */ +/* Return the number of records in this table. */ +/************************************************************************/ + +int DBFGetRecordCount( DBFHandle psDBF ) + +{ + return( psDBF->nRecords ); +} + +/************************************************************************/ +/* DBFGetFieldInfo() */ +/* */ +/* Return any requested information about the field. */ +/************************************************************************/ + +DBFFieldType DBFGetFieldInfo( DBFHandle psDBF, int iField, char * pszFieldName, + int * pnWidth, int * pnDecimals ) + +{ + if( iField < 0 || iField >= psDBF->nFields ) + return( FTInvalid ); + + if( pnWidth != NULL ) + *pnWidth = psDBF->panFieldSize[iField]; + + if( pnDecimals != NULL ) + *pnDecimals = psDBF->panFieldDecimals[iField]; + + if( pszFieldName != NULL ) + { + int i; + + strncpy( pszFieldName, (char *) psDBF->pszHeader+iField*32, 11 ); + pszFieldName[11] = '\0'; + for( i = 10; i > 0 && pszFieldName[i] == ' '; i-- ) + pszFieldName[i] = '\0'; + } + + if( psDBF->pachFieldType[iField] == 'N' + || psDBF->pachFieldType[iField] == 'F' + || psDBF->pachFieldType[iField] == 'D' ) + { + if( psDBF->panFieldDecimals[iField] > 0 ) + return( FTDouble ); + else + return( FTInteger ); + } + else + { + return( FTString ); + } +} + +/************************************************************************/ +/* DBFWriteAttribute() */ +/* */ +/* Write an attribute record to the file. */ +/************************************************************************/ + +static int DBFWriteAttribute(DBFHandle psDBF, int hEntity, int iField, + void * pValue ) + +{ + int nRecordOffset, i, j; + uchar *pabyRec; + char szSField[40], szFormat[12]; + +/* -------------------------------------------------------------------- */ +/* Is this a valid record? */ +/* -------------------------------------------------------------------- */ + if( hEntity < 0 || hEntity > psDBF->nRecords ) + return( FALSE ); + + if( psDBF->bNoHeader ) + DBFWriteHeader(psDBF); + +/* -------------------------------------------------------------------- */ +/* Is this a brand new record? */ +/* -------------------------------------------------------------------- */ + if( hEntity == psDBF->nRecords ) + { + DBFFlushRecord( psDBF ); + + psDBF->nRecords++; + for( i = 0; i < psDBF->nRecordLength; i++ ) + psDBF->pszCurrentRecord[i] = ' '; + + psDBF->nCurrentRecord = hEntity; + + } + +/* -------------------------------------------------------------------- */ +/* Is this an existing record, but different than the last one */ +/* we accessed? */ +/* -------------------------------------------------------------------- */ + if( psDBF->nCurrentRecord != hEntity ) + { + DBFFlushRecord( psDBF ); + + nRecordOffset = psDBF->nRecordLength * hEntity + psDBF->nHeaderLength; + + fseek( psDBF->fp, nRecordOffset, 0 ); + fread( psDBF->pszCurrentRecord, psDBF->nRecordLength, 1, psDBF->fp ); + + psDBF->nCurrentRecord = hEntity; + } + + pabyRec = (uchar *) psDBF->pszCurrentRecord; + +/* -------------------------------------------------------------------- */ +/* Assign all the record fields. */ +/* -------------------------------------------------------------------- */ + switch( psDBF->pachFieldType[iField] ) + { + case 'D': + case 'N': + case 'F': + if( psDBF->panFieldDecimals[iField] == 0 ) + { + sprintf( szFormat, "%%%dd", psDBF->panFieldSize[iField] ); + sprintf(szSField, szFormat, (int) *((double *) pValue) ); + if( strlen(szSField) > psDBF->panFieldSize[iField] ) + szSField[psDBF->panFieldSize[iField]] = '\0'; + strncpy((char *) (pabyRec+psDBF->panFieldOffset[iField]), + szSField, strlen(szSField) ); + } + else + { + sprintf( szFormat, "%%%d.%df", + psDBF->panFieldSize[iField], + psDBF->panFieldDecimals[iField] ); + sprintf(szSField, szFormat, *((double *) pValue) ); + if( strlen(szSField) > psDBF->panFieldSize[iField] ) + szSField[psDBF->panFieldSize[iField]] = '\0'; + strncpy((char *) (pabyRec+psDBF->panFieldOffset[iField]), + szSField, strlen(szSField) ); + } + + break; + + default: + if( strlen((char *) pValue) > psDBF->panFieldSize[iField] ) + j = psDBF->panFieldSize[iField]; + else + { + memset( pabyRec+psDBF->panFieldOffset[iField], ' ', + psDBF->panFieldSize[iField] ); + j = strlen((char *) pValue); + } + + strncpy((char *) (pabyRec+psDBF->panFieldOffset[iField]), + (char *) pValue, j ); + break; + } + + psDBF->bCurrentRecordModified = TRUE; + psDBF->bUpdated = TRUE; + + return( TRUE ); +} + +/************************************************************************/ +/* DBFWriteDoubleAttribute() */ +/* */ +/* Write a double attribute. */ +/************************************************************************/ + +int DBFWriteDoubleAttribute( DBFHandle psDBF, int iRecord, int iField, + double dValue ) + +{ + return( DBFWriteAttribute( psDBF, iRecord, iField, (void *) &dValue ) ); +} + +/************************************************************************/ +/* DBFWriteIntegerAttribute() */ +/* */ +/* Write a integer attribute. */ +/************************************************************************/ + +int DBFWriteIntegerAttribute( DBFHandle psDBF, int iRecord, int iField, + int nValue ) + +{ + double dValue = nValue; + + return( DBFWriteAttribute( psDBF, iRecord, iField, (void *) &dValue ) ); +} + +/************************************************************************/ +/* DBFWriteStringAttribute() */ +/* */ +/* Write a string attribute. */ +/************************************************************************/ + +int DBFWriteStringAttribute( DBFHandle psDBF, int iRecord, int iField, + const char * pszValue ) + +{ + return( DBFWriteAttribute( psDBF, iRecord, iField, (void *) pszValue ) ); +} + +/************************************************************************/ +/* DBFWriteTuple() */ +/* */ +/* Write an attribute record to the file. */ +/************************************************************************/ + +int DBFWriteTuple(DBFHandle psDBF, int hEntity, void * pRawTuple ) + +{ + int nRecordOffset, i; + uchar *pabyRec; + +/* -------------------------------------------------------------------- */ +/* Is this a valid record? */ +/* -------------------------------------------------------------------- */ + if( hEntity < 0 || hEntity > psDBF->nRecords ) + return( FALSE ); + + if( psDBF->bNoHeader ) + DBFWriteHeader(psDBF); + +/* -------------------------------------------------------------------- */ +/* Is this a brand new record? */ +/* -------------------------------------------------------------------- */ + if( hEntity == psDBF->nRecords ) + { + DBFFlushRecord( psDBF ); + + psDBF->nRecords++; + for( i = 0; i < psDBF->nRecordLength; i++ ) + psDBF->pszCurrentRecord[i] = ' '; + + psDBF->nCurrentRecord = hEntity; + } + +/* -------------------------------------------------------------------- */ +/* Is this an existing record, but different than the last one */ +/* we accessed? */ +/* -------------------------------------------------------------------- */ + if( psDBF->nCurrentRecord != hEntity ) + { + DBFFlushRecord( psDBF ); + + nRecordOffset = psDBF->nRecordLength * hEntity + psDBF->nHeaderLength; + + fseek( psDBF->fp, nRecordOffset, 0 ); + fread( psDBF->pszCurrentRecord, psDBF->nRecordLength, 1, psDBF->fp ); + + psDBF->nCurrentRecord = hEntity; + } + + pabyRec = (uchar *) psDBF->pszCurrentRecord; + + memcpy ( pabyRec, pRawTuple, psDBF->nRecordLength ); + + psDBF->bCurrentRecordModified = TRUE; + psDBF->bUpdated = TRUE; + + return( TRUE ); +} + +/************************************************************************/ +/* DBFReadTuple() */ +/* */ +/* Read one of the attribute fields of a record. */ +/************************************************************************/ + +const char *DBFReadTuple(DBFHandle psDBF, int hEntity ) + +{ + int nRecordOffset; + uchar *pabyRec; + static char *pReturnTuple = NULL; + + static int nTupleLen = 0; + +/* -------------------------------------------------------------------- */ +/* Have we read the record? */ +/* -------------------------------------------------------------------- */ + if( hEntity < 0 || hEntity >= psDBF->nRecords ) + return( NULL ); + + if( psDBF->nCurrentRecord != hEntity ) + { + DBFFlushRecord( psDBF ); + + nRecordOffset = psDBF->nRecordLength * hEntity + psDBF->nHeaderLength; + + fseek( psDBF->fp, nRecordOffset, 0 ); + fread( psDBF->pszCurrentRecord, psDBF->nRecordLength, 1, psDBF->fp ); + + psDBF->nCurrentRecord = hEntity; + } + + pabyRec = (uchar *) psDBF->pszCurrentRecord; + + if ( nTupleLen < psDBF->nRecordLength) { + nTupleLen = psDBF->nRecordLength; + pReturnTuple = (char *) SfRealloc(pReturnTuple, psDBF->nRecordLength); + } + + memcpy ( pReturnTuple, pabyRec, psDBF->nRecordLength ); + + return( pReturnTuple ); +} + +/************************************************************************/ +/* DBFCloneEmpty() */ +/* */ +/* Read one of the attribute fields of a record. */ +/************************************************************************/ + +DBFHandle DBFCloneEmpty(DBFHandle psDBF, const char * pszFilename ) +{ + DBFHandle newDBF; + + newDBF = DBFCreate ( pszFilename ); + if ( newDBF == NULL ) return ( NULL ); + + newDBF->pszHeader = (void *) malloc ( 32 * psDBF->nFields ); + memcpy ( newDBF->pszHeader, psDBF->pszHeader, 32 * psDBF->nFields ); + + newDBF->nFields = psDBF->nFields; + newDBF->nRecordLength = psDBF->nRecordLength; + newDBF->nHeaderLength = psDBF->nHeaderLength; + + newDBF->panFieldOffset = (void *) malloc ( sizeof(int) * psDBF->nFields ); + memcpy ( newDBF->panFieldOffset, psDBF->panFieldOffset, sizeof(int) * psDBF->nFields ); + newDBF->panFieldSize = (void *) malloc ( sizeof(int) * psDBF->nFields ); + memcpy ( newDBF->panFieldSize, psDBF->panFieldSize, sizeof(int) * psDBF->nFields ); + newDBF->panFieldDecimals = (void *) malloc ( sizeof(int) * psDBF->nFields ); + memcpy ( newDBF->panFieldDecimals, psDBF->panFieldDecimals, sizeof(int) * psDBF->nFields ); + newDBF->pachFieldType = (void *) malloc ( sizeof(int) * psDBF->nFields ); + memcpy ( newDBF->pachFieldType, psDBF->pachFieldType, sizeof(int) * psDBF->nFields ); + + newDBF->bNoHeader = TRUE; + newDBF->bUpdated = TRUE; + + DBFWriteHeader ( newDBF ); + DBFClose ( newDBF ); + + newDBF = DBFOpen ( pszFilename, "rb+" ); + + return ( newDBF ); +} diff --git a/loader/shapefil.h b/loader/shapefil.h new file mode 100644 index 000000000..5b2d781a6 --- /dev/null +++ b/loader/shapefil.h @@ -0,0 +1,351 @@ +#ifndef _SHAPEFILE_H_INCLUDED +#define _SHAPEFILE_H_INCLUDED + +/****************************************************************************** + * $Id$ + * + * Project: Shapelib + * Purpose: Primary include file for Shapelib. + * Author: Frank Warmerdam, warmerda@home.com + * + ****************************************************************************** + * Copyright (c) 1999, Frank Warmerdam + * + * This software is available under the following "MIT Style" license, + * or at the option of the licensee under the LGPL (see LICENSE.LGPL). This + * option is discussed in more detail in shapelib.html. + * + * -- + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************** + * + * $Log$ + * Revision 1.1 2001/07/18 21:42:12 pramsey + * Initial add of the data loader code. + * + * Revision 1.15 2000/02/16 16:03:51 warmerda + * added null shape support + * + * Revision 1.14 1999/11/05 14:12:05 warmerda + * updated license terms + * + * Revision 1.13 1999/06/02 18:24:21 warmerda + * added trimming code + * + * Revision 1.12 1999/06/02 17:56:12 warmerda + * added quad'' subnode support for trees + * + * Revision 1.11 1999/05/18 19:11:11 warmerda + * Added example searching capability + * + * Revision 1.10 1999/05/18 17:49:38 warmerda + * added initial quadtree support + * + * Revision 1.9 1999/05/11 03:19:28 warmerda + * added new Tuple api, and improved extension handling - add from candrsn + * + * Revision 1.8 1999/03/23 17:22:27 warmerda + * Added extern "C" protection for C++ users of shapefil.h. + * + * Revision 1.7 1998/12/31 15:31:07 warmerda + * Added the TRIM_DBF_WHITESPACE and DISABLE_MULTIPATCH_MEASURE options. + * + * Revision 1.6 1998/12/03 15:48:15 warmerda + * Added SHPCalculateExtents(). + * + * Revision 1.5 1998/11/09 20:57:16 warmerda + * Altered SHPGetInfo() call. + * + * Revision 1.4 1998/11/09 20:19:33 warmerda + * Added 3D support, and use of SHPObject. + * + * Revision 1.3 1995/08/23 02:24:05 warmerda + * Added support for reading bounds. + * + * Revision 1.2 1995/08/04 03:17:39 warmerda + * Added header. + * + */ + +#include + +#ifdef USE_DBMALLOC +#include +#endif + +#ifdef __cplusplus +extern "C" { +#endif + +/************************************************************************/ +/* Configuration options. */ +/************************************************************************/ + +/* -------------------------------------------------------------------- */ +/* Should the DBFReadStringAttribute() strip leading and */ +/* trailing white space? */ +/* -------------------------------------------------------------------- */ +#define TRIM_DBF_WHITESPACE + +/* -------------------------------------------------------------------- */ +/* Should we write measure values to the Multipatch object? */ +/* Reportedly ArcView crashes if we do write it, so for now it */ +/* is disabled. */ +/* -------------------------------------------------------------------- */ +#define DISABLE_MULTIPATCH_MEASURE + +/************************************************************************/ +/* SHP Support. */ +/************************************************************************/ +typedef struct +{ + FILE *fpSHP; + FILE *fpSHX; + + int nShapeType; /* SHPT_* */ + + int nFileSize; /* SHP file */ + + int nRecords; + int nMaxRecords; + int *panRecOffset; + int *panRecSize; + + double adBoundsMin[4]; + double adBoundsMax[4]; + + int bUpdated; +} SHPInfo; + +typedef SHPInfo * SHPHandle; + +/* -------------------------------------------------------------------- */ +/* Shape types (nSHPType) */ +/* -------------------------------------------------------------------- */ +#define SHPT_NULL 0 +#define SHPT_POINT 1 +#define SHPT_ARC 3 +#define SHPT_POLYGON 5 +#define SHPT_MULTIPOINT 8 +#define SHPT_POINTZ 11 +#define SHPT_ARCZ 13 +#define SHPT_POLYGONZ 15 +#define SHPT_MULTIPOINTZ 18 +#define SHPT_POINTM 21 +#define SHPT_ARCM 23 +#define SHPT_POLYGONM 25 +#define SHPT_MULTIPOINTM 28 +#define SHPT_MULTIPATCH 31 + + +/* -------------------------------------------------------------------- */ +/* Part types - everything but SHPT_MULTIPATCH just uses */ +/* SHPP_RING. */ +/* -------------------------------------------------------------------- */ + +#define SHPP_TRISTRIP 0 +#define SHPP_TRIFAN 1 +#define SHPP_OUTERRING 2 +#define SHPP_INNERRING 3 +#define SHPP_FIRSTRING 4 +#define SHPP_RING 5 + +/* -------------------------------------------------------------------- */ +/* SHPObject - represents on shape (without attributes) read */ +/* from the .shp file. */ +/* -------------------------------------------------------------------- */ +typedef struct +{ + int nSHPType; + + int nShapeId; /* -1 is unknown/unassigned */ + + int nParts; + int *panPartStart; + int *panPartType; + + int nVertices; + double *padfX; + double *padfY; + double *padfZ; + double *padfM; + + double dfXMin; + double dfYMin; + double dfZMin; + double dfMMin; + + double dfXMax; + double dfYMax; + double dfZMax; + double dfMMax; +} SHPObject; + +/* -------------------------------------------------------------------- */ +/* SHP API Prototypes */ +/* -------------------------------------------------------------------- */ +SHPHandle SHPOpen( const char * pszShapeFile, const char * pszAccess ); +SHPHandle SHPCreate( const char * pszShapeFile, int nShapeType ); +void SHPGetInfo( SHPHandle hSHP, int * pnEntities, int * pnShapeType, + double * padfMinBound, double * padfMaxBound ); + +SHPObject *SHPReadObject( SHPHandle hSHP, int iShape ); +int SHPWriteObject( SHPHandle hSHP, int iShape, SHPObject * psObject ); + +void SHPDestroyObject( SHPObject * psObject ); +void SHPComputeExtents( SHPObject * psObject ); +SHPObject *SHPCreateObject( int nSHPType, int nShapeId, + int nParts, int * panPartStart, int * panPartType, + int nVertices, double * padfX, double * padfY, + double * padfZ, double * padfM ); +SHPObject *SHPCreateSimpleObject( int nSHPType, int nVertices, + double * padfX, double * padfY, double * padfZ ); + +void SHPClose( SHPHandle hSHP ); + +const char *SHPTypeName( int nSHPType ); +const char *SHPPartTypeName( int nPartType ); + +/* -------------------------------------------------------------------- */ +/* Shape quadtree indexing API. */ +/* -------------------------------------------------------------------- */ + +/* this can be two or four for binary or quad tree */ +#define MAX_SUBNODE 4 + +typedef struct shape_tree_node +{ + /* region covered by this node */ + double adfBoundsMin[4]; + double adfBoundsMax[4]; + + /* list of shapes stored at this node. The papsShapeObj pointers + or the whole list can be NULL */ + int nShapeCount; + int *panShapeIds; + SHPObject **papsShapeObj; + + int nSubNodes; + struct shape_tree_node *apsSubNode[MAX_SUBNODE]; + +} SHPTreeNode; + +typedef struct +{ + SHPHandle hSHP; + + int nMaxDepth; + int nDimension; + + SHPTreeNode *psRoot; +} SHPTree; + +SHPTree *SHPCreateTree( SHPHandle hSHP, int nDimension, int nMaxDepth, + double *padfBoundsMin, double *padfBoundsMax ); +void SHPDestroyTree( SHPTree * hTree ); + +int SHPWriteTree( SHPTree *hTree, const char * pszFilename ); +SHPTree SHPReadTree( const char * pszFilename ); + +int SHPTreeAddObject( SHPTree * hTree, SHPObject * psObject ); +int SHPTreeAddShapeId( SHPTree * hTree, SHPObject * psObject ); +int SHPTreeRemoveShapeId( SHPTree * hTree, int nShapeId ); + +void SHPTreeTrimExtraNodes( SHPTree * hTree ); + +int *SHPTreeFindLikelyShapes( SHPTree * hTree, + double * padfBoundsMin, + double * padfBoundsMax, + int * ); +int SHPCheckBoundsOverlap( double *, double *, double *, double *, int ); + +/************************************************************************/ +/* DBF Support. */ +/************************************************************************/ +typedef struct +{ + FILE *fp; + + int nRecords; + + int nRecordLength; + int nHeaderLength; + int nFields; + int *panFieldOffset; + int *panFieldSize; + int *panFieldDecimals; + char *pachFieldType; + + char *pszHeader; + + int nCurrentRecord; + int bCurrentRecordModified; + char *pszCurrentRecord; + + int bNoHeader; + int bUpdated; +} DBFInfo; + +typedef DBFInfo * DBFHandle; + +typedef enum { + FTString, + FTInteger, + FTDouble, + FTInvalid, +} DBFFieldType; + +#define XBASE_FLDHDR_SZ 32 + +DBFHandle DBFOpen( const char * pszDBFFile, const char * pszAccess ); +DBFHandle DBFCreate( const char * pszDBFFile ); + +int DBFGetFieldCount( DBFHandle psDBF ); +int DBFGetRecordCount( DBFHandle psDBF ); +int DBFAddField( DBFHandle hDBF, const char * pszFieldName, + DBFFieldType eType, int nWidth, int nDecimals ); + +DBFFieldType DBFGetFieldInfo( DBFHandle psDBF, int iField, + char * pszFieldName, + int * pnWidth, int * pnDecimals ); + +int DBFReadIntegerAttribute( DBFHandle hDBF, int iShape, int iField ); +double DBFReadDoubleAttribute( DBFHandle hDBF, int iShape, int iField ); +const char *DBFReadStringAttribute( DBFHandle hDBF, int iShape, int iField ); + +int DBFWriteIntegerAttribute( DBFHandle hDBF, int iShape, int iField, + int nFieldValue ); +int DBFWriteDoubleAttribute( DBFHandle hDBF, int iShape, int iField, + double dFieldValue ); +int DBFWriteStringAttribute( DBFHandle hDBF, int iShape, int iField, + const char * pszFieldValue ); + +const char *DBFReadTuple(DBFHandle psDBF, int hEntity ); +int DBFWriteTuple(DBFHandle psDBF, int hEntity, void * pRawTuple ); + +DBFHandle DBFCloneEmpty(DBFHandle psDBF, const char * pszFilename ); + +void DBFClose( DBFHandle hDBF ); + +#ifdef __cplusplus +} +#endif + +#endif /* ndef _SHAPEFILE_H_INCLUDED */ diff --git a/loader/shp2pgsql.c b/loader/shp2pgsql.c new file mode 100644 index 000000000..ef429ae4c --- /dev/null +++ b/loader/shp2pgsql.c @@ -0,0 +1,858 @@ +//compile line for Solaris +//gcc -g pop.c ../shapelib-1.2.8/shpopen.o ../shapelib-1.2.8/dbfopen.o -o pop + +// usage: pop [-d || -a || -c] [-dump] | psql -h -d -p ... +// -d: drops the table , then recreates it and populates it with current shape file data +// -a: appends shape file into current table, must be excatly the same table schema +// -c: creates a new table and populates it, this is the default if you don't specify any options + +// -dump : use postgresql dump format (defaults to sql insert statements) + +// Using shapelib 1.2.8, this program reads in shape files and processes it's contents +// into a Insert statements which can be easily piped into a database frontend. +// Specifically designed to insert type 'geometry' (a custom written PostgreSQL type) +// for the shape files and PostgreSQL standard types for all attributes of the entity. +// +// Basically the program determines which type of shape (currently supports: 2d points,2d lines,2d +// polygons,3d points, and 3d lines) is in the file and takes appropriate action to read out the attributes. + + + +// BUGS: +// possible: no row # for polygons? + +#include "shapefil.h" +#include +#include + +typedef struct {double x, y;} Point; + +typedef struct Ring{ + Point *list; //list of points + struct Ring *next; + int n; //number of points in list +} Ring; + +int dump_format = 0; //0=insert statements, 1 = dump + +int Insert_attributes(DBFHandle hDBFHandle, char **ARGV, int row); + +char *make_good_string(char *str) +{ + //find all the tabs and make them \s + // + // 1. find # of tabs + // 2. make new string + // + // we dont escape already escaped tabs + + char *result; + char *str2; + char *start,*end; + int num_tabs = 0; + + str2 = str; + + while (str2 = strchr(str2, '\t') ) + { + if ( (str2 == str) || (str2[-1] != '\\') ) //the previous char isnt a \ + num_tabs ++; + str2++; + } + if (num_tabs == 0) + return str; + + result =(char *) malloc ( strlen(str) + num_tabs+1); + memset(result,0, strlen(str) + num_tabs+1 ); + start = str; + + while(end = strchr(start,'\t')) + { + if ( (end == str) || (end[-1] != '\\' ) ) //the previous char isnt a \ + { + strncat(result, start, (end-start)); + strcat(result,"\\\t"); + start = end +1; + } + else + { + strncat(result, start, (end-start)); + strcat(result,"\t"); + start = end +1; + } + } + strcat(result,start); + return result; + +} + +char *protect_quotes_string(char *str) +{ + //find all the tabs and make them \s + // + // 1. find # of quotes + // 2. make new string + + char *result; + char *str2; + char *start,*end; + int num_tabs = 0; + + str2 = str; + + while (str2 = strchr(str2, '\'') ) + { + if ( (str2 == str) || (str2[-1] != '\\') ) //the previous char isnt a \ + num_tabs ++; + str2++; + } + if (num_tabs == 0) + return str; + + result =(char *) malloc ( strlen(str) + num_tabs+1); + memset(result,0, strlen(str) + num_tabs+1 ); + start = str; + + while(end = strchr(start,'\'')) + { + if ( (end == str) || (end[-1] != '\\' ) ) //the previous char isnt a \ + { + + strncat(result, start, (end-start)); + strcat(result,"\\'"); + start = end +1; + } + else + { + strncat(result, start, (end-start)); + strcat(result,"'"); + start = end +1; + } + + } + strcat(result,start); + return result; +} + + + +// PIP(): crossing number test for a point in a polygon +// input: P = a point, +// V[] = vertex points of a polygon V[n+1] with V[n]=V[0] +// returns: 0 = outside, 1 = inside +int PIP( Point P, Point* V, int n ) +{ + int cn = 0; // the crossing number counter + int i; + // loop through all edges of the polygon + for (i=0; i P.y)) // an upward crossing + || ((V[i].y > P.y) && (V[i+1].y <= P.y))) { // a downward crossing + double vt = (float)(P.y - V[i].y) / (V[i+1].y - V[i].y); + if (P.x < V[i].x + vt * (V[i+1].x - V[i].x)) // P.x < intersect + ++cn; // a valid crossing of y=P.y right of P.x + } + } + return (cn&1); // 0 if even (out), and 1 if odd (in) + +} + + +//This function basically deals with the polygon case. + //it sorts the polys in order of outer,inner,inner, so that inners always come after outers they are within + //return value is the number of rings seen so far, used to keep id's unique. +int ring_check(SHPObject* obj, char **ARGV, int rings, DBFHandle hDBFHandle){ + Point pt,pt2; + Ring *Poly; + Ring *temp; + Ring **Outer; //pointer to a list of outer polygons + Ring **Inner; //pointer to a list of inner polygons + int out_index,in_index,indet_index; //indexes to the lists **Outer and **Inner + int in,temp2; + int u,i,N,n,new_outer; + int next_ring; //the index of the panPartStart list + double area; + + //initialize all counters/indexes + out_index=0; + in_index=0; + indet_index=0; + area=0; + n=0; + i=0; + N = obj->nVertices; + + if(obj->nParts >1){ + next_ring = 1;//if there is more than one part, set the next_ring index to one + }else{ + next_ring = -99; + } + + + //allocate initial pointer memory + Outer = (Ring**)malloc(sizeof(Ring*)*obj->nParts); + Inner = (Ring**)malloc(sizeof(Ring*)*obj->nParts); + Poly = (Ring*)malloc(sizeof(Ring)); + Poly->list = (Point*)malloc(sizeof(Point)*N); + Poly->next = NULL; + + + for (u=0;upanPartStart[next_ring] )) || u==N-1){ + //check if a ring is clockwise(outer) or not(inner) by getting positive(inner) or negative(outer) area. + //'area' is actually twice actual polygon area so divide by 2, not that it matters but in case we use it latter... + area = area/2.0; + if(area < 0.0 || obj->nParts ==1){ + + //An outer ring so fill in the last point then put it in the 'Outer' list + Poly->list[n].x = obj->padfX[u]; //the polygon is ended with it's first co-ordinates reused + Poly->list[n].y = obj->padfY[u]; + Poly->n = n+1; + Outer[out_index] = Poly; + out_index++; + + //allocate memory to start building the next ring + Poly = (Ring*)malloc(sizeof(Ring)); + + //temp2 is the number of points in the list of the next ring + //determined so that we can allocate the right amount of mem 6 lines down + if((next_ring + 1) == obj->nParts){ + temp2 = N; + }else{ + temp2 = obj->panPartStart[next_ring+1] - obj->panPartStart[next_ring]; + } + Poly->list = (Point*)malloc(sizeof(Point)*temp2); + Poly->next = NULL;//make sure to make to initiale next to null or you never know when the list ends + //this never used to be here and was a pain in the ass bug to find... + + n=0;//set your count of what point you are at in the current ring back to 0 + + }else{ + + Poly->list[n].x = obj->padfX[u]; //the polygon ends with it's first co-ordinates reused + Poly->list[n].y = obj->padfY[u]; + Poly->n = n+1; + + Inner[in_index] = Poly; + in_index++; + + Poly = (Ring*)malloc(sizeof(Ring)); + temp2 = N; + if((next_ring + 1) == obj->nParts){ + }else{ + temp2 = obj->panPartStart[next_ring+1] - obj->panPartStart[next_ring]; + } + + + //printf("temp2 -> %d for the list in the loop parts = %d N= %d\n",temp2,obj->nParts,N); + Poly->list = (Point*)malloc(sizeof(Point)*temp2); + Poly->next = NULL; + + n=0; + //printf("pt is ( %g %g)\n Poly starts with (%e %e)\n",obj->padfX[u],obj->padfY[u],Outer[0].x,Outer[0].y); + + } + area=0.0; + if((next_ring + 1) == obj->nParts){ + //printf("go to end of N\n"); + }else{ + next_ring++; + } + }else{ + + //printf(" x--- %g , y--- %g u=%d\n",obj->padfX[u],obj->padfY[u],u); + Poly->list[n].x = obj->padfX[u]; + Poly->list[n].y = obj->padfY[u]; + n++; + area += (obj->padfX[u] * obj->padfY[u+1]) - (obj->padfY[u] * obj->padfX[u+1]); //calculate the area + + } + } + + + +//Put the inner rings into the list of the outer rings of which they are within + for(u=0; u < in_index; u++){ + pt.x = Inner[u]->list[0].x; + pt.y = Inner[u]->list[0].y; + + pt2.x = Inner[u]->list[1].x; + pt2.y = Inner[u]->list[1].y; + for(i=0;i< out_index; i++){ + in = PIP(pt,Outer[i]->list,Outer[i]->n); + if(in==1 && PIP(pt2,Outer[i]->list,Outer[i]->n)){ + Poly = Outer[i]; + while(Poly->next != NULL){ + Poly = Poly->next; + } + Poly->next = Inner[u]; + break; + } + } + //if the ring wasn't within any outer rings, assume it is a new outer ring + + if(i == out_index){ + Outer[out_index] = Inner[u]; + out_index++; + } + } + + //start spitting out the sql for ordered entities now. + + if (dump_format) + { + printf("%i\tMULTIPOLYGON(",rings ); + } + else + printf("\nInsert into %s values('%i','MULTIPOLYGON(",ARGV[2],rings); + rings++; + for(u=0; u < out_index; u++){ + Poly = Outer[u]; + if(u==0){ + printf("("); + }else{ + printf(",("); + } + while(Poly != NULL){ + for(i=0;in;i++){ + if(i==0){ + if(Poly != Outer[u]){ + printf(","); + } + printf("(%.15g %.15g ",Poly->list[i].x,Poly->list[i].y); + }else{ + printf(",%.15g %.15g ",Poly->list[i].x,Poly->list[i].y); + } + } + printf(")"); + temp = Poly; + Poly = Poly->next; + free(temp->list); + free(temp); + } + printf(")"); + } + if (dump_format) + { + printf(")"); + } + else + printf(")'"); + Insert_attributes(hDBFHandle, ARGV,rings-1); + + + free(Outer); + free(Inner); + free(Poly); + + + return rings; +} + + + +//Insert the attributes from the correct row of dbf file + +int Insert_attributes(DBFHandle hDBFHandle, char **ARGV, int row){ + + int i,num_fields; + + + num_fields = DBFGetFieldCount( hDBFHandle ); + + + for( i = 0; i < num_fields; i++ ){ + if (dump_format) + printf("\t%s",make_good_string(DBFReadStringAttribute( hDBFHandle, row, i )) ); + else + printf(",'%s'",protect_quotes_string(DBFReadStringAttribute( hDBFHandle, row, i )) ); + } + if (!(dump_format) ) + printf (");\n"); + else + printf("\n"); + + return 1; +} + + + + +// main() +// usage: pop
[-d || -a || -c] [-dump]| psql -h -d -p ... +// -d: drops the table , then recreates it and populates it with current shape file data +// -a: appends shape file into current table, must be excatly the same table schema +// -c: creates a new table and populates it, this is the default if you don't specify any options + +// -dump : use postgresql dump format (defaults to sql insert statements) + +// Using shapelib 1.2.8, this program reads in shape files and processes it's contents +// into a Insert statements which can be easily piped into a database frontend. +// Specifically designed to insert type 'geometry' (a custom written PostgreSQL type) +// for the shape files and PostgreSQL standard types for all attributes of the entity. +// +// Basically the program determines which type of shape (currently supports: 2d points,2d lines,2d +// polygons,3d points, and 3d lines) is in the file and takes appropriate action to read out the attributes. + + +main (int ARGC, char **ARGV) +{ + SHPHandle hSHPHandle; + DBFHandle hDBFHandle; + int num_fields,num_records,begin,trans; + int num_entities, phnshapetype,next_ring; + double padminbound[8], padmaxbound[8]; + int u,j,tot_rings; + SHPObject *obj; + char name[12]; + char opt; + DBFFieldType type; + + //display proper usage if incorrect number of arguments given + if (ARGC != 3 && ARGC != 4 && ARGC != 5) + { + printf ("usage: pop
[-d || -a || -c] [-dump]\n"); + printf(" usage: pop
[-d || -a || -c] [-dump]| psql -h -d -p ...\n"); + printf(" -d: drops the table , then recreates it and populates it with current shape file data\n"); + printf(" -a: appends shape file into current table, must be excatly the same table schema\n"); + printf(" -c: creates a new table and populates it, this is the default if you don't specify any options\n"); + printf(" -dump : use postgresql dump format (defaults to sql insert statements)\n"); + + exit (-1); + } + + opt = 'c'; + + if(ARGC ==4){ + if(strcmp(ARGV[3], "-d")==0){ + opt = 'd'; + }else if(strcmp(ARGV[3], "-c")==0){ + opt = 'c'; + }else if(strcmp(ARGV[3], "-a")==0){ + opt = 'a'; + }else if (strcmp(ARGV[3],"-dump") == 0) { + dump_format = 1; + }else{ + printf("option %s is not a valid option, use -a, -d, or -c\n",ARGV[3]); + exit(-1); + } + }else{ + opt = 'c'; + } + + if (ARGC==5) + { + if (strcmp(ARGV[4],"-dump") == 0) + { + dump_format = 1; + } + } + + + //Open the shp and dbf files + hSHPHandle = SHPOpen( ARGV[1], "rb" ); + hDBFHandle = DBFOpen( ARGV[1], "rb" ); + if (hSHPHandle == NULL || hDBFHandle == NULL){ + printf ("shape is null\n"); + exit(-1); + } + + if(opt == 'd'){ + //-------------------------Drop the table-------------------------------- + //drop the table given + printf("drop table %s;",ARGV[2]); + + } + + + + if(opt == 'c' || opt == 'd'){ //if opt is 'a' do nothing, go straight to making inserts + + //-------------------------Create the table-------------------------------- + //create a table for inserting the shapes into with appropriate columns and types + + printf("create table %s (geoid int4, geo_value geometry ",ARGV[2]); + + num_fields = DBFGetFieldCount( hDBFHandle ); + num_records = DBFGetRecordCount(hDBFHandle); + + for(j=0;jpachFieldType[j] == 'D' ){ + printf ("varchar(8)");//date type is not supported in API so check for it explicity before the api call. + }else{ + if(type == FTString){ + printf ("varchar"); + }else if(type == FTInteger){ + printf ("int4"); + }else if(type == FTDouble){ + printf ("float8"); + }else{ + printf ("Invalid type in DBF file"); + } + } + } + printf (");\n"); + //finished creating the table + } + + if (dump_format) + { + printf("COPY \"%s\" from stdin;\n",ARGV[2]); + } + + SHPGetInfo( hSHPHandle, &num_entities, &phnshapetype, &padminbound[0], &padmaxbound[0]); + obj = SHPReadObject(hSHPHandle,0); + trans=0; + //Determine what type of shape is in the file and do appropriate processing + if (obj->nVertices == 0) + { + if (dump_format) + { + printf("\\N"); + Insert_attributes(hDBFHandle, ARGV,j);//add the attributes of each shape to the insert statement + SHPDestroyObject(obj); + } + else + { + printf("NULL"); + Insert_attributes(hDBFHandle, ARGV,j);//add the attributes of each shape to the insert statement + SHPDestroyObject(obj); + } + } + else + { + if( obj->nSHPType == 5 ){ + //--------------------------------------------------------------------------------- + //---------POLYGONS---------------------------------------------------------------- + + // sorts of all the rings so that they are outer,inner,iner,outer,inner... + // with the inner ones coming after the outer ones they are within spatially + + tot_rings = 0; + + + //go through each entity and call ring_check() to sort the rings and print out the sql statement + // keep track of total number of inserts in tot_rings so + // you can pass it to the function for the next entity + for (j=0;jnSHPType == 1){ + //--------------------------------------------------------------------- + //----------POINTS----------------------------------------------------- + + for (j=0;jnVertices; u++){ + if (u>0){ + printf(",%.15g %.15g",obj->padfX[u],obj->padfY[u]); + }else{ + printf("%.15g %.15g",obj->padfX[u],obj->padfY[u]); + } + } + + if (dump_format) + { + printf(")"); + + } + else + printf(")'"); + Insert_attributes(hDBFHandle, ARGV,j); //add the attributes for each entity to the insert statement + + SHPDestroyObject(obj); + } + if (!(dump_format) ) + printf("end;"); //End the last transaction + + }else if( obj->nSHPType == 3){ + //------------------------------------------------------------------------ + //--------ARCs / LINES---------------------------------------------------- + + begin=0;//initialize the begin flag + + for (j=0;jnParts >1){ + next_ring = 1;//if there is more than one part, set the next_ring index to one + }else{ + next_ring = -99; + } + + //for each vertice write out the coordinates in the insert statement, when there is a new line + //you must end the brackets and start new ones etc. + for (u=0;unVertices; u++){ + + //check if the next vertice is the start of a new line + // printf("\n\nu+1 = %d, next_ring = %d index = %d\n",u+1,next_ring,obj->panPartStart[next_ring]); + + if(next_ring==-99 && obj->nVertices ==1){ + printf("(%.15g %.15g )",obj->padfX[u],obj->padfY[u]); + }else if((next_ring != -99)&& (begin==1) && (u+1 == obj->panPartStart[next_ring]) ){ + printf("(%.15g %.15g )",obj->padfX[u],obj->padfY[u]); + next_ring++; + begin=1; + }else if(((next_ring != -99) && (u+1 == obj->panPartStart[next_ring] )) || u==(obj->nVertices-1) ){ + printf(",%.15g %.15g ",obj->padfX[u],obj->padfY[u]); + printf(")"); + + next_ring++; + begin=1;//flag the fact that you area at a new line next time through the loop + }else{ + if (u==0 || begin==1){ //if you are at the begging of a new line add comma and brackets + if(u!=0) printf(","); + printf("(%.15g %.15g ",obj->padfX[u],obj->padfY[u]); + begin=0; + }else{ + printf(",%.15g %.15g ",obj->padfX[u],obj->padfY[u]); + } + } + } + + if (dump_format) + { + printf(")"); + } + else + printf(")'"); + Insert_attributes(hDBFHandle, ARGV,j); //add the attributes of each shape to the insert statement + + + SHPDestroyObject(obj); + } + + if (!(dump_format) ) + printf("end;");//end the last transaction + + + + }else if( obj->nSHPType == 13){ + //--------------------------------------------------------------------- + //------PolyLineZ(3D lines)-------------------------------------------- + + begin=0;//initialize the begin flag + + for (j=0;jnParts >1){ + next_ring = 1;//if there is more than one part, set the next_ring index to one + }else{ + next_ring = -99; + } + + //for each vertice write out the coordinates in the insert statement, when there is a new line + //you must end the brackets and start new ones etc. + for (u=0;unVertices; u++){ + + //check if the next vertice is the start of a new line + if(((next_ring != -99) && (u+1 == obj->panPartStart[next_ring] )) || u==(obj->nVertices-1) ){ + printf(",%.15g %.15g ",obj->padfX[u],obj->padfY[u]); + printf(")"); + next_ring++; + begin =1;//flag the fact that you area at a new line next time through the loop + }else{ + if (u==0 || begin==1){ + if(u!=0) printf(","); + printf("(%.15g %.15g %.15g ",obj->padfX[u],obj->padfY[u],obj->padfZ[u]); + begin=0; + }else{ + printf(",%.15g %.15g %.15g ",obj->padfX[u],obj->padfY[u],obj->padfZ[u]); + } + } + } + + if (dump_format) + { + printf(")"); + } + else + printf(")'"); + Insert_attributes(hDBFHandle, ARGV,j);//add the attributes of each shape to the insert statement + + + SHPDestroyObject(obj); + } + + if (!(dump_format) ) + printf("end;");//close last transaction + + + }else if( obj->nSHPType == 11){ + //--------------------------------------------------------------------------- + //------POINTZ (3D POINTS)--------------------------------------------------- + + for (j=0;jnVertices; u++){ + if (u>0){ + printf(",%.15g %.15g %.15g",obj->padfX[u],obj->padfY[u],obj->padfZ[u]); + }else{ + printf("%.15g %.15g %.15g",obj->padfX[u],obj->padfY[u],obj->padfZ[u]); + } + } + + if (dump_format) + { + printf(")"); + } + else + printf(")'"); + Insert_attributes(hDBFHandle, ARGV,j);//add the attributes of each shape to the insert statement + + + SHPDestroyObject(obj); + } + + + + if (!(dump_format) ) + printf("end;");//end the last transaction + }else{ + printf (""); + printf ("\n\n**** Type is NOT SUPPORTED, type id = %d ****\n\n",obj->nSHPType); + //print out what type the file is and that it is not supported + + }//end the if statement for shape types + + } + +if ((dump_format) ) + printf("\\.\n"); +}//end main() diff --git a/loader/shpopen.c b/loader/shpopen.c new file mode 100644 index 000000000..d7c475f55 --- /dev/null +++ b/loader/shpopen.c @@ -0,0 +1,1623 @@ +/****************************************************************************** + * $Id$ + * + * Project: Shapelib + * Purpose: Implementation of core Shapefile read/write functions. + * Author: Frank Warmerdam, warmerda@home.com + * + ****************************************************************************** + * Copyright (c) 1999, Frank Warmerdam + * + * This software is available under the following "MIT Style" license, + * or at the option of the licensee under the LGPL (see LICENSE.LGPL). This + * option is discussed in more detail in shapelib.html. + * + * -- + * + * Permission is hereby granted, free of charge, to any person obtaining a + * copy of this software and associated documentation files (the "Software"), + * to deal in the Software without restriction, including without limitation + * the rights to use, copy, modify, merge, publish, distribute, sublicense, + * and/or sell copies of the Software, and to permit persons to whom the + * Software is furnished to do so, subject to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS + * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER + * DEALINGS IN THE SOFTWARE. + ****************************************************************************** + * + * $Log$ + * Revision 1.1 2001/07/18 21:42:12 pramsey + * Initial add of the data loader code. + * + * Revision 1.26 2000/02/16 16:03:51 warmerda + * added null shape support + * + * Revision 1.25 1999/12/15 13:47:07 warmerda + * Fixed record size settings in .shp file (was 4 words too long) + * Added stdlib.h. + * + * Revision 1.24 1999/11/05 14:12:04 warmerda + * updated license terms + * + * Revision 1.23 1999/07/27 00:53:46 warmerda + * added support for rewriting shapes + * + * Revision 1.22 1999/06/11 19:19:11 warmerda + * Cleanup pabyRec static buffer on SHPClose(). + * + * Revision 1.21 1999/06/02 14:57:56 kshih + * Remove unused variables + * + * Revision 1.20 1999/04/19 21:04:17 warmerda + * Fixed syntax error. + * + * Revision 1.19 1999/04/19 21:01:57 warmerda + * Force access string to binary in SHPOpen(). + * + * Revision 1.18 1999/04/01 18:48:07 warmerda + * Try upper case extensions if lower case doesn't work. + * + * Revision 1.17 1998/12/31 15:29:39 warmerda + * Disable writing measure values to multipatch objects if + * DISABLE_MULTIPATCH_MEASURE is defined. + * + * Revision 1.16 1998/12/16 05:14:33 warmerda + * Added support to write MULTIPATCH. Fixed reading Z coordinate of + * MULTIPATCH. Fixed record size written for all feature types. + * + * Revision 1.15 1998/12/03 16:35:29 warmerda + * r+b is proper binary access string, not rb+. + * + * Revision 1.14 1998/12/03 15:47:56 warmerda + * Fixed setting of nVertices in SHPCreateObject(). + * + * Revision 1.13 1998/12/03 15:33:54 warmerda + * Made SHPCalculateExtents() separately callable. + * + * Revision 1.12 1998/11/11 20:01:50 warmerda + * Fixed bug writing ArcM/Z, and PolygonM/Z for big endian machines. + * + * Revision 1.11 1998/11/09 20:56:44 warmerda + * Fixed up handling of file wide bounds. + * + * Revision 1.10 1998/11/09 20:18:51 warmerda + * Converted to support 3D shapefiles, and use of SHPObject. + * + * Revision 1.9 1998/02/24 15:09:05 warmerda + * Fixed memory leak. + * + * Revision 1.8 1997/12/04 15:40:29 warmerda + * Fixed byte swapping of record number, and record length fields in the + * .shp file. + * + * Revision 1.7 1995/10/21 03:15:58 warmerda + * Added support for binary file access, the magic cookie 9997 + * and tried to improve the int32 selection logic for 16bit systems. + * + * Revision 1.6 1995/09/04 04:19:41 warmerda + * Added fix for file bounds. + * + * Revision 1.5 1995/08/25 15:16:44 warmerda + * Fixed a couple of problems with big endian systems ... one with bounds + * and the other with multipart polygons. + * + * Revision 1.4 1995/08/24 18:10:17 warmerda + * Switch to use SfRealloc() to avoid problems with pre-ANSI realloc() + * functions (such as on the Sun). + * + * Revision 1.3 1995/08/23 02:23:15 warmerda + * Added support for reading bounds, and fixed up problems in setting the + * file wide bounds. + * + * Revision 1.2 1995/08/04 03:16:57 warmerda + * Added header. + * + */ + +static char rcsid[] = + "$Id$"; + +#include "shapefil.h" + +#include +#include +#include +#include + +typedef unsigned char uchar; + +#if UINT_MAX == 65535 +typedef long int32; +#else +typedef int int32; +#endif + +#ifndef FALSE +# define FALSE 0 +# define TRUE 1 +#endif + +#define ByteCopy( a, b, c ) memcpy( b, a, c ) +#ifndef MAX +# define MIN(a,b) ((ab) ? a : b) +#endif + +static int bBigEndian; +static uchar *pabyRec = NULL; +static int nBufSize = 0; + + +/************************************************************************/ +/* SwapWord() */ +/* */ +/* Swap a 2, 4 or 8 byte word. */ +/************************************************************************/ + +static void SwapWord( int length, void * wordP ) + +{ + int i; + uchar temp; + + for( i=0; i < length/2; i++ ) + { + temp = ((uchar *) wordP)[i]; + ((uchar *)wordP)[i] = ((uchar *) wordP)[length-i-1]; + ((uchar *) wordP)[length-i-1] = temp; + } +} + +/************************************************************************/ +/* SfRealloc() */ +/* */ +/* A realloc cover function that will access a NULL pointer as */ +/* a valid input. */ +/************************************************************************/ + +static void * SfRealloc( void * pMem, int nNewSize ) + +{ + if( pMem == NULL ) + return( (void *) malloc(nNewSize) ); + else + return( (void *) realloc(pMem,nNewSize) ); +} + +/************************************************************************/ +/* SHPWriteHeader() */ +/* */ +/* Write out a header for the .shp and .shx files as well as the */ +/* contents of the index (.shx) file. */ +/************************************************************************/ + +static void SHPWriteHeader( SHPHandle psSHP ) + +{ + uchar abyHeader[100]; + int i; + int32 i32; + double dValue; + int32 *panSHX; + +/* -------------------------------------------------------------------- */ +/* Prepare header block for .shp file. */ +/* -------------------------------------------------------------------- */ + for( i = 0; i < 100; i++ ) + abyHeader[i] = 0; + + abyHeader[2] = 0x27; /* magic cookie */ + abyHeader[3] = 0x0a; + + i32 = psSHP->nFileSize/2; /* file size */ + ByteCopy( &i32, abyHeader+24, 4 ); + if( !bBigEndian ) SwapWord( 4, abyHeader+24 ); + + i32 = 1000; /* version */ + ByteCopy( &i32, abyHeader+28, 4 ); + if( bBigEndian ) SwapWord( 4, abyHeader+28 ); + + i32 = psSHP->nShapeType; /* shape type */ + ByteCopy( &i32, abyHeader+32, 4 ); + if( bBigEndian ) SwapWord( 4, abyHeader+32 ); + + dValue = psSHP->adBoundsMin[0]; /* set bounds */ + ByteCopy( &dValue, abyHeader+36, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+36 ); + + dValue = psSHP->adBoundsMin[1]; + ByteCopy( &dValue, abyHeader+44, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+44 ); + + dValue = psSHP->adBoundsMax[0]; + ByteCopy( &dValue, abyHeader+52, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+52 ); + + dValue = psSHP->adBoundsMax[1]; + ByteCopy( &dValue, abyHeader+60, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+60 ); + + dValue = psSHP->adBoundsMin[2]; /* z */ + ByteCopy( &dValue, abyHeader+68, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+68 ); + + dValue = psSHP->adBoundsMax[2]; + ByteCopy( &dValue, abyHeader+76, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+76 ); + + dValue = psSHP->adBoundsMin[3]; /* m */ + ByteCopy( &dValue, abyHeader+84, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+84 ); + + dValue = psSHP->adBoundsMax[3]; + ByteCopy( &dValue, abyHeader+92, 8 ); + if( bBigEndian ) SwapWord( 8, abyHeader+92 ); + +/* -------------------------------------------------------------------- */ +/* Write .shp file header. */ +/* -------------------------------------------------------------------- */ + fseek( psSHP->fpSHP, 0, 0 ); + fwrite( abyHeader, 100, 1, psSHP->fpSHP ); + +/* -------------------------------------------------------------------- */ +/* Prepare, and write .shx file header. */ +/* -------------------------------------------------------------------- */ + i32 = (psSHP->nRecords * 2 * sizeof(int32) + 100)/2; /* file size */ + ByteCopy( &i32, abyHeader+24, 4 ); + if( !bBigEndian ) SwapWord( 4, abyHeader+24 ); + + fseek( psSHP->fpSHX, 0, 0 ); + fwrite( abyHeader, 100, 1, psSHP->fpSHX ); + +/* -------------------------------------------------------------------- */ +/* Write out the .shx contents. */ +/* -------------------------------------------------------------------- */ + panSHX = (int32 *) malloc(sizeof(int32) * 2 * psSHP->nRecords); + + for( i = 0; i < psSHP->nRecords; i++ ) + { + panSHX[i*2 ] = psSHP->panRecOffset[i]/2; + panSHX[i*2+1] = psSHP->panRecSize[i]/2; + if( !bBigEndian ) SwapWord( 4, panSHX+i*2 ); + if( !bBigEndian ) SwapWord( 4, panSHX+i*2+1 ); + } + + fwrite( panSHX, sizeof(int32) * 2, psSHP->nRecords, psSHP->fpSHX ); + + free( panSHX ); +} + +/************************************************************************/ +/* SHPOpen() */ +/* */ +/* Open the .shp and .shx files based on the basename of the */ +/* files or either file name. */ +/************************************************************************/ + +SHPHandle SHPOpen( const char * pszLayer, const char * pszAccess ) + +{ + char *pszFullname, *pszBasename; + SHPHandle psSHP; + + uchar *pabyBuf; + int i; + double dValue; + +/* -------------------------------------------------------------------- */ +/* Ensure the access string is one of the legal ones. We */ +/* ensure the result string indicates binary to avoid common */ +/* problems on Windows. */ +/* -------------------------------------------------------------------- */ + if( strcmp(pszAccess,"rb+") == 0 || strcmp(pszAccess,"r+b") == 0 + || strcmp(pszAccess,"r+") == 0 ) + pszAccess = "r+b"; + else + pszAccess = "rb"; + +/* -------------------------------------------------------------------- */ +/* Establish the byte order on this machine. */ +/* -------------------------------------------------------------------- */ + i = 1; + if( *((uchar *) &i) == 1 ) + bBigEndian = FALSE; + else + bBigEndian = TRUE; + +/* -------------------------------------------------------------------- */ +/* Initialize the info structure. */ +/* -------------------------------------------------------------------- */ + psSHP = (SHPHandle) malloc(sizeof(SHPInfo)); + + psSHP->bUpdated = FALSE; + +/* -------------------------------------------------------------------- */ +/* Compute the base (layer) name. If there is any extension */ +/* on the passed in filename we will strip it off. */ +/* -------------------------------------------------------------------- */ + pszBasename = (char *) malloc(strlen(pszLayer)+5); + strcpy( pszBasename, pszLayer ); + for( i = strlen(pszBasename)-1; + i > 0 && pszBasename[i] != '.' && pszBasename[i] != '/' + && pszBasename[i] != '\\'; + i-- ) {} + + if( pszBasename[i] == '.' ) + pszBasename[i] = '\0'; + +/* -------------------------------------------------------------------- */ +/* Open the .shp and .shx files. Note that files pulled from */ +/* a PC to Unix with upper case filenames won't work! */ +/* -------------------------------------------------------------------- */ + pszFullname = (char *) malloc(strlen(pszBasename) + 5); + sprintf( pszFullname, "%s.shp", pszBasename ); + psSHP->fpSHP = fopen(pszFullname, pszAccess ); + if( psSHP->fpSHP == NULL ) + { + sprintf( pszFullname, "%s.SHP", pszBasename ); + psSHP->fpSHP = fopen(pszFullname, pszAccess ); + } + + if( psSHP->fpSHP == NULL ) + return( NULL ); + + sprintf( pszFullname, "%s.shx", pszBasename ); + psSHP->fpSHX = fopen(pszFullname, pszAccess ); + if( psSHP->fpSHX == NULL ) + { + sprintf( pszFullname, "%s.SHX", pszBasename ); + psSHP->fpSHX = fopen(pszFullname, pszAccess ); + } + + if( psSHP->fpSHX == NULL ) + return( NULL ); + + free( pszFullname ); + free( pszBasename ); + +/* -------------------------------------------------------------------- */ +/* Read the file size from the SHP file. */ +/* -------------------------------------------------------------------- */ + pabyBuf = (uchar *) malloc(100); + fread( pabyBuf, 100, 1, psSHP->fpSHP ); + + psSHP->nFileSize = (pabyBuf[24] * 256 * 256 * 256 + + pabyBuf[25] * 256 * 256 + + pabyBuf[26] * 256 + + pabyBuf[27]) * 2; + +/* -------------------------------------------------------------------- */ +/* Read SHX file Header info */ +/* -------------------------------------------------------------------- */ + fread( pabyBuf, 100, 1, psSHP->fpSHX ); + + if( pabyBuf[0] != 0 + || pabyBuf[1] != 0 + || pabyBuf[2] != 0x27 + || (pabyBuf[3] != 0x0a && pabyBuf[3] != 0x0d) ) + { + fclose( psSHP->fpSHP ); + fclose( psSHP->fpSHX ); + free( psSHP ); + + return( NULL ); + } + + psSHP->nRecords = pabyBuf[27] + pabyBuf[26] * 256 + + pabyBuf[25] * 256 * 256 + pabyBuf[24] * 256 * 256 * 256; + psSHP->nRecords = (psSHP->nRecords*2 - 100) / 8; + + psSHP->nShapeType = pabyBuf[32]; + +/* -------------------------------------------------------------------- */ +/* Read the bounds. */ +/* -------------------------------------------------------------------- */ + if( bBigEndian ) SwapWord( 8, pabyBuf+36 ); + memcpy( &dValue, pabyBuf+36, 8 ); + psSHP->adBoundsMin[0] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+44 ); + memcpy( &dValue, pabyBuf+44, 8 ); + psSHP->adBoundsMin[1] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+52 ); + memcpy( &dValue, pabyBuf+52, 8 ); + psSHP->adBoundsMax[0] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+60 ); + memcpy( &dValue, pabyBuf+60, 8 ); + psSHP->adBoundsMax[1] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+68 ); /* z */ + memcpy( &dValue, pabyBuf+68, 8 ); + psSHP->adBoundsMin[2] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+76 ); + memcpy( &dValue, pabyBuf+76, 8 ); + psSHP->adBoundsMax[2] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+84 ); /* z */ + memcpy( &dValue, pabyBuf+84, 8 ); + psSHP->adBoundsMin[3] = dValue; + + if( bBigEndian ) SwapWord( 8, pabyBuf+92 ); + memcpy( &dValue, pabyBuf+92, 8 ); + psSHP->adBoundsMax[3] = dValue; + + free( pabyBuf ); + +/* -------------------------------------------------------------------- */ +/* Read the .shx file to get the offsets to each record in */ +/* the .shp file. */ +/* -------------------------------------------------------------------- */ + psSHP->nMaxRecords = psSHP->nRecords; + + psSHP->panRecOffset = + (int *) malloc(sizeof(int) * MAX(1,psSHP->nMaxRecords) ); + psSHP->panRecSize = + (int *) malloc(sizeof(int) * MAX(1,psSHP->nMaxRecords) ); + + pabyBuf = (uchar *) malloc(8 * MAX(1,psSHP->nRecords) ); + fread( pabyBuf, 8, psSHP->nRecords, psSHP->fpSHX ); + + for( i = 0; i < psSHP->nRecords; i++ ) + { + int32 nOffset, nLength; + + memcpy( &nOffset, pabyBuf + i * 8, 4 ); + if( !bBigEndian ) SwapWord( 4, &nOffset ); + + memcpy( &nLength, pabyBuf + i * 8 + 4, 4 ); + if( !bBigEndian ) SwapWord( 4, &nLength ); + + psSHP->panRecOffset[i] = nOffset*2; + psSHP->panRecSize[i] = nLength*2; + } + free( pabyBuf ); + + return( psSHP ); +} + +/************************************************************************/ +/* SHPClose() */ +/* */ +/* Close the .shp and .shx files. */ +/************************************************************************/ + +void SHPClose(SHPHandle psSHP ) + +{ +/* -------------------------------------------------------------------- */ +/* Update the header if we have modified anything. */ +/* -------------------------------------------------------------------- */ + if( psSHP->bUpdated ) + { + SHPWriteHeader( psSHP ); + } + +/* -------------------------------------------------------------------- */ +/* Free all resources, and close files. */ +/* -------------------------------------------------------------------- */ + free( psSHP->panRecOffset ); + free( psSHP->panRecSize ); + + fclose( psSHP->fpSHX ); + fclose( psSHP->fpSHP ); + + free( psSHP ); + + if( pabyRec != NULL ) + { + free( pabyRec ); + pabyRec = NULL; + nBufSize = 0; + } +} + +/************************************************************************/ +/* SHPGetInfo() */ +/* */ +/* Fetch general information about the shape file. */ +/************************************************************************/ + +void SHPGetInfo(SHPHandle psSHP, int * pnEntities, int * pnShapeType, + double * padfMinBound, double * padfMaxBound ) + +{ + int i; + + if( pnEntities != NULL ) + *pnEntities = psSHP->nRecords; + + if( pnShapeType != NULL ) + *pnShapeType = psSHP->nShapeType; + + for( i = 0; i < 4; i++ ) + { + if( padfMinBound != NULL ) + padfMinBound[i] = psSHP->adBoundsMin[i]; + if( padfMaxBound != NULL ) + padfMaxBound[i] = psSHP->adBoundsMax[i]; + } +} + +/************************************************************************/ +/* SHPCreate() */ +/* */ +/* Create a new shape file and return a handle to the open */ +/* shape file with read/write access. */ +/************************************************************************/ + +SHPHandle SHPCreate( const char * pszLayer, int nShapeType ) + +{ + char *pszBasename, *pszFullname; + int i; + FILE *fpSHP, *fpSHX; + uchar abyHeader[100]; + int32 i32; + double dValue; + +/* -------------------------------------------------------------------- */ +/* Establish the byte order on this system. */ +/* -------------------------------------------------------------------- */ + i = 1; + if( *((uchar *) &i) == 1 ) + bBigEndian = FALSE; + else + bBigEndian = TRUE; + +/* -------------------------------------------------------------------- */ +/* Compute the base (layer) name. If there is any extension */ +/* on the passed in filename we will strip it off. */ +/* -------------------------------------------------------------------- */ + pszBasename = (char *) malloc(strlen(pszLayer)+5); + strcpy( pszBasename, pszLayer ); + for( i = strlen(pszBasename)-1; + i > 0 && pszBasename[i] != '.' && pszBasename[i] != '/' + && pszBasename[i] != '\\'; + i-- ) {} + + if( pszBasename[i] == '.' ) + pszBasename[i] = '\0'; + +/* -------------------------------------------------------------------- */ +/* Open the two files so we can write their headers. */ +/* -------------------------------------------------------------------- */ + pszFullname = (char *) malloc(strlen(pszBasename) + 5); + sprintf( pszFullname, "%s.shp", pszBasename ); + fpSHP = fopen(pszFullname, "wb" ); + if( fpSHP == NULL ) + return( NULL ); + + sprintf( pszFullname, "%s.shx", pszBasename ); + fpSHX = fopen(pszFullname, "wb" ); + if( fpSHX == NULL ) + return( NULL ); + + free( pszFullname ); + free( pszBasename ); + +/* -------------------------------------------------------------------- */ +/* Prepare header block for .shp file. */ +/* -------------------------------------------------------------------- */ + for( i = 0; i < 100; i++ ) + abyHeader[i] = 0; + + abyHeader[2] = 0x27; /* magic cookie */ + abyHeader[3] = 0x0a; + + i32 = 50; /* file size */ + ByteCopy( &i32, abyHeader+24, 4 ); + if( !bBigEndian ) SwapWord( 4, abyHeader+24 ); + + i32 = 1000; /* version */ + ByteCopy( &i32, abyHeader+28, 4 ); + if( bBigEndian ) SwapWord( 4, abyHeader+28 ); + + i32 = nShapeType; /* shape type */ + ByteCopy( &i32, abyHeader+32, 4 ); + if( bBigEndian ) SwapWord( 4, abyHeader+32 ); + + dValue = 0.0; /* set bounds */ + ByteCopy( &dValue, abyHeader+36, 8 ); + ByteCopy( &dValue, abyHeader+44, 8 ); + ByteCopy( &dValue, abyHeader+52, 8 ); + ByteCopy( &dValue, abyHeader+60, 8 ); + +/* -------------------------------------------------------------------- */ +/* Write .shp file header. */ +/* -------------------------------------------------------------------- */ + fwrite( abyHeader, 100, 1, fpSHP ); + +/* -------------------------------------------------------------------- */ +/* Prepare, and write .shx file header. */ +/* -------------------------------------------------------------------- */ + i32 = 50; /* file size */ + ByteCopy( &i32, abyHeader+24, 4 ); + if( !bBigEndian ) SwapWord( 4, abyHeader+24 ); + + fwrite( abyHeader, 100, 1, fpSHX ); + +/* -------------------------------------------------------------------- */ +/* Close the files, and then open them as regular existing files. */ +/* -------------------------------------------------------------------- */ + fclose( fpSHP ); + fclose( fpSHX ); + + return( SHPOpen( pszLayer, "r+b" ) ); +} + +/************************************************************************/ +/* _SHPSetBounds() */ +/* */ +/* Compute a bounds rectangle for a shape, and set it into the */ +/* indicated location in the record. */ +/************************************************************************/ + +static void _SHPSetBounds( uchar * pabyRec, SHPObject * psShape ) + +{ + ByteCopy( &(psShape->dfXMin), pabyRec + 0, 8 ); + ByteCopy( &(psShape->dfYMin), pabyRec + 8, 8 ); + ByteCopy( &(psShape->dfXMax), pabyRec + 16, 8 ); + ByteCopy( &(psShape->dfYMax), pabyRec + 24, 8 ); + + if( bBigEndian ) + { + SwapWord( 8, pabyRec + 0 ); + SwapWord( 8, pabyRec + 8 ); + SwapWord( 8, pabyRec + 16 ); + SwapWord( 8, pabyRec + 24 ); + } +} + +/************************************************************************/ +/* SHPComputeExtents() */ +/* */ +/* Recompute the extents of a shape. Automatically done by */ +/* SHPCreateObject(). */ +/************************************************************************/ + +void SHPComputeExtents( SHPObject * psObject ) + +{ + int i; + +/* -------------------------------------------------------------------- */ +/* Build extents for this object. */ +/* -------------------------------------------------------------------- */ + if( psObject->nVertices > 0 ) + { + psObject->dfXMin = psObject->dfXMax = psObject->padfX[0]; + psObject->dfYMin = psObject->dfYMax = psObject->padfY[0]; + psObject->dfZMin = psObject->dfZMax = psObject->padfZ[0]; + psObject->dfMMin = psObject->dfMMax = psObject->padfM[0]; + } + + for( i = 0; i < psObject->nVertices; i++ ) + { + psObject->dfXMin = MIN(psObject->dfXMin, psObject->padfX[i]); + psObject->dfYMin = MIN(psObject->dfYMin, psObject->padfY[i]); + psObject->dfZMin = MIN(psObject->dfZMin, psObject->padfZ[i]); + psObject->dfMMin = MIN(psObject->dfMMin, psObject->padfM[i]); + + psObject->dfXMax = MAX(psObject->dfXMax, psObject->padfX[i]); + psObject->dfYMax = MAX(psObject->dfYMax, psObject->padfY[i]); + psObject->dfZMax = MAX(psObject->dfZMax, psObject->padfZ[i]); + psObject->dfMMax = MAX(psObject->dfMMax, psObject->padfM[i]); + } +} + +/************************************************************************/ +/* SHPCreateObject() */ +/* */ +/* Create a shape object. It should be freed with */ +/* SHPDestroyObject(). */ +/************************************************************************/ + +SHPObject *SHPCreateObject( int nSHPType, int nShapeId, int nParts, + int * panPartStart, int * panPartType, + int nVertices, double * padfX, double * padfY, + double * padfZ, double * padfM ) + +{ + SHPObject *psObject; + int i, bHasM, bHasZ; + + psObject = (SHPObject *) calloc(1,sizeof(SHPObject)); + psObject->nSHPType = nSHPType; + psObject->nShapeId = nShapeId; + +/* -------------------------------------------------------------------- */ +/* Establish whether this shape type has M, and Z values. */ +/* -------------------------------------------------------------------- */ + if( nSHPType == SHPT_ARCM + || nSHPType == SHPT_POINTM + || nSHPType == SHPT_POLYGONM + || nSHPType == SHPT_MULTIPOINTM ) + { + bHasM = TRUE; + bHasZ = FALSE; + } + else if( nSHPType == SHPT_ARCZ + || nSHPType == SHPT_POINTZ + || nSHPType == SHPT_POLYGONZ + || nSHPType == SHPT_MULTIPOINTZ + || nSHPType == SHPT_MULTIPATCH ) + { + bHasM = TRUE; + bHasZ = TRUE; + } + else + { + bHasM = FALSE; + bHasZ = FALSE; + } + +/* -------------------------------------------------------------------- */ +/* Capture parts. Note that part type is optional, and */ +/* defaults to ring. */ +/* -------------------------------------------------------------------- */ + if( nSHPType == SHPT_ARC || nSHPType == SHPT_POLYGON + || nSHPType == SHPT_ARCM || nSHPType == SHPT_POLYGONM + || nSHPType == SHPT_ARCZ || nSHPType == SHPT_POLYGONZ + || nSHPType == SHPT_MULTIPATCH ) + { + psObject->nParts = MAX(1,nParts); + + psObject->panPartStart = (int *) + malloc(sizeof(int) * psObject->nParts); + psObject->panPartType = (int *) + malloc(sizeof(int) * psObject->nParts); + + psObject->panPartStart[0] = 0; + psObject->panPartType[0] = SHPP_RING; + + for( i = 0; i < nParts; i++ ) + { + psObject->panPartStart[i] = panPartStart[i]; + if( panPartType != NULL ) + psObject->panPartType[i] = panPartType[i]; + else + psObject->panPartType[i] = SHPP_RING; + } + } + +/* -------------------------------------------------------------------- */ +/* Capture vertices. Note that Z and M are optional, but X and */ +/* Y are not. */ +/* -------------------------------------------------------------------- */ + psObject->padfX = (double *) calloc(sizeof(double),nVertices); + psObject->padfY = (double *) calloc(sizeof(double),nVertices); + psObject->padfZ = (double *) calloc(sizeof(double),nVertices); + psObject->padfM = (double *) calloc(sizeof(double),nVertices); + + assert( padfX != NULL ); + assert( padfY != NULL ); + + for( i = 0; i < nVertices; i++ ) + { + psObject->padfX[i] = padfX[i]; + psObject->padfY[i] = padfY[i]; + if( padfZ != NULL && bHasZ ) + psObject->padfZ[i] = padfZ[i]; + if( padfM != NULL && bHasM ) + psObject->padfM[i] = padfM[i]; + } + +/* -------------------------------------------------------------------- */ +/* Compute the extents. */ +/* -------------------------------------------------------------------- */ + psObject->nVertices = nVertices; + + SHPComputeExtents( psObject ); + + return( psObject ); +} + +/************************************************************************/ +/* SHPCreateSimpleObject() */ +/* */ +/* Create a simple (common) shape object. Destroy with */ +/* SHPDestroyObject(). */ +/************************************************************************/ + +SHPObject *SHPCreateSimpleObject( int nSHPType, int nVertices, + double * padfX, double * padfY, + double * padfZ ) + +{ + return( SHPCreateObject( nSHPType, -1, 0, NULL, NULL, + nVertices, padfX, padfY, padfZ, NULL ) ); +} + +/************************************************************************/ +/* SHPWriteObject() */ +/* */ +/* Write out the vertices of a new structure. Note that it is */ +/* only possible to write vertices at the end of the file. */ +/************************************************************************/ + +int SHPWriteObject(SHPHandle psSHP, int nShapeId, SHPObject * psObject ) + +{ + int nRecordOffset, i, nRecordSize; + uchar *pabyRec; + int32 i32; + + psSHP->bUpdated = TRUE; + +/* -------------------------------------------------------------------- */ +/* Ensure that shape object matches the type of the file it is */ +/* being written to. */ +/* -------------------------------------------------------------------- */ + assert( psObject->nSHPType == psSHP->nShapeType ); + +/* -------------------------------------------------------------------- */ +/* Add the new entity to the in memory index. */ +/* -------------------------------------------------------------------- */ + if( nShapeId == -1 && psSHP->nRecords+1 > psSHP->nMaxRecords ) + { + psSHP->nMaxRecords =(int) ( psSHP->nMaxRecords * 1.3 + 100); + + psSHP->panRecOffset = (int *) + SfRealloc(psSHP->panRecOffset,sizeof(int) * psSHP->nMaxRecords ); + psSHP->panRecSize = (int *) + SfRealloc(psSHP->panRecSize,sizeof(int) * psSHP->nMaxRecords ); + } + +/* -------------------------------------------------------------------- */ +/* Initialize record. */ +/* -------------------------------------------------------------------- */ + pabyRec = (uchar *) malloc(psObject->nVertices * 4 * sizeof(double) + + psObject->nParts * 8 + 128); + +/* -------------------------------------------------------------------- */ +/* Extract vertices for a Polygon or Arc. */ +/* -------------------------------------------------------------------- */ + if( psSHP->nShapeType == SHPT_POLYGON + || psSHP->nShapeType == SHPT_POLYGONZ + || psSHP->nShapeType == SHPT_POLYGONM + || psSHP->nShapeType == SHPT_ARC + || psSHP->nShapeType == SHPT_ARCZ + || psSHP->nShapeType == SHPT_ARCM + || psSHP->nShapeType == SHPT_MULTIPATCH ) + { + int32 nPoints, nParts; + int i; + + nPoints = psObject->nVertices; + nParts = psObject->nParts; + + _SHPSetBounds( pabyRec + 12, psObject ); + + if( bBigEndian ) SwapWord( 4, &nPoints ); + if( bBigEndian ) SwapWord( 4, &nParts ); + + ByteCopy( &nPoints, pabyRec + 40 + 8, 4 ); + ByteCopy( &nParts, pabyRec + 36 + 8, 4 ); + + nRecordSize = 52; + + /* + * Write part start positions. + */ + ByteCopy( psObject->panPartStart, pabyRec + 44 + 8, + 4 * psObject->nParts ); + for( i = 0; i < psObject->nParts; i++ ) + { + if( bBigEndian ) SwapWord( 4, pabyRec + 44 + 8 + 4*i ); + nRecordSize += 4; + } + + /* + * Write multipatch part types if needed. + */ + if( psSHP->nShapeType == SHPT_MULTIPATCH ) + { + memcpy( pabyRec + nRecordSize, psObject->panPartType, + 4*psObject->nParts ); + for( i = 0; i < psObject->nParts; i++ ) + { + if( bBigEndian ) SwapWord( 4, pabyRec + nRecordSize ); + nRecordSize += 4; + } + } + + /* + * Write the (x,y) vertex values. + */ + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfX + i, pabyRec + nRecordSize, 8 ); + ByteCopy( psObject->padfY + i, pabyRec + nRecordSize + 8, 8 ); + + if( bBigEndian ) + SwapWord( 8, pabyRec + nRecordSize ); + + if( bBigEndian ) + SwapWord( 8, pabyRec + nRecordSize + 8 ); + + nRecordSize += 2 * 8; + } + + /* + * Write the Z coordinates (if any). + */ + if( psSHP->nShapeType == SHPT_POLYGONZ + || psSHP->nShapeType == SHPT_ARCZ + || psSHP->nShapeType == SHPT_MULTIPATCH ) + { + ByteCopy( &(psObject->dfZMin), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + ByteCopy( &(psObject->dfZMax), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfZ + i, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + } + + /* + * Write the M values, if any. + */ + if( psSHP->nShapeType == SHPT_POLYGONM + || psSHP->nShapeType == SHPT_ARCM +#ifndef DISABLE_MULTIPATCH_MEASURE + || psSHP->nShapeType == SHPT_MULTIPATCH +#endif + || psSHP->nShapeType == SHPT_POLYGONZ + || psSHP->nShapeType == SHPT_ARCZ ) + { + ByteCopy( &(psObject->dfMMin), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + ByteCopy( &(psObject->dfMMax), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfM + i, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + } + } + +/* -------------------------------------------------------------------- */ +/* Extract vertices for a MultiPoint. */ +/* -------------------------------------------------------------------- */ + else if( psSHP->nShapeType == SHPT_MULTIPOINT + || psSHP->nShapeType == SHPT_MULTIPOINTZ + || psSHP->nShapeType == SHPT_MULTIPOINTM ) + { + int32 nPoints; + int i; + + nPoints = psObject->nVertices; + + _SHPSetBounds( pabyRec + 12, psObject ); + + if( bBigEndian ) SwapWord( 4, &nPoints ); + ByteCopy( &nPoints, pabyRec + 44, 4 ); + + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfX + i, pabyRec + 48 + i*16, 8 ); + ByteCopy( psObject->padfY + i, pabyRec + 48 + i*16 + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, pabyRec + 48 + i*16 ); + if( bBigEndian ) SwapWord( 8, pabyRec + 48 + i*16 + 8 ); + } + + nRecordSize = 48 + 16 * psObject->nVertices; + + if( psSHP->nShapeType == SHPT_MULTIPOINTZ ) + { + ByteCopy( &(psObject->dfZMin), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + ByteCopy( &(psObject->dfZMax), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfZ + i, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + } + + if( psSHP->nShapeType == SHPT_MULTIPOINTZ + || psSHP->nShapeType == SHPT_MULTIPOINTM ) + { + ByteCopy( &(psObject->dfMMin), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + ByteCopy( &(psObject->dfMMax), pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + + for( i = 0; i < psObject->nVertices; i++ ) + { + ByteCopy( psObject->padfM + i, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + } + } + +/* -------------------------------------------------------------------- */ +/* Extract vertices for a point. */ +/* -------------------------------------------------------------------- */ + else if( psSHP->nShapeType == SHPT_POINT + || psSHP->nShapeType == SHPT_POINTZ + || psSHP->nShapeType == SHPT_POINTM ) + { + ByteCopy( psObject->padfX, pabyRec + 12, 8 ); + ByteCopy( psObject->padfY, pabyRec + 20, 8 ); + + if( bBigEndian ) SwapWord( 8, pabyRec + 12 ); + if( bBigEndian ) SwapWord( 8, pabyRec + 20 ); + + nRecordSize = 28; + + if( psSHP->nShapeType == SHPT_POINTZ ) + { + ByteCopy( psObject->padfZ, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + + if( psSHP->nShapeType == SHPT_POINTZ + || psSHP->nShapeType == SHPT_POINTM ) + { + ByteCopy( psObject->padfM, pabyRec + nRecordSize, 8 ); + if( bBigEndian ) SwapWord( 8, pabyRec + nRecordSize ); + nRecordSize += 8; + } + } + + else + { + /* unknown type */ + assert( FALSE ); + } + +/* -------------------------------------------------------------------- */ +/* Establish where we are going to put this record. If we are */ +/* rewriting and existing record, and it will fit, then put it */ +/* back where the original came from. Otherwise write at the end. */ +/* -------------------------------------------------------------------- */ + if( nShapeId == -1 || psSHP->panRecSize[nShapeId] < nRecordSize-8 ) + { + if( nShapeId == -1 ) + nShapeId = psSHP->nRecords++; + + psSHP->panRecOffset[nShapeId] = nRecordOffset = psSHP->nFileSize; + psSHP->panRecSize[nShapeId] = nRecordSize-8; + psSHP->nFileSize += nRecordSize; + } + else + { + nRecordOffset = psSHP->panRecOffset[nShapeId]; + } + +/* -------------------------------------------------------------------- */ +/* Set the shape type, record number, and record size. */ +/* -------------------------------------------------------------------- */ + i32 = nShapeId+1; /* record # */ + if( !bBigEndian ) SwapWord( 4, &i32 ); + ByteCopy( &i32, pabyRec, 4 ); + + i32 = (nRecordSize-8)/2; /* record size */ + if( !bBigEndian ) SwapWord( 4, &i32 ); + ByteCopy( &i32, pabyRec + 4, 4 ); + + i32 = psSHP->nShapeType; /* shape type */ + if( bBigEndian ) SwapWord( 4, &i32 ); + ByteCopy( &i32, pabyRec + 8, 4 ); + +/* -------------------------------------------------------------------- */ +/* Write out record. */ +/* -------------------------------------------------------------------- */ + if( fseek( psSHP->fpSHP, nRecordOffset, 0 ) != 0 + || fwrite( pabyRec, nRecordSize, 1, psSHP->fpSHP ) < 1 ) + { + printf( "Error in fseek() or fwrite().\n" ); + free( pabyRec ); + return -1; + } + + free( pabyRec ); + +/* -------------------------------------------------------------------- */ +/* Expand file wide bounds based on this shape. */ +/* -------------------------------------------------------------------- */ + if( psSHP->nRecords == 1 ) + { + psSHP->adBoundsMin[0] = psSHP->adBoundsMax[0] = psObject->padfX[0]; + psSHP->adBoundsMin[1] = psSHP->adBoundsMax[1] = psObject->padfY[0]; + psSHP->adBoundsMin[2] = psSHP->adBoundsMax[2] = psObject->padfZ[0]; + psSHP->adBoundsMin[3] = psSHP->adBoundsMax[3] = psObject->padfM[0]; + } + + for( i = 0; i < psObject->nVertices; i++ ) + { + psSHP->adBoundsMin[0] = MIN(psSHP->adBoundsMin[0],psObject->padfX[i]); + psSHP->adBoundsMin[1] = MIN(psSHP->adBoundsMin[1],psObject->padfY[i]); + psSHP->adBoundsMin[2] = MIN(psSHP->adBoundsMin[2],psObject->padfZ[i]); + psSHP->adBoundsMin[3] = MIN(psSHP->adBoundsMin[3],psObject->padfM[i]); + psSHP->adBoundsMax[0] = MAX(psSHP->adBoundsMax[0],psObject->padfX[i]); + psSHP->adBoundsMax[1] = MAX(psSHP->adBoundsMax[1],psObject->padfY[i]); + psSHP->adBoundsMax[2] = MAX(psSHP->adBoundsMax[2],psObject->padfZ[i]); + psSHP->adBoundsMax[3] = MAX(psSHP->adBoundsMax[3],psObject->padfM[i]); + } + + return( nShapeId ); +} + +/************************************************************************/ +/* SHPReadObject() */ +/* */ +/* Read the vertices, parts, and other non-attribute information */ +/* for one shape. */ +/************************************************************************/ + +SHPObject *SHPReadObject( SHPHandle psSHP, int hEntity ) + +{ + SHPObject *psShape; + +/* -------------------------------------------------------------------- */ +/* Validate the record/entity number. */ +/* -------------------------------------------------------------------- */ + if( hEntity < 0 || hEntity >= psSHP->nRecords ) + return( NULL ); + +/* -------------------------------------------------------------------- */ +/* Ensure our record buffer is large enough. */ +/* -------------------------------------------------------------------- */ + if( psSHP->panRecSize[hEntity]+8 > nBufSize ) + { + nBufSize = psSHP->panRecSize[hEntity]+8; + pabyRec = (uchar *) SfRealloc(pabyRec,nBufSize); + } + +/* -------------------------------------------------------------------- */ +/* Read the record. */ +/* -------------------------------------------------------------------- */ + fseek( psSHP->fpSHP, psSHP->panRecOffset[hEntity], 0 ); + fread( pabyRec, psSHP->panRecSize[hEntity]+8, 1, psSHP->fpSHP ); + +/* -------------------------------------------------------------------- */ +/* Allocate and minimally initialize the object. */ +/* -------------------------------------------------------------------- */ + psShape = (SHPObject *) calloc(1,sizeof(SHPObject)); + psShape->nShapeId = hEntity; + + memcpy( &psShape->nSHPType, pabyRec + 8, 4 ); + if( bBigEndian ) SwapWord( 4, &(psShape->nSHPType) ); + +/* ==================================================================== */ +/* Extract vertices for a Polygon or Arc. */ +/* ==================================================================== */ + if( psShape->nSHPType == SHPT_POLYGON || psShape->nSHPType == SHPT_ARC + || psShape->nSHPType == SHPT_POLYGONZ + || psShape->nSHPType == SHPT_POLYGONM + || psShape->nSHPType == SHPT_ARCZ + || psShape->nSHPType == SHPT_ARCM + || psShape->nSHPType == SHPT_MULTIPATCH ) + { + int32 nPoints, nParts; + int i, nOffset; + +/* -------------------------------------------------------------------- */ +/* Get the X/Y bounds. */ +/* -------------------------------------------------------------------- */ + memcpy( &(psShape->dfXMin), pabyRec + 8 + 4, 8 ); + memcpy( &(psShape->dfYMin), pabyRec + 8 + 12, 8 ); + memcpy( &(psShape->dfXMax), pabyRec + 8 + 20, 8 ); + memcpy( &(psShape->dfYMax), pabyRec + 8 + 28, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfXMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfYMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfXMax) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfYMax) ); + +/* -------------------------------------------------------------------- */ +/* Extract part/point count, and build vertex and part arrays */ +/* to proper size. */ +/* -------------------------------------------------------------------- */ + memcpy( &nPoints, pabyRec + 40 + 8, 4 ); + memcpy( &nParts, pabyRec + 36 + 8, 4 ); + + if( bBigEndian ) SwapWord( 4, &nPoints ); + if( bBigEndian ) SwapWord( 4, &nParts ); + + psShape->nVertices = nPoints; + psShape->padfX = (double *) calloc(nPoints,sizeof(double)); + psShape->padfY = (double *) calloc(nPoints,sizeof(double)); + psShape->padfZ = (double *) calloc(nPoints,sizeof(double)); + psShape->padfM = (double *) calloc(nPoints,sizeof(double)); + + psShape->nParts = nParts; + psShape->panPartStart = (int *) calloc(nParts,sizeof(int)); + psShape->panPartType = (int *) calloc(nParts,sizeof(int)); + + for( i = 0; i < nParts; i++ ) + psShape->panPartType[i] = SHPP_RING; + +/* -------------------------------------------------------------------- */ +/* Copy out the part array from the record. */ +/* -------------------------------------------------------------------- */ + memcpy( psShape->panPartStart, pabyRec + 44 + 8, 4 * nParts ); + for( i = 0; i < nParts; i++ ) + { + if( bBigEndian ) SwapWord( 4, psShape->panPartStart+i ); + } + + nOffset = 44 + 8 + 4*nParts; + +/* -------------------------------------------------------------------- */ +/* If this is a multipatch, we will also have parts types. */ +/* -------------------------------------------------------------------- */ + if( psShape->nSHPType == SHPT_MULTIPATCH ) + { + memcpy( psShape->panPartType, pabyRec + nOffset, 4*nParts ); + for( i = 0; i < nParts; i++ ) + { + if( bBigEndian ) SwapWord( 4, psShape->panPartType+i ); + } + + nOffset += 4*nParts; + } + +/* -------------------------------------------------------------------- */ +/* Copy out the vertices from the record. */ +/* -------------------------------------------------------------------- */ + for( i = 0; i < nPoints; i++ ) + { + memcpy(psShape->padfX + i, + pabyRec + nOffset + i * 16, + 8 ); + + memcpy(psShape->padfY + i, + pabyRec + nOffset + i * 16 + 8, + 8 ); + + if( bBigEndian ) SwapWord( 8, psShape->padfX + i ); + if( bBigEndian ) SwapWord( 8, psShape->padfY + i ); + } + + nOffset += 16*nPoints; + +/* -------------------------------------------------------------------- */ +/* If we have a Z coordinate, collect that now. */ +/* -------------------------------------------------------------------- */ + if( psShape->nSHPType == SHPT_POLYGONZ + || psShape->nSHPType == SHPT_ARCZ + || psShape->nSHPType == SHPT_MULTIPATCH ) + { + memcpy( &(psShape->dfZMin), pabyRec + nOffset, 8 ); + memcpy( &(psShape->dfZMax), pabyRec + nOffset + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfZMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfZMax) ); + + for( i = 0; i < nPoints; i++ ) + { + memcpy( psShape->padfZ + i, + pabyRec + nOffset + 16 + i*8, 8 ); + if( bBigEndian ) SwapWord( 8, psShape->padfZ + i ); + } + + nOffset += 16 + 8*nPoints; + } + +/* -------------------------------------------------------------------- */ +/* If we have a M measure value, then read it now. We assume */ +/* that the measure can be present for any shape if the size is */ +/* big enough, but really it will only occur for the Z shapes */ +/* (options), and the M shapes. */ +/* -------------------------------------------------------------------- */ + if( psSHP->panRecSize[hEntity]+8 >= nOffset + 16 + 8*nPoints ) + { + memcpy( &(psShape->dfMMin), pabyRec + nOffset, 8 ); + memcpy( &(psShape->dfMMax), pabyRec + nOffset + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfMMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfMMax) ); + + for( i = 0; i < nPoints; i++ ) + { + memcpy( psShape->padfM + i, + pabyRec + nOffset + 16 + i*8, 8 ); + if( bBigEndian ) SwapWord( 8, psShape->padfM + i ); + } + } + + } + +/* ==================================================================== */ +/* Extract vertices for a MultiPoint. */ +/* ==================================================================== */ + else if( psShape->nSHPType == SHPT_MULTIPOINT + || psShape->nSHPType == SHPT_MULTIPOINTM + || psShape->nSHPType == SHPT_MULTIPOINTZ ) + { + int32 nPoints; + int i, nOffset; + + memcpy( &nPoints, pabyRec + 44, 4 ); + if( bBigEndian ) SwapWord( 4, &nPoints ); + + psShape->nVertices = nPoints; + psShape->padfX = (double *) calloc(nPoints,sizeof(double)); + psShape->padfY = (double *) calloc(nPoints,sizeof(double)); + psShape->padfZ = (double *) calloc(nPoints,sizeof(double)); + psShape->padfM = (double *) calloc(nPoints,sizeof(double)); + + for( i = 0; i < nPoints; i++ ) + { + memcpy(psShape->padfX+i, pabyRec + 48 + 16 * i, 8 ); + memcpy(psShape->padfY+i, pabyRec + 48 + 16 * i + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, psShape->padfX + i ); + if( bBigEndian ) SwapWord( 8, psShape->padfY + i ); + } + + nOffset = 48 + 16*nPoints; + +/* -------------------------------------------------------------------- */ +/* Get the X/Y bounds. */ +/* -------------------------------------------------------------------- */ + memcpy( &(psShape->dfXMin), pabyRec + 8 + 4, 8 ); + memcpy( &(psShape->dfYMin), pabyRec + 8 + 12, 8 ); + memcpy( &(psShape->dfXMax), pabyRec + 8 + 20, 8 ); + memcpy( &(psShape->dfYMax), pabyRec + 8 + 28, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfXMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfYMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfXMax) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfYMax) ); + +/* -------------------------------------------------------------------- */ +/* If we have a Z coordinate, collect that now. */ +/* -------------------------------------------------------------------- */ + if( psShape->nSHPType == SHPT_MULTIPOINTZ ) + { + memcpy( &(psShape->dfZMin), pabyRec + nOffset, 8 ); + memcpy( &(psShape->dfZMax), pabyRec + nOffset + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfZMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfZMax) ); + + for( i = 0; i < nPoints; i++ ) + { + memcpy( psShape->padfZ + i, + pabyRec + nOffset + 16 + i*8, 8 ); + if( bBigEndian ) SwapWord( 8, psShape->padfZ + i ); + } + + nOffset += 16 + 8*nPoints; + } + +/* -------------------------------------------------------------------- */ +/* If we have a M measure value, then read it now. We assume */ +/* that the measure can be present for any shape if the size is */ +/* big enough, but really it will only occur for the Z shapes */ +/* (options), and the M shapes. */ +/* -------------------------------------------------------------------- */ + if( psSHP->panRecSize[hEntity]+8 >= nOffset + 16 + 8*nPoints ) + { + memcpy( &(psShape->dfMMin), pabyRec + nOffset, 8 ); + memcpy( &(psShape->dfMMax), pabyRec + nOffset + 8, 8 ); + + if( bBigEndian ) SwapWord( 8, &(psShape->dfMMin) ); + if( bBigEndian ) SwapWord( 8, &(psShape->dfMMax) ); + + for( i = 0; i < nPoints; i++ ) + { + memcpy( psShape->padfM + i, + pabyRec + nOffset + 16 + i*8, 8 ); + if( bBigEndian ) SwapWord( 8, psShape->padfM + i ); + } + } + } + +/* ==================================================================== */ +/* Extract vertices for a point. */ +/* ==================================================================== */ + else if( psShape->nSHPType == SHPT_POINT + || psShape->nSHPType == SHPT_POINTM + || psShape->nSHPType == SHPT_POINTZ ) + { + int nOffset; + + psShape->nVertices = 1; + psShape->padfX = (double *) calloc(1,sizeof(double)); + psShape->padfY = (double *) calloc(1,sizeof(double)); + psShape->padfZ = (double *) calloc(1,sizeof(double)); + psShape->padfM = (double *) calloc(1,sizeof(double)); + + memcpy( psShape->padfX, pabyRec + 12, 8 ); + memcpy( psShape->padfY, pabyRec + 20, 8 ); + + if( bBigEndian ) SwapWord( 8, psShape->padfX ); + if( bBigEndian ) SwapWord( 8, psShape->padfY ); + + nOffset = 20 + 8; + +/* -------------------------------------------------------------------- */ +/* If we have a Z coordinate, collect that now. */ +/* -------------------------------------------------------------------- */ + if( psShape->nSHPType == SHPT_POINTZ ) + { + memcpy( psShape->padfZ, pabyRec + nOffset, 8 ); + + if( bBigEndian ) SwapWord( 8, psShape->padfZ ); + + nOffset += 8; + } + +/* -------------------------------------------------------------------- */ +/* If we have a M measure value, then read it now. We assume */ +/* that the measure can be present for any shape if the size is */ +/* big enough, but really it will only occur for the Z shapes */ +/* (options), and the M shapes. */ +/* -------------------------------------------------------------------- */ + if( psSHP->panRecSize[hEntity]+8 >= nOffset + 8 ) + { + memcpy( psShape->padfM, pabyRec + nOffset, 8 ); + + if( bBigEndian ) SwapWord( 8, psShape->padfM ); + } + +/* -------------------------------------------------------------------- */ +/* Since no extents are supplied in the record, we will apply */ +/* them from the single vertex. */ +/* -------------------------------------------------------------------- */ + psShape->dfXMin = psShape->dfXMax = psShape->padfX[0]; + psShape->dfYMin = psShape->dfYMax = psShape->padfY[0]; + psShape->dfZMin = psShape->dfZMax = psShape->padfZ[0]; + psShape->dfMMin = psShape->dfMMax = psShape->padfM[0]; + } + + return( psShape ); +} + +/************************************************************************/ +/* SHPTypeName() */ +/************************************************************************/ + +const char *SHPTypeName( int nSHPType ) + +{ + switch( nSHPType ) + { + case 0: + return "NullShape"; + + case SHPT_POINT: + return "Point"; + + case SHPT_ARC: + return "Arc"; + + case SHPT_POLYGON: + return "Polygon"; + + case SHPT_MULTIPOINT: + return "MultiPoint"; + + case SHPT_POINTZ: + return "PointZ"; + + case SHPT_ARCZ: + return "ArcZ"; + + case SHPT_POLYGONZ: + return "PolygonZ"; + + case SHPT_MULTIPOINTZ: + return "MultiPointZ"; + + case SHPT_POINTM: + return "PointM"; + + case SHPT_ARCM: + return "ArcM"; + + case SHPT_POLYGONM: + return "PolygonM"; + + case SHPT_MULTIPOINTM: + return "MultiPointM"; + + case SHPT_MULTIPATCH: + return "MultiPatch"; + + default: + return "UnknownShapeType"; + } +} + +/************************************************************************/ +/* SHPPartTypeName() */ +/************************************************************************/ + +const char *SHPPartTypeName( int nPartType ) + +{ + switch( nPartType ) + { + case SHPP_TRISTRIP: + return "TriangleStrip"; + + case SHPP_TRIFAN: + return "TriangleFan"; + + case SHPP_OUTERRING: + return "OuterRing"; + + case SHPP_INNERRING: + return "InnerRing"; + + case SHPP_FIRSTRING: + return "FirstRing"; + + case SHPP_RING: + return "Ring"; + + default: + return "UnknownPartType"; + } +} + +/************************************************************************/ +/* SHPDestroyObject() */ +/************************************************************************/ + +void SHPDestroyObject( SHPObject * psShape ) + +{ + if( psShape == NULL ) + return; + + if( psShape->padfX != NULL ) + free( psShape->padfX ); + if( psShape->padfY != NULL ) + free( psShape->padfY ); + if( psShape->padfZ != NULL ) + free( psShape->padfZ ); + if( psShape->padfM != NULL ) + free( psShape->padfM ); + + if( psShape->panPartStart != NULL ) + free( psShape->panPartStart ); + if( psShape->panPartType != NULL ) + free( psShape->panPartType ); + + free( psShape ); +} diff --git a/postgis.h b/postgis.h index c54aee45a..c034e77b9 100644 --- a/postgis.h +++ b/postgis.h @@ -20,6 +20,22 @@ #define COLLECTIONTYPE 7 #define BBOXONLYTYPE 99 + +//standard definition of an ellipsoid (what wkt calls a spheroid) +// f = (a-b)/a +// e_sq = (a*a - b*b)/(a*a) +// b = a - fa +typedef struct +{ + double a; //semimajor axis + double b; //semiminor axis + double f; //flattening + double e; //eccentricity (first) + double e_sq; //eccentricity (first), squared + char name[20]; //name of ellipse +} SPHEROID; + + /*--------------------------------------------------------------------- * POINT3D - (x,y,z) * Base type for all geometries @@ -27,7 +43,7 @@ *-------------------------------------------------------------------*/ typedef struct { - double x,y,z; + double x,y,z; //for lat/long x=long, y=lat } POINT3D; /*--------------------------------------------------------------------- @@ -189,6 +205,10 @@ typedef struct geomkey { int isspace(int c); +/* constructors*/ +POLYGON3D *make_polygon(int nrings, int *pts_per_ring, POINT3D *pts, int npoints, int *size); +void set_point( POINT3D *pt,double x, double y, double z); +GEOMETRY *make_oneobj_geometry(int sub_obj_size, char *sub_obj, int type, bool is3d); void print_box(BOX3D *box); void print_box_oneline(BOX3D *box); @@ -241,7 +261,10 @@ bool linestring_inside_box(POINT3D *pts, int npoints, BOX3D *box); bool polygon_truely_inside(POLYGON3D *poly, BOX3D *box); bool line_truely_inside( LINE3D *line, BOX3D *box); void translate_points(POINT3D *pt, int npoints,double x_off, double y_off, double z_off); - +int size_subobject (char *sub_obj, int type); +GEOMETRY *add_to_geometry(GEOMETRY *geom,int sub_obj_size, char *sub_obj, int type); +LINE3D *make_line(int npoints, POINT3D *pts, int *size); +char *print_geometry(GEOMETRY *geom); void swap_char(char *a, char*b); void flip_endian_double(char *dd); @@ -262,6 +285,18 @@ void decode_wkb_collection(char *wkb,int *size); void decode_wkb(char *wkb, int *size); void dump_bytes( char *a, int numb); +double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere); +double bigA(double u2); +double bigB(double u2); +double distance_ellipse(double lat1, double long1, + double lat2, double long2, + SPHEROID *sphere); +double mu2(double azimuth,SPHEROID *sphere); + +double length2d_ellipse_linestring(LINE3D *line, SPHEROID *sphere); +double length3d_ellipse_linestring(LINE3D *line, SPHEROID *sphere); + + //exposed to psql Datum box3d_in(PG_FUNCTION_ARGS); @@ -304,14 +339,40 @@ Datum numb_sub_objs(PG_FUNCTION_ARGS); Datum summary(PG_FUNCTION_ARGS); Datum translate(PG_FUNCTION_ARGS); -Datum wkb_XDR(PG_FUNCTION_ARGS); -Datum wkb_NDR(PG_FUNCTION_ARGS); +Datum asbinary_specify(PG_FUNCTION_ARGS); +Datum asbinary_simple(PG_FUNCTION_ARGS); Datum force_2d(PG_FUNCTION_ARGS); Datum force_3d(PG_FUNCTION_ARGS); +Datum force_collection(PG_FUNCTION_ARGS); Datum combine_bbox(PG_FUNCTION_ARGS); +Datum dimension(PG_FUNCTION_ARGS); +Datum geometrytype(PG_FUNCTION_ARGS); +Datum envelope(PG_FUNCTION_ARGS); + +Datum x_point(PG_FUNCTION_ARGS); +Datum y_point(PG_FUNCTION_ARGS); +Datum z_point(PG_FUNCTION_ARGS); + +Datum numpoints_linestring(PG_FUNCTION_ARGS); +Datum pointn_linestring(PG_FUNCTION_ARGS); + +Datum exteriorring_polygon(PG_FUNCTION_ARGS); +Datum numinteriorrings_polygon(PG_FUNCTION_ARGS); +Datum interiorringn_polygon(PG_FUNCTION_ARGS); + +Datum numgeometries_collection(PG_FUNCTION_ARGS); +Datum geometryn_collection(PG_FUNCTION_ARGS); + +Datum ellipsoid_out(PG_FUNCTION_ARGS); +Datum ellipsoid_in(PG_FUNCTION_ARGS); +Datum length_ellipsoid(PG_FUNCTION_ARGS); +Datum length3d_ellipsoid(PG_FUNCTION_ARGS); + +Datum point_inside_circle(PG_FUNCTION_ARGS); + //for GIST index diff --git a/postgis.sql.in b/postgis.sql.in index c2b0553a6..5b797d812 100644 --- a/postgis.sql.in +++ b/postgis.sql.in @@ -5,18 +5,35 @@ BEGIN TRANSACTION; CREATE FUNCTION BOX3D_in(opaque) RETURNS BOX3D AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION BOX3D_out(opaque) RETURNS opaque AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION SPHEROID_in(opaque) + RETURNS SPHEROID + AS '@MODULE_FILENAME@','ellipsoid_in' + LANGUAGE 'c' with (isstrict,iscachable); + +CREATE FUNCTION SPHEROID_out(opaque) + RETURNS opaque + AS '@MODULE_FILENAME@','ellipsoid_out' + LANGUAGE 'c' with (isstrict); + +CREATE TYPE SPHEROID ( + alignment = double, + internallength = 65, + input = SPHEROID_in, + output = SPHEROID_out +); CREATE TYPE BOX3D ( alignment = double, internallength = 48, - input = BOX3D_in, - output = BOX3D_out + input = BOX3D_in, + output = BOX3D_out ); CREATE TYPE WKB ( @@ -31,12 +48,12 @@ CREATE TYPE WKB ( create function geometry_in(opaque) RETURNS GEOMETRY AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); create function geometry_out(opaque) RETURNS opaque AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE TYPE GEOMETRY ( alignment = double, @@ -50,24 +67,24 @@ CREATE TYPE GEOMETRY ( CREATE FUNCTION box3d(GEOMETRY) RETURNS BOX3D AS '@MODULE_FILENAME@','get_bbox_of_geometry' - LANGUAGE 'c' WITH (iscachable); + LANGUAGE 'c' WITH (iscachable,isstrict); CREATE FUNCTION geometry(BOX3D) RETURNS GEOMETRY AS '@MODULE_FILENAME@','get_geometry_of_bbox' - LANGUAGE 'c' WITH (iscachable); + LANGUAGE 'c' WITH (iscachable,isstrict); --------- functions for converting to wkb -CREATE FUNCTION wkb_NDR(GEOMETRY) +CREATE FUNCTION asbinary(GEOMETRY) RETURNS WKB - AS '@MODULE_FILENAME@','wkb_NDR' - LANGUAGE 'c' WITH (iscachable); + AS '@MODULE_FILENAME@','asbinary_simple' + LANGUAGE 'c' WITH (iscachable,isstrict); -CREATE FUNCTION wkb_XDR(GEOMETRY) +CREATE FUNCTION asbinary(GEOMETRY,TEXT) RETURNS WKB - AS '@MODULE_FILENAME@','wkb_XDR' - LANGUAGE 'c' WITH (iscachable); + AS '@MODULE_FILENAME@','asbinary_specify' + LANGUAGE 'c' WITH (iscachable,isstrict); ---- Debug (info) functions @@ -75,32 +92,110 @@ CREATE FUNCTION wkb_XDR(GEOMETRY) CREATE FUNCTION npoints(GEOMETRY) RETURNS INT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION nrings(GEOMETRY) RETURNS INT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict) ; CREATE FUNCTION mem_size(GEOMETRY) RETURNS INT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION numb_sub_objs(GEOMETRY) RETURNS INT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION summary(GEOMETRY) RETURNS text AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION translate(GEOMETRY,float8,float8,float8) RETURNS GEOMETRY AS '@MODULE_FILENAME@' - LANGUAGE 'c' ; + LANGUAGE 'c' with (isstrict) ; + +CREATE FUNCTION dimension(GEOMETRY) + RETURNS INT4 + AS '@MODULE_FILENAME@' + LANGUAGE 'c' with (isstrict) ; + +CREATE FUNCTION geometrytype(GEOMETRY) + RETURNS text + AS '@MODULE_FILENAME@' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION envelope(GEOMETRY) + RETURNS geometry + AS '@MODULE_FILENAME@' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION x(GEOMETRY) + RETURNS float4 + AS '@MODULE_FILENAME@','x_point' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION y(GEOMETRY) + RETURNS float4 + AS '@MODULE_FILENAME@','y_point' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION z(GEOMETRY) + RETURNS float4 + AS '@MODULE_FILENAME@','z_point' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION numpoints(GEOMETRY) + RETURNS integer + AS '@MODULE_FILENAME@','numpoints_linestring' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION pointn(GEOMETRY,INTEGER) + RETURNS GEOMETRY + AS '@MODULE_FILENAME@','pointn_linestring' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION exteriorring(GEOMETRY) + RETURNS GEOMETRY + AS '@MODULE_FILENAME@','exteriorring_polygon' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION numinteriorrings(GEOMETRY) + RETURNS INTEGER + AS '@MODULE_FILENAME@','numinteriorrings_polygon' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION interiorringn(GEOMETRY,INTEGER) + RETURNS GEOMETRY + AS '@MODULE_FILENAME@','interiorringn_polygon' + LANGUAGE 'c' with (isstrict); + + +CREATE FUNCTION numgeometries(GEOMETRY) + RETURNS INTEGER + AS '@MODULE_FILENAME@','numgeometries_collection' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION geometryn(GEOMETRY,INTEGER) + RETURNS GEOMETRY + AS '@MODULE_FILENAME@','geometryn_collection' + LANGUAGE 'c' with (isstrict); + +------- spheroid calcs + +CREATE FUNCTION length_spheroid(GEOMETRY,SPHEROID) + RETURNS FLOAT4 + AS '@MODULE_FILENAME@','length_ellipsoid' + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION length3d_spheroid(GEOMETRY,SPHEROID) + RETURNS FLOAT4 + AS '@MODULE_FILENAME@','length3d_ellipsoid' + LANGUAGE 'c' with (isstrict); ------- generic operations @@ -108,32 +203,38 @@ CREATE FUNCTION translate(GEOMETRY,float8,float8,float8) CREATE FUNCTION length3d(GEOMETRY) RETURNS FLOAT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION length2d(GEOMETRY) RETURNS FLOAT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION area2d(GEOMETRY) RETURNS FLOAT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION perimeter3d(GEOMETRY) RETURNS FLOAT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION perimeter2d(GEOMETRY) RETURNS FLOAT4 AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); CREATE FUNCTION truly_inside(GEOMETRY,GEOMETRY) RETURNS bool AS '@MODULE_FILENAME@' - LANGUAGE 'c'; + LANGUAGE 'c' with (isstrict); + +CREATE FUNCTION point_inside_circle(GEOMETRY,float8,float8,float8) + RETURNS bool + AS '@MODULE_FILENAME@' + LANGUAGE 'c' with (isstrict); + ------- Aggregate @@ -153,48 +254,52 @@ CREATE AGGREGATE extent( ------- OPERATOR functions CREATE FUNCTION geometry_overleft(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_overright(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_left(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_right(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_contain(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_contained(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_overlap(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_same(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); --------- functions for doing sorting-like things (not very usefull) CREATE FUNCTION geometry_lt(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_gt(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION geometry_eq(GEOMETRY, GEOMETRY) RETURNS bool - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); --------- functions for forcing geometry to be 2d or 3d CREATE FUNCTION force_2d(GEOMETRY) RETURNS GEOMETRY - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); CREATE FUNCTION force_3d(GEOMETRY) RETURNS GEOMETRY - AS '@MODULE_FILENAME@' LANGUAGE 'c'; + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); + +-------- cause geometry to be represented as a geometry collection +CREATE FUNCTION force_collection(GEOMETRY) RETURNS GEOMETRY + AS '@MODULE_FILENAME@' LANGUAGE 'c' with (isstrict); -------- GiST support functions diff --git a/postgis_debug.c b/postgis_debug.c index 420b3f10b..d5afd2276 100644 --- a/postgis_debug.c +++ b/postgis_debug.c @@ -714,3 +714,12 @@ void decode_wkb(char *wkb, int *size) } +char *print_geometry(GEOMETRY *geom) +{ + char *text; + + text = (char *) DatumGetPointer(DirectFunctionCall1(geometry_out,PointerGetDatum(geom) ) ) ; + + printf( "%s", text); + return text; +} diff --git a/postgis_fn.c b/postgis_fn.c index 5052bc1c0..dbed0cc0c 100644 --- a/postgis_fn.c +++ b/postgis_fn.c @@ -891,9 +891,13 @@ PG_FUNCTION_INFO_V1(force_2d); Datum force_2d(PG_FUNCTION_ARGS) { GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GEOMETRY *result; - geom->is3d = FALSE; - PG_RETURN_POINTER(geom); + result = (GEOMETRY *) palloc(geom->size); + memcpy(result,geom, geom->size); + + result->is3d = FALSE; + PG_RETURN_POINTER(result); } @@ -901,9 +905,13 @@ PG_FUNCTION_INFO_V1(force_3d); Datum force_3d(PG_FUNCTION_ARGS) { GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GEOMETRY *result; + + result = (GEOMETRY *) palloc(geom->size); + memcpy(result,geom, geom->size); - geom->is3d = TRUE; - PG_RETURN_POINTER(geom); + result->is3d = TRUE; + PG_RETURN_POINTER(result); } @@ -929,12 +937,20 @@ Datum combine_bbox(PG_FUNCTION_ARGS) if (box3d_ptr == NULL) { geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1)); - PG_RETURN_POINTER(&(geom->bvol)); // combine_bbox(null, geometry) => geometry->bvol + + box = (BOX3D *) palloc(sizeof(BOX3D)); + memcpy(box, &(geom->bvol), sizeof(BOX3D) ); + + PG_RETURN_POINTER(box); // combine_bbox(null, geometry) => geometry->bvol } if (geom_ptr == NULL) { - PG_RETURN_POINTER( PG_GETARG_DATUM(0) ); // combine_bbox(BOX3D, null) => BOX3D + box = (BOX3D *) palloc(sizeof(BOX3D)); + memcpy(box, (char *) PG_GETARG_DATUM(0) , sizeof(BOX3D) ); + + + PG_RETURN_POINTER( box ); // combine_bbox(BOX3D, null) => BOX3D } //combine_bbox(BOX3D, geometry) => union(BOX3D, geometry->bvol) @@ -970,5 +986,511 @@ Datum combine_bbox(PG_FUNCTION_ARGS) } +//return 2 for 2d geometries or 3 for 3d geometries +PG_FUNCTION_INFO_V1(dimension); +Datum dimension(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + if (geom->is3d) + PG_RETURN_INT32(3); + else + PG_RETURN_INT32(2); +} + +//returns a string representation of this geometry's type +PG_FUNCTION_INFO_V1(geometrytype); +Datum geometrytype(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *text_ob = palloc(20+4); + char *result = text_ob+4; + int32 size; + + + memset(result,0,20); + + if (geom->type == POINTTYPE) + strcpy(result,"POINT"); + if (geom->type == MULTIPOINTTYPE) + strcpy(result,"MULTIPOINT"); + + if (geom->type == LINETYPE) + strcpy(result,"LINESTRING"); + if (geom->type == MULTILINETYPE) + strcpy(result,"MULTILINESTRING"); + + if (geom->type == POLYGONTYPE) + strcpy(result,"POLYGON"); + if (geom->type == MULTIPOLYGONTYPE) + strcpy(result,"MULTIPOLYGON"); + + if (geom->type == COLLECTIONTYPE) + strcpy(result,"GEOMETRYCOLLECTION"); + + if (strlen(result) == 0) + strcpy(result,"UNKNOWN"); + + size = strlen(result) +4 ; + + memcpy(text_ob, &size,4); // size of string + + + PG_RETURN_POINTER(text_ob); + +} + + +//makes a polygon of a features bvol - 1st point = LL 3rd=UR +// 2d only +// +// create new geometry of type polygon, 1 ring, 5 points + +PG_FUNCTION_INFO_V1(envelope); +Datum envelope(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GEOMETRY *result; + POLYGON3D *poly; + POINT3D pts[5]; //5 points around box + int pts_per_ring[1]; + int poly_size; + + //use LLB's z value (we're going to set is3d to false) + + set_point( &pts[0], geom->bvol.LLB.x , geom->bvol.LLB.y , geom->bvol.LLB.z ); + set_point( &pts[1], geom->bvol.URT.x , geom->bvol.LLB.y , geom->bvol.LLB.z ); + set_point( &pts[2], geom->bvol.URT.x , geom->bvol.URT.y , geom->bvol.LLB.z ); + set_point( &pts[3], geom->bvol.LLB.x , geom->bvol.URT.y , geom->bvol.LLB.z ); + set_point( &pts[4], geom->bvol.LLB.x , geom->bvol.LLB.y , geom->bvol.LLB.z ); + + pts_per_ring[0] = 5; //ring has 5 points + + //make a polygon + poly = make_polygon(1, pts_per_ring, pts, 5, &poly_size); + + result = make_oneobj_geometry(poly_size, (char *)poly, POLYGONTYPE, FALSE); + + PG_RETURN_POINTER(result); + +} + +//X(GEOMETRY) -- find the first POINT(..) in GEOMETRY, returns its X value. +//Return NULL if there is no POINT(..) in GEOMETRY +PG_FUNCTION_INFO_V1(x_point); +Datum x_point(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *o; + int type1,j; + int32 *offsets1; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POINTTYPE) //point + { + PG_RETURN_FLOAT4( ((POINT3D *)o)->x) ; + } + + } + PG_RETURN_NULL(); +} + +//Y(GEOMETRY) -- find the first POINT(..) in GEOMETRY, returns its Y value. +//Return NULL if there is no POINT(..) in GEOMETRY +PG_FUNCTION_INFO_V1(y_point); +Datum y_point(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *o; + int type1,j; + int32 *offsets1; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POINTTYPE) //point + { + PG_RETURN_FLOAT4( ((POINT3D *)o)->y) ; + } + + } + PG_RETURN_NULL(); +} + +//Z(GEOMETRY) -- find the first POINT(..) in GEOMETRY, returns its Z value. +//Return NULL if there is no POINT(..) in GEOMETRY +PG_FUNCTION_INFO_V1(z_point); +Datum z_point(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *o; + int type1,j; + int32 *offsets1; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POINTTYPE) //point + { + PG_RETURN_FLOAT4( ((POINT3D *)o)->z) ; + } + + } + PG_RETURN_NULL(); +} + +//numpoints(GEOMETRY) -- find the first linestring in GEOMETRY, return +//the number of points in it. Return NULL if there is no LINESTRING(..) +//in GEOMETRY +PG_FUNCTION_INFO_V1(numpoints_linestring); +Datum numpoints_linestring(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *o; + int type1,j; + int32 *offsets1; + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == LINETYPE) //linestring + { + PG_RETURN_INT32( ((LINE3D *)o)->npoints) ; + } + + } + PG_RETURN_NULL(); +} + +//pointn(GEOMETRY,INTEGER) -- find the first linestring in GEOMETRY, +//return the point at index INTEGER (0 is 1st point). Return NULL if +//there is no LINESTRING(..) in GEOMETRY or INTEGER is out of bounds. +// keeps is3d flag + +PG_FUNCTION_INFO_V1(pointn_linestring); +Datum pointn_linestring(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + int32 wanted_index =PG_GETARG_INT32(1); + char *o; + int type1,j; + int32 *offsets1; + LINE3D *line; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == LINETYPE) //linestring + { + line = (LINE3D *)o; + if ( (wanted_index<0) || (wanted_index> (line->npoints-1) ) ) + PG_RETURN_NULL(); //index out of range + //get the point, make a new geometry + PG_RETURN_POINTER( + make_oneobj_geometry(sizeof(POINT3D), + (char *)&line->points[wanted_index], + POINTTYPE, geom->is3d) + ); + } + + } + PG_RETURN_NULL(); +} + +// exteriorRing(GEOMETRY) -- find the first polygon in GEOMETRY, return +//its exterior ring (as a linestring). Return NULL if there is no +//POLYGON(..) in GEOMETRY. + +PG_FUNCTION_INFO_V1(exteriorring_polygon); +Datum exteriorring_polygon(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + char *o; + int type1,j; + int32 *offsets1; + LINE3D *line; + POLYGON3D *poly; + POINT3D *pts; + int size_line; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POLYGONTYPE) //polygon object + { + poly = (POLYGON3D *)o; + pts = (POINT3D *) ( (char *)&(poly->npoints[poly->nrings] ) ); + pts = (POINT3D *) MAXALIGN(pts); + + line = make_line(poly->npoints[0], pts, &size_line); + + //get the ring, make a new geometry + PG_RETURN_POINTER( + make_oneobj_geometry(size_line, + (char *) line, + LINETYPE, geom->is3d) + ); + } + + } + PG_RETURN_NULL(); +} + +//NumInteriorRings(GEOMETRY) -- find the first polygon in GEOMETRY, +//return the number of interior rings. Return NULL if there is no +//POLYGON(..) in GEOMETRY. + +PG_FUNCTION_INFO_V1(numinteriorrings_polygon); +Datum numinteriorrings_polygon(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + char *o; + int type1,j; + int32 *offsets1; + POLYGON3D *poly; + + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POLYGONTYPE) //polygon object + { + poly = (POLYGON3D *)o; + PG_RETURN_INT32( poly->nrings -1 ) ; + } + } + PG_RETURN_NULL(); +} + +// InteriorRingN(GEOMETRY,INTEGER) -- find the first polygon in GEOMETRY, +//return the interior ring at index INTEGER (as a linestring). Return +//NULL if there is no POLYGON(..) in GEOMETRY or INTEGER is out of bounds. +// 1st ring = exerior ring + +PG_FUNCTION_INFO_V1(interiorringn_polygon); +Datum interiorringn_polygon(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + int32 wanted_index =PG_GETARG_INT32(1); + char *o; + int type1,j,t,point_offset; + int32 *offsets1; + POLYGON3D *poly; + LINE3D *line; + POINT3D *pts; + int size_line; + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + + if (type1 == POLYGONTYPE) //polygon object + { + poly = (POLYGON3D *)o; + + if ( (wanted_index<0) || (wanted_index> (poly->nrings-2) ) ) + PG_RETURN_NULL(); //index out of range + + + pts = (POINT3D *) ( (char *)&(poly->npoints[poly->nrings] ) ); + pts = (POINT3D *) MAXALIGN(pts); + + //find the 1st point in wanted ring + point_offset=0; + for (t=0; t< (wanted_index+1); t++) + { + point_offset += poly->npoints[t]; + } + + line = make_line(poly->npoints[wanted_index+1], &pts[point_offset], &size_line); + + //get the ring, make a new geometry + PG_RETURN_POINTER( + make_oneobj_geometry(size_line, + (char *) line, + LINETYPE, geom->is3d) + ); + } + } + PG_RETURN_NULL(); +} + + +// numgeometries(GEOMETRY) -- if GEOMETRY is a GEOMETRYCOLLECTION, return +//the number of geometries in it, otherwise return NULL. + +PG_FUNCTION_INFO_V1(numgeometries_collection); +Datum numgeometries_collection(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + + if (geom->type == COLLECTIONTYPE) + PG_RETURN_INT32( geom->nobjs ) ; + else + PG_RETURN_NULL(); +} + + +//geometryN(GEOMETRY, INTEGER) -- if GEOMETRY is a GEOMETRYCOLLECTION, +//return the sub-geometry at index INTEGER (0=first geometry), otherwise +//return NULL. NOTE: MULTIPOINT, MULTILINESTRING,MULTIPOLYGON are +//converted to sets of POINT,LINESTRING, and POLYGON so the index may +//change. + + +PG_FUNCTION_INFO_V1(geometryn_collection); +Datum geometryn_collection(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + int32 wanted_index =PG_GETARG_INT32(1); + int type; + int32 *offsets1; + char *o; + + + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + if (geom->type != COLLECTIONTYPE) + PG_RETURN_NULL(); + + if ( (wanted_index <0) || (wanted_index > (geom->nobjs-1) ) ) + PG_RETURN_NULL(); //bad index + + type = geom->objType[wanted_index]; + o = (char *) geom +offsets1[wanted_index] ; + + if (type == POINTTYPE) + { + PG_RETURN_POINTER( + make_oneobj_geometry(sizeof(POINT3D), + o, + POINTTYPE, geom->is3d) + ); + } + if (type == LINETYPE) + { + PG_RETURN_POINTER( + make_oneobj_geometry( size_subobject (o, LINETYPE), + o, + LINETYPE, geom->is3d) + ); + } + if (type == POLYGONTYPE) + { + PG_RETURN_POINTER( + make_oneobj_geometry( size_subobject (o, POLYGONTYPE), + o, + POLYGONTYPE, geom->is3d) + ); + } + + PG_RETURN_NULL(); +} + + +//force the geometry to be a geometrycollection type + +PG_FUNCTION_INFO_V1(force_collection); +Datum force_collection(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GEOMETRY *result; + + result = (GEOMETRY *) palloc(geom->size); + memcpy(result,geom, geom->size); + + result->type = COLLECTIONTYPE; + + PG_RETURN_POINTER(result); +} + + + +// point_inside_circle(geometry, Px, Py, d) +// returns true if there is a point in geometry whose distance to (Px,Py) is < d + +PG_FUNCTION_INFO_V1(point_inside_circle); +Datum point_inside_circle(PG_FUNCTION_ARGS) +{ + + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + char *o; + int type1,j; + int32 *offsets1; + POINT3D *pt; + double Px = PG_GETARG_FLOAT8(1); + double Py = PG_GETARG_FLOAT8(2); + double d = PG_GETARG_FLOAT8(3); + double dd = d*d; //d squared + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom + { + o = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + if (type1 == POINTTYPE) //point + { + //see if its within d of Px,Py + pt = (POINT3D *) o; + if ( ( (pt->x - Px)*(pt->x - Px) + (pt->y - Py)*(pt->y - Py) ) < dd) + { + PG_RETURN_BOOL(TRUE); + } + } + + } + PG_RETURN_BOOL(FALSE); +} \ No newline at end of file diff --git a/postgis_inout.c b/postgis_inout.c index ac814ff71..72b1f4e25 100644 --- a/postgis_inout.c +++ b/postgis_inout.c @@ -1358,6 +1358,7 @@ BOX3D *bbox_of_geometry(GEOMETRY *geom) if (geom->objType[i] == POINTTYPE) { +//printf("box of a point\n"); a_box = bbox_of_point( (POINT3D *) obj); result= union_box3d(a_box ,result); @@ -1366,6 +1367,7 @@ BOX3D *bbox_of_geometry(GEOMETRY *geom) } if (geom->objType[i] == LINETYPE) { +//printf("box of a line, # points = %i\n",((LINE3D *) obj)->npoints ); a_box = bbox_of_line( (LINE3D *) obj); result = union_box3d(a_box ,result); if (a_box != NULL) @@ -1373,6 +1375,7 @@ BOX3D *bbox_of_geometry(GEOMETRY *geom) } if (geom->objType[i] == POLYGONTYPE) { +//printf("box of a polygon\n"); a_box = bbox_of_polygon( (POLYGON3D *) obj); result =union_box3d(a_box ,result); if (a_box != NULL) @@ -2487,7 +2490,7 @@ char *to_wkb_collection(GEOMETRY *geom, bool flip_endian, int32 *end_size) pfree( sizes); //total size of the wkb - *end_size = total_size; + *end_size = total_size+9; @@ -2608,39 +2611,304 @@ char *to_wkb_sub(GEOMETRY *geom, bool flip_endian, int32 *wkb_size) } -//convert binary geometry into OGIS well know binary format with XDR (big endian) formating +//convert binary geometry into OGIS well know binary format with NDR (little endian) formating // see http://www.opengis.org/techno/specs/99-049.rtf page 3-24 for specification // // 3d geometries are encode as in OGR by adding 32768 to the type. Points are then 24bytes (X,Y,Z) // instead of 16 bytes (X,Y) - -PG_FUNCTION_INFO_V1(wkb_XDR); -Datum wkb_XDR(PG_FUNCTION_ARGS) +// +//dont do any flipping of endian asbinary_simple(GEOMETRY) +PG_FUNCTION_INFO_V1(asbinary_simple); +Datum asbinary_simple(PG_FUNCTION_ARGS) { GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - if (BYTE_ORDER == BIG_ENDIAN) - PG_RETURN_POINTER(to_wkb(geom, FALSE)); - else - PG_RETURN_POINTER(to_wkb(geom, TRUE)); + PG_RETURN_POINTER(to_wkb(geom, FALSE)); } - //convert binary geometry into OGIS well know binary format with NDR (little endian) formating // see http://www.opengis.org/techno/specs/99-049.rtf page 3-24 for specification // // 3d geometries are encode as in OGR by adding 32768 to the type. Points are then 24bytes (X,Y,Z) // instead of 16 bytes (X,Y) - -PG_FUNCTION_INFO_V1(wkb_NDR); -Datum wkb_NDR(PG_FUNCTION_ARGS) +// +//flip if required asbinary_specify(GEOMETRY,'xdr') or asbinary_specify(GEOMETRY,'ndr') +PG_FUNCTION_INFO_V1(asbinary_specify); +Datum asbinary_specify(PG_FUNCTION_ARGS) { - GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + text *type = PG_GETARG_TEXT_P(1); - if (BYTE_ORDER == LITTLE_ENDIAN) - PG_RETURN_POINTER(to_wkb(geom, FALSE)); + if ( ( strcmp(VARDATA(type) ,"xdr") == 0 ) || (strcmp(VARDATA(type) ,"XDR") == 0) ) + { +printf("requested XDR\n"); + if (BYTE_ORDER == BIG_ENDIAN) + PG_RETURN_POINTER(to_wkb(geom, FALSE)); + else + PG_RETURN_POINTER(to_wkb(geom, TRUE)); + } else - PG_RETURN_POINTER(to_wkb(geom, TRUE)); + { +printf("requested NDR\n"); + if (BYTE_ORDER == LITTLE_ENDIAN) + PG_RETURN_POINTER(to_wkb(geom, FALSE)); + else + PG_RETURN_POINTER(to_wkb(geom, TRUE)); + } +} + + + + + +//make a geometry with one obj in it (of size new_obj_size) +// type should be POINTTYPE, LINETYPE or POLYGONTYPE +// if you want to change the object's type to something else (ie GEOMETRYCOLLECTION), do +// that after with geom->type = GEOMETRYCOLLECTION +// this does calculate the bvol +// +// do not call this with type = GEOMETRYCOLLECTION + +GEOMETRY *make_oneobj_geometry(int sub_obj_size, char *sub_obj, int type, bool is3d) +{ + int size = sizeof(GEOMETRY) + 4 + sub_obj_size ; + GEOMETRY *result = palloc(size); + char *sub_obj_loc; + BOX3D *bbox; + + result->size = size; + result->nobjs = 1; + result->type = type; + result->is3d = is3d; + + result->objType[0] = type; + if (type == MULTIPOINTTYPE) + result->objType[0] = POINTTYPE; + if (type == MULTILINETYPE) + result->objType[0] = LINETYPE; + if (type == MULTIPOLYGONTYPE) + result->objType[0] = POLYGONTYPE; + + if (type == COLLECTIONTYPE) + { + pfree(result); + return(NULL); //error + } + + sub_obj_loc = (char *) &result->objType[2]; + sub_obj_loc = (char *) MAXALIGN(sub_obj_loc); + + result->objType[1] = sub_obj_loc - (char *) result; //result->objType[1] is where objOffset is + + //copy in the subobject + memcpy(sub_obj_loc , sub_obj, sub_obj_size); + + bbox = bbox_of_geometry(result); + memcpy(&result->bvol,bbox, sizeof(BOX3D) ); //make bounding box + + return result; +} + + +//find the size of the subobject and return it +int size_subobject (char *sub_obj, int type) +{ + if (type == POINTTYPE) + { + return (sizeof(POINT3D)); + } + if (type == LINETYPE) + { + return(sizeof(LINE3D) + sizeof(POINT3D) * ( ((LINE3D *)sub_obj)->npoints )); + } + if (type==POLYGONTYPE) + { + POLYGON3D *poly = (POLYGON3D *) sub_obj; + int t,points=0; + + for (t=0;tnrings;t++) + { + points += poly->npoints[t]; + } + if ( ( (long) ( &poly->npoints[poly->nrings] )) == (MAXALIGN(&poly->npoints[poly->nrings] ) ) ) + return (sizeof(POLYGON3D) + 4*(poly->nrings-1) + sizeof(POINT3D)*(points-1) ); //no extra align + else + return (sizeof(POLYGON3D) + 4*(poly->nrings-1) + sizeof(POINT3D)*(points-1) +4 ); + } + + return (-1);//unknown sub-object type +} + + +//produce a new geometry, which is the old geometry with another object stuck in it +// This will try to make the geometry's type is correct (move POINTTYPE to MULTIPOINTTYPE or +// change to GEOMETRYCOLLECTION) +// +// this does NOT calculate the bvol - you should set it with "bbox_of_geometry" +// +// type is the type of the subobject +// do not call this as with type = GEOMETRYCOLLECTION +// +// doesnt change the is3d flag + +GEOMETRY *add_to_geometry(GEOMETRY *geom,int sub_obj_size, char *sub_obj, int type) +{ + int size,t; + int size_obj,next_offset; + GEOMETRY *result; + int32 *old_offsets, *new_offsets; + BOX3D *bbox; + + //all the offsets could cause re-alignment problems, so need to deal with each on + size = geom->size +(4*geom->nobjs +1) /*byte align*/ + +sub_obj_size + 4 /*ObjType[]*/ +4 /*new offset*/; + + result = (GEOMETRY *) palloc(size); + result->size = size; + result->is3d = geom->is3d; + + //accidently sent in a single-entity type but gave it a multi-entity type + // re-type it as single-entity + if (type == MULTIPOINTTYPE) + type = POINTTYPE; + if (type == MULTILINETYPE) + type = LINETYPE; + if (type == MULTIPOLYGONTYPE) + type = POLYGONTYPE; + + + //simple conversion + if (geom->type == POINTTYPE) + { + if (type == POINTTYPE) + result->type = MULTIPOINTTYPE; + else + result->type = COLLECTIONTYPE; + } + if (geom->type == LINETYPE) + { + if (type == LINETYPE) + result->type = MULTILINETYPE; + else + result->type = COLLECTIONTYPE; + } + if (geom->type == POLYGONTYPE) + { + if (type == POLYGONTYPE) + result->type = MULTIPOLYGONTYPE; + else + result->type = COLLECTIONTYPE; + } + if (geom->type == COLLECTIONTYPE) + result->type = COLLECTIONTYPE; + + // now result geometry's type and sub-object's type is okay + // we have to setup the geometry + + result->nobjs = geom->nobjs+1; + + for (t=0; t< geom->nobjs; t++) + { + result->objType[t] = geom->objType[t]; + } + +//printf("about to copy geomes\n"); +//printf("result is at %p and is %i bytes long\n",result,result->size); +//printf("geom is at %p and is %i bytes long\n",geom,geom->size); + + old_offsets =& geom->objType[geom->nobjs] ; + new_offsets =& result->objType[result->nobjs] ; + next_offset = ( (char *) &new_offsets[result->nobjs] ) - ( (char *) result) ; + next_offset = MAXALIGN(next_offset); + + //have to re-set the offsets and copy in the sub-object data + for (t=0; t< geom->nobjs; t++) + { + //where is this going to go? + new_offsets[t] = next_offset; + + size_obj = size_subobject ((char *) (((char *) geom) + old_offsets[t] ), geom->objType[t]); + + next_offset += size_obj; + next_offset = MAXALIGN(next_offset); // make sure its aligned properly + + +//printf("coping %i bytes from %p to %p\n", size_obj,( (char *) geom) + old_offsets[t],((char *) result) + new_offsets[t] ); + memcpy( ((char *) result) + new_offsets[t] , ( (char *) geom) + old_offsets[t], size_obj); +//printf("copying done\n"); + + } + +//printf("copying in new object\n"); + + //now, put in the new data + result->objType[ result->nobjs -1 ] = type; + new_offsets[ result->nobjs -1 ] = next_offset; + memcpy( ((char *) result) + new_offsets[result->nobjs-1] ,sub_obj , sub_obj_size); + +//printf("calculating bbox\n"); + + bbox = bbox_of_geometry(result); + memcpy(&result->bvol,bbox, sizeof(BOX3D) ); //make bounding box + +//printf("returning\n"); + + return result; +} + + + +//make a polygon obj +// size is return in arg "size" +POLYGON3D *make_polygon(int nrings, int *pts_per_ring, POINT3D *pts, int npoints, int *size) +{ + POLYGON3D *result; + int t; + POINT3D *inside_poly_pts; + + + *size = sizeof(POLYGON3D) + 4 /*align*/ + + 4*(nrings-1)/*npoints struct*/ + + sizeof(POINT3D) *(npoints-1) /*points struct*/ ; + + result= (POLYGON3D *) palloc (*size); + result->nrings = nrings; + + + for (t=0;tnpoints[t] = pts_per_ring[t]; + } + + inside_poly_pts = (POINT3D *) ( (char *)&(result->npoints[result->nrings] ) ); + inside_poly_pts = (POINT3D *) MAXALIGN(inside_poly_pts); + + memcpy(inside_poly_pts, pts, npoints *sizeof(POINT3D) ); + + return result; +} + +void set_point( POINT3D *pt,double x, double y, double z) +{ + pt->x = x; + pt->y = y; + pt->z = z; } + +//make a line3d object +// return size of structure in 'size' +LINE3D *make_line(int npoints, POINT3D *pts, int *size) +{ + LINE3D *result; + + *size = sizeof(LINE3D) + (npoints-1)*sizeof(POINT3D); + + result= (LINE3D *) palloc (*size); + + result->npoints = npoints; + memcpy( result->points, pts, npoints*sizeof(POINT3D) ); + + return result; +} \ No newline at end of file diff --git a/postgis_proj.c b/postgis_proj.c new file mode 100644 index 000000000..6a80d8411 --- /dev/null +++ b/postgis_proj.c @@ -0,0 +1,338 @@ +#include "postgres.h" + + +#include +#include +#include +#include +#include + +#include "access/gist.h" +#include "access/itup.h" +#include "access/rtree.h" + + +#include "fmgr.h" + + +#include "postgis.h" +#include "utils/elog.h" + +#define SHOW_DIGS_DOUBLE 15 +#define MAX_DIGS_DOUBLE (SHOW_DIGS_DOUBLE + 6 + 1 + 3 +1) + + +// distance from -126 49 to -126 49.011096139863 in 'SPHEROID["GRS_1980",6378137,298.257222101]' is 1234.000 + + + + + +//use the WKT definition of an ellipsoid +// ie. SPHEROID["name",A,rf] or SPHEROID("name",A,rf) +// SPHEROID["GRS_1980",6378137,298.257222101] +// wkt says you can use "(" or "[" + +PG_FUNCTION_INFO_V1(ellipsoid_in); +Datum ellipsoid_in(PG_FUNCTION_ARGS) +{ + char *str = PG_GETARG_CSTRING(0); + SPHEROID *sphere = (SPHEROID *) palloc(sizeof(SPHEROID)); + int nitems; + double rf; + + + memset(sphere,0, sizeof(SPHEROID)); + + if (strstr(str,"SPHEROID") != str ) + { + elog(ERROR,"SPHEROID parser - doesnt start with SPHEROID"); + pfree(sphere); + PG_RETURN_NULL(); + } + + nitems = sscanf(str,"SPHEROID[\"%19[^\"]\",%lf,%lf]",sphere->name,&sphere->a,&rf); + + if ( nitems==0) + nitems = sscanf(str,"SPHEROID(\"%19[^\"]\",%lf,%lf)",sphere->name,&sphere->a,&rf); + + if (nitems != 3) + { + elog(ERROR,"SPHEROID parser - couldnt parse the spheroid"); + pfree(sphere); + PG_RETURN_NULL(); + } + + sphere->f = 1.0/rf; + sphere->b = sphere->a - (1.0/rf)*sphere->a; + sphere->e_sq = ((sphere->a*sphere->a) - (sphere->b*sphere->b) )/ (sphere->a*sphere->a); + sphere->e = sqrt(sphere->e_sq); + + PG_RETURN_POINTER(sphere); + +} + + + +PG_FUNCTION_INFO_V1(ellipsoid_out); +Datum ellipsoid_out(PG_FUNCTION_ARGS) +{ + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(0); + char *result; + + result = palloc(MAX_DIGS_DOUBLE + MAX_DIGS_DOUBLE + 20 + 9 + 2); + + sprintf(result,"SPHEROID(\"%s\",%.15g,%.15g)", sphere->name,sphere->a, 1.0/sphere->f); + + PG_RETURN_CSTRING(result); +} + +//support function for distance calc + //code is taken from David Skae + //Geographic Data BC, Province of British Columbia, Canada. + // Thanks to GDBC and David Skae for allowing this to be + // put in PostGIS. +double deltaLongitude(double azimuth, double sigma, double tsm,SPHEROID *sphere) +{ + // compute the expansion C + double das,C; + double ctsm,DL; + + das = cos(azimuth)*cos(azimuth); + C = sphere->f/16.0 * das * (4.0 + sphere->f * (4.0 - 3.0 * das)); + // compute the difference in longitude + + ctsm = cos(tsm); + DL = ctsm + C * cos(sigma) * (-1.0 + 2.0 * ctsm*ctsm); + DL = sigma + C * sin(sigma) * DL; + return (1.0 - C) * sphere->f * sin(azimuth) * DL; +} + + +//support function for distance calc + //code is taken from David Skae + //Geographic Data BC, Province of British Columbia, Canada. + // Thanks to GDBC and David Skae for allowing this to be + // put in PostGIS. +double mu2(double azimuth,SPHEROID *sphere) +{ + double e2; + + e2 = sqrt(sphere->a*sphere->a-sphere->b*sphere->b)/sphere->b; + return cos(azimuth)*cos(azimuth) * e2*e2; +} + + +//support function for distance calc + //code is taken from David Skae + //Geographic Data BC, Province of British Columbia, Canada. + // Thanks to GDBC and David Skae for allowing this to be + // put in PostGIS. +double bigA(double u2) +{ + return 1.0 + u2/256.0 * (64.0 + u2 * (-12.0 + 5.0 * u2)); +} + + +//support function for distance calc + //code is taken from David Skae + //Geographic Data BC, Province of British Columbia, Canada. + // Thanks to GDBC and David Skae for allowing this to be + // put in PostGIS. +double bigB(double u2) +{ + return u2/512.0 * (128.0 + u2 * (-64.0 + 37.0 * u2)); +} + + +//given 2 lat/longs and ellipse, find the distance +// note original r = 1st, s=2nd location +double distance_ellipse(double lat1, double long1, + double lat2, double long2, + SPHEROID *sphere) +{ + //code is taken from David Skae + //Geographic Data BC, Province of British Columbia, Canada. + // Thanks to GDBC and David Skae for allowing this to be + // put in PostGIS. + + double L1,L2,sinU1,sinU2,cosU1,cosU2; + double dl,dl1,dl2,dl3,cosdl1,sindl1; + double cosSigma,sigma,azimuthEQ,tsm; + double u2,A,B; + double dsigma; + + L1 = atan((1.0 - sphere->f ) * tan( lat1) ); + L2 = atan((1.0 - sphere->f ) * tan( lat2) ); + sinU1 = sin(L1); + sinU2 = sin(L2); + cosU1 = cos(L1); + cosU2 = cos(L2); + + dl = long2- long1; + dl1 = dl; + cosdl1 = cos(dl); + sindl1 = sin(dl); + do { + cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosdl1; + sigma = acos(cosSigma); + azimuthEQ = asin((cosU1 * cosU2 * sindl1)/sin(sigma)); + tsm = acos(cosSigma - (2.0 * sinU1 * sinU2)/(cos(azimuthEQ)*cos(azimuthEQ))); + dl2 = deltaLongitude(azimuthEQ, sigma, tsm,sphere); + dl3 = dl1 - (dl + dl2); + dl1 = dl + dl2; + cosdl1 = cos(dl1); + sindl1 = sin(dl1); + } while (fabs(dl3) > 1.0e-032); + + // compute expansions A and B + u2 = mu2(azimuthEQ,sphere); + A = bigA(u2); + B = bigB(u2); + + // compute length of geodesic + dsigma = B * sin(sigma) * (cos(tsm) + (B*cosSigma*(-1.0 + 2.0 * (cos(tsm)*cos(tsm))))/4.0); + return sphere->b * (A * (sigma - dsigma)); +} + + +double length2d_ellipse_linestring(LINE3D *line, SPHEROID *sphere) +{ + + int i; + POINT3D *frm,*to; + double dist = 0.0; + + if (line->npoints <2) + return 0.0; //must have >1 point to make sense + + frm = &line->points[0]; + + for (i=1; inpoints;i++) + { + to = &line->points[i]; + + dist += distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0, + to->y*M_PI/180.0 , to->x*M_PI/180.0, + sphere); + + frm = to; + } + return dist; +} + +double length3d_ellipse_linestring(LINE3D *line, SPHEROID *sphere) +{ + int i; + POINT3D *frm,*to; + double dist = 0.0; + double dist_ellipse; + + if (line->npoints <2) + return 0.0; //must have >1 point to make sense + + frm = &line->points[0]; + + for (i=1; inpoints;i++) + { + to = &line->points[i]; + + dist_ellipse = distance_ellipse(frm->y*M_PI/180.0 , frm->x*M_PI/180.0, + to->y*M_PI/180.0 , to->x*M_PI/180.0, + sphere); + + dist += sqrt(dist_ellipse*dist_ellipse + (frm->z*frm->z) ); + + frm = to; + } + return dist; + + +} + + +// length_ellipsoid(GEOMETRY, SPHEROID) +// find the "length of a geometry" +// length2d(point) = 0 +// length2d(line) = length of line +// length2d(polygon) = 0 +// uses ellipsoidal math to find the distance +//// x's are longitude, and y's are latitude - both in decimal degrees + +PG_FUNCTION_INFO_V1(length_ellipsoid); +Datum length_ellipsoid(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); + + int32 *offsets1; + char *o1; + int32 type1,j; + LINE3D *line; + double dist = 0.0; + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o1 = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + if (type1 == LINETYPE) //LINESTRING + { + line = (LINE3D *) o1; + dist += length2d_ellipse_linestring(line,sphere); + } + } + PG_RETURN_FLOAT4(dist); + + + +} + +// length3d_ellipsoid(GEOMETRY, SPHEROID) +// find the "length of a geometry" +// length3d(point) = 0 +// length3d(line) = length of line +// length3d(polygon) = 0 +// uses ellipsoidal math to find the distance on the XY plane, then +// uses simple Pythagoras theorm to find the 3d distance on each segment +// x's are longitude, and y's are latitude - both in decimal degrees + +PG_FUNCTION_INFO_V1(length3d_ellipsoid); +Datum length3d_ellipsoid(PG_FUNCTION_ARGS) +{ + GEOMETRY *geom = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + SPHEROID *sphere = (SPHEROID *) PG_GETARG_POINTER(1); + + int32 *offsets1; + char *o1; + int32 type1,j; + LINE3D *line; + double dist = 0.0; + + offsets1 = (int32 *) ( ((char *) &(geom->objType[0] ))+ sizeof(int32) * geom->nobjs ) ; + + + //now have to do a scan of each object + + for (j=0; j< geom->nobjs; j++) //for each object in geom1 + { + o1 = (char *) geom +offsets1[j] ; + type1= geom->objType[j]; + if (type1 == LINETYPE) //LINESTRING + { + line = (LINE3D *) o1; + dist += length3d_ellipse_linestring(line,sphere); + } + } + PG_RETURN_FLOAT4(dist); +} + + + + + + -- 2.49.0