Commit 8a7bc7f3 authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

suite implementation class SHPReader ds ShapeLib, Reza

parent 277526f1
=========================================================================
S O P H Y A
SOftware for PHYsics Analysis
=========================================================================
ShapeLib module - Created in April 2022
o Main directory content ( components of libShapeLib.a )
---------------------------------------------------------
Files from shapelib c-library:
shapefil.h , dbfopen.c, safileio.c, sbnsearch.c, shpopen.c, shptree.c
and the c++ sophya interface :
sopshapelintf.h sopshapelintf.cc
The gbsopshapelintf.honly file is the original, inline only version of sopshapelintf.h
Kept for Guy Barrand
- shapelib web ressources:
http://shapelib.maptools.org
https://github.com/OSGeo/shapelib
- Possible data sources:
https://www.naturalearthdata.com
https://www.eea.europa.eu/data-and-maps/
o Tests/ subdirectory
----------------------
shapelib c-library executable programs and the gbread.cc test program,
and associated makefiles. o_makefile is the original, stand-alone makefile
......@@ -21,62 +21,80 @@ namespace SOPHYA {
\brief a class with static methods to read shapefiles and fill GeoShapeCollection object
\warning the path for .shp / .dbf file should be specified without the extension (.shp , .dbf)
\sa SOPHYA::GeoShapeCollection
\sa SOPHYA::SLI_SHP_Interface
\sa SOPHYA::SLI_DBF_Interface
*/
/* --Methode-- */
GeoShapeType SHPReader::convert2GST(int type)
SHPReader::SHPReader(const std::string& a_path, bool fgdbf)
: dbf_p_(nullptr), shp_p_(nullptr)
{
GeoShapeType gst=GST_unknown;
switch (type) {
case SHPT_POINT:
case SHPT_POINTZ:
case SHPT_POINTM:
gst=GST_point;
break;
case SHPT_ARC:
case SHPT_ARCZ:
case SHPT_ARCM:
gst=GST_Arc;
break;
case SHPT_POLYGON:
case SHPT_POLYGONZ:
case SHPT_POLYGONM:
gst=GST_Polygon;
break;
case SHPT_MULTIPOINT:
case SHPT_MULTIPOINTZ:
case SHPT_MULTIPOINTM:
gst=GST_MultiPoint;
break;
default:
gst=GST_unknown;
break;
}
return gst;
Init(a_path, fgdbf);
}
/* --Methode-- */
size_t SHPReader::Read(GeoShapeCollection & gsc, const std::string& shp_path)
SHPReader::SHPReader(const char* a_path, bool fgdbf)
{
if (prtlev_>0) std::cout<<"SHPReader::Read()/Info opening file:"<<shp_path<<std::endl;
SLI_SHP_Interface shp(shp_path);
if (prtlev_>2) shp.printHeader(std::cout);
std::string spath(a_path);
Init(spath, fgdbf);
}
/* --Methode-- */
SHPReader::~SHPReader()
{
if (dbf_p_) delete dbf_p_;
if (shp_p_) delete shp_p_;
}
/* --Methode-- */
void SHPReader::Init(const std::string& a_path, bool fgdbf)
{
shp_p_ = new SLI_SHP_Interface(a_path);
if (fgdbf) dbf_p_ = new SLI_DBF_Interface(a_path);
SetPrintLevel();
}
/* --Methode-- */
size_t SHPReader::Fill(GeoShapeCollection & gsc)
{
SLI_SHP_Interface& shp=(*shp_p_);
if (prtlev_>1) shp.printHeader(std::cout);
int jfname=-1;
int jfid=-1;
std::vector<int> jflist(list_att_.size());
if (dbf_p_) { // we have also to read an attribute file
if (prtlev_>2) dbf_p_->listFields(std::cout);
if (name_att_.length()>0)
jfname = dbf_p_->findFieldNum(name_att_);
if (id_att_.length()>0)
jfid = dbf_p_->findFieldNum(id_att_);
for(size_t k=0; k<list_att_.size(); k++)
jflist[k]=dbf_p_->findFieldNum(list_att_[k]);
}
for(int i=0;i<shp.getNShapes();i++) {
shp.selectShape(i);
int sht=shp.getShapeType();
GeoShapeType gst=convert2GST(sht);
if (sht == GST_unknown) {
std::cout << "SHPReader::Read() Shape type="<<sht<<" ("<<shp.ShapeType2Name(sht)<<") not handled "<<std::endl;
std::cout << "SHPReader::Fill()/Skipping Shape type="<<sht<<" ("<<shp.ShapeType2Name(sht)<<") not handled "<<std::endl;
continue;
}
std::string myname;
int myid=0;
if (dbf_p_) {
if (jfname>=0) myname = (std::string)(dbf_p_->getFieldValue(i,jfname));
if (jfid>=0) myid = (int)dbf_p_->getFieldValue(i,jfid);
}
if (prtlev_>0)
std::cout << "SHPReader::Read() shape["<<i<<"] Type="<<shp.ShapeType2Name(sht)<<std::endl;
std::cout << "SHPReader::Fill() shape["<<i<<"] Type="<<shp.ShapeType2Name(sht)<<" nVertices="
<<shp.getNVertices()<<" nParts="<<shp.getNParts()<<" Name="<<myname<<" Id="<<myid<<std::endl;
for(int iPart=0;iPart<shp.getNParts();iPart++) {
if (shp.getPartType(iPart) != SHPP_RING) {
int ptyp=shp.getPartType(iPart);
std::cout << " Part["<<iPart<<"] partType="<<ptyp<<" -> "<<shp.ShapePartType2Name(ptyp)<<std::endl;
std::cout << "SHPReader::Fill()/Skipping Part["<<iPart<<"] partType="<<ptyp
<<" -> "<<shp.ShapePartType2Name(ptyp)<<std::endl;
continue;
}
int jbeg, jend;
......@@ -89,11 +107,143 @@ size_t SHPReader::Read(GeoShapeCollection & gsc, const std::string& shp_path)
GeoCoord gc(ll, shp.getZ(j));
mgs.append(gc);
}
if (!dbf_p_) continue;
if ((jfname>=0)||(jfid>=0)) mgs.setNameId(myname, myid);
for(size_t k=0; k<list_att_.size(); k++) {
if ((jflist[k]<0)||(list_att_[k].length()<1)) continue;
MuTyV mtv=dbf_p_->getFieldValue(i, jflist[k]);
mgs.Info()[list_att_[k]]=mtv;
}
}
}
return gsc.size();
}
/* --Methode-- */
DataTable SHPReader::FillDT(SphProjectionBase* proj)
{
DataTable dt(64);
dt.AddIntegerColumn("num");
dt.AddStringColumn("name");
dt.AddStringColumn("continent");
dt.AddIntegerColumn("contId");
dt.AddIntegerColumn("stype");
dt.AddIntegerColumn("part");
dt.AddFloatColumn("x");
dt.AddFloatColumn("y");
dt.AddFloatColumn("z");
dt.AddFloatColumn("xproj");
dt.AddFloatColumn("yproj");
MuTyV rec[15];
SLI_SHP_Interface& shp=(*shp_p_);
if (prtlev_>1) shp.printHeader(std::cout);
int jfname=-1;
int jfcontinent=-1;
std::vector<int> jflist(list_att_.size());
if (dbf_p_) { // we have also to read an attribute file
if (prtlev_>2) dbf_p_->listFields(std::cout);
jfname = dbf_p_->findFieldNum("NAME");
jfcontinent = dbf_p_->findFieldNum("CONTINENT");
}
for(int i=0;i<shp.getNShapes();i++) {
shp.selectShape(i);
int sht=shp.getShapeType();
GeoShapeType gst=convert2GST(sht);
if (sht == GST_unknown) {
std::cout << "SHPReader::FillDT() Shape type="<<sht<<" ("<<shp.ShapeType2Name(sht)<<") not handled "<<std::endl;
continue;
}
std::string myname = (std::string)dbf_p_->getFieldValue(i,jfname);
std::string mycontinent = (std::string)dbf_p_->getFieldValue(i,jfcontinent);
std::cout << " Lecture Shape No "<<i<<" Type="<<shp.ShapeType2Name(shp.getShapeType())<<" nVertices="
<<shp.getNVertices()<<" nParts="<<shp.getNParts()<<" Name="<<myname<<" Continent="<<mycontinent<<endl;
int contId=0;
if (mycontinent == "Africa") contId=1;
else if (mycontinent == "Asia") contId=2;
else if (mycontinent == "Europe") contId=3;
else if (mycontinent == "North America") contId=4;
else if (mycontinent == "South America") contId=5;
else if (mycontinent == "Oceania") contId=6;
if (prtlev_>0)
std::cout << "SHPReader::FillDT() shape["<<i<<"] Type="<<shp.ShapeType2Name(sht)
<<" nVertices="<<shp.getNVertices()<<" nParts="<<shp.getNParts()
<<" Name="<<myname<<" Continent="<<mycontinent<<std::endl;
rec[0]=(int_4)i;
rec[1]=myname;
rec[2]=mycontinent;
rec[3]=contId;
rec[4]=sht;
for(int iPart=0;iPart<shp.getNParts();iPart++) {
if (shp.getPartType(iPart) != SHPP_RING) {
int ptyp=shp.getPartType(iPart);
std::cout << "SHPReader::FillDT()/Skipping Part["<<iPart<<"] partType="<<ptyp
<<" -> "<<shp.ShapePartType2Name(ptyp)<<std::endl;
continue;
}
int jbeg, jend;
shp.getVertexRange(iPart,jbeg,jend);
for(int j=jbeg;j<jend;j++) {
rec[6]=(float)shp.getX(j);
rec[7]=(float)shp.getY(j);
rec[8]=(float)shp.getZ(j);
rec[9]=rec[10]=-999999.;
if (proj) {
LongitudeLatitude ll(Angle(shp.getX(j), Angle::Degree), Angle(shp.getY(j), Angle::Degree), true);
double xp, yp;
bool fgokproj=proj->SphToPlane(ll,xp,yp);
if (fgokproj) {
rec[9]=xp; rec[10]=yp;
}
}
dt.AddLine(rec);
}
}
}
return dt;
}
/* --Methode-- */
GeoShapeType SHPReader::convert2GST(int type)
{
GeoShapeType gst=GST_unknown;
switch (type) {
case SHPT_POINT:
case SHPT_POINTZ:
case SHPT_POINTM:
gst=GST_point;
break;
case SHPT_ARC:
case SHPT_ARCZ:
case SHPT_ARCM:
gst=GST_Arc;
break;
case SHPT_POLYGON:
case SHPT_POLYGONZ:
case SHPT_POLYGONM:
gst=GST_Polygon;
break;
case SHPT_MULTIPOINT:
case SHPT_MULTIPOINTZ:
case SHPT_MULTIPOINTM:
gst=GST_MultiPoint;
break;
default:
gst=GST_unknown;
break;
}
return gst;
}
/*!
\class SLI_SHP_Interface
\ingroup ShapeLib
......
......@@ -11,26 +11,54 @@
#include <string>
#include <vector>
#include "shapefil.h"
#include "datatable.h"
#include "geoshape.h"
#include "projsph.h"
namespace SOPHYA {
//--------------------------------------------------------------
//--- reads in shapefile (and attribute file) and fills a GeoShapeCollection objet
//--- reads in shapefile (and attribute file) and fills a GeoShapeCollection objet
// Useful forward declarations
class SLI_SHP_Interface;
class SLI_DBF_Interface;
class SHPReader {
public:
//! constructor , path to the .shp file, and possibly the associated .dbf file if fgdbf==true
SHPReader(const std::string& a_path, bool fgdbf);
//! constructor , path to the .shp file, and possibly the associated .dbf file if fgdbf==true
SHPReader(const char* a_path, bool fgdbf);
//! destructor , closes the opened .shp and .dbf files
virtual ~SHPReader();
//! define the print level
inline void SetPrintLevel(int prtlev=0) { prtlev_=prtlev; }
//! define the attributes corresponding to the name and id of GeoShape objects
inline void setNameIdAtt(std::string& name, std::string& id)
{ name_att_=name; id_att_=id; }
//! define the list of attributes to be extracted from .dbf attribute file
inline void setAttributeList(std::vector<std::string> & list)
{ list_att_=list; }
/*! \brief read the shapefile content and fill the GeoShapeCollection object gsc
Note that new GeoShape objects are appended to the provided GeoShapeCollection object
*/
size_t Fill(GeoShapeCollection & gsc);
//! read the shapefile content and fill a DataTable object gsc. If provided, proj is used to compute the projected coordinates
DataTable FillDT(SphProjectionBase* proj=nullptr);
//! converts a shapelib type to SOPHYA GeoShapeType
static GeoShapeType convert2GST(int typ);
//! default constructor -
explicit SHPReader(int prtlev=0) : prtlev_(prtlev) { }
//! read the shape file specified shapefile (shp_path) and fills the GeoShapeCollection object gsc
size_t Read(GeoShapeCollection & gsc, const std::string& shp_path);
inline size_t Read(GeoShapeCollection & gsc, const char * path)
{
std::string spath=path;
return Read(gsc, spath);
}
protected:
void Init(const std::string& a_path, bool fgdbf);
SLI_SHP_Interface* shp_p_;
SLI_DBF_Interface* dbf_p_;
int prtlev_;
std::string name_att_; // attribute used for GeoShape name
std::string id_att_; // attribute used for GeoShape Id
std::vector<std::string> list_att_; // list of attributes to be extracted from .dbf file
};
//--------------------------------------------------------------
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment