Commit fc41e13f authored by Reza  ANSARI's avatar Reza ANSARI
Browse files

correction et amelioration ds l'interface a shapelib et ajout option proj= ds rdshp.cc, Reza

parent 7ac3a5ab
......@@ -151,7 +151,7 @@ public:
//! return the Y-ccordinate for vertex k
inline double getY(int k) { return psShape_->padfY[k]; }
//! return the Z-ccordinate for vertex k
inline double getZ(int k) { return psShape_->padfY[k]; }
inline double getZ(int k) { return psShape_->padfZ[k]; }
//! prints the file header information (number of shapes, type and bounds)
ostream& printHeader(ostream& os)
......
......@@ -135,8 +135,11 @@ DataTable SHPReader::FillDT(SphProjectionBase* proj)
dt.AddFloatColumn("z");
dt.AddFloatColumn("xproj");
dt.AddFloatColumn("yproj");
MuTyV rec[15];
dt.AddFloatColumn("delx");
dt.AddFloatColumn("dely");
dt.AddFloatColumn("delxproj");
dt.AddFloatColumn("delyproj");
MuTyV rec[20];
SLI_SHP_Interface& shp=(*shp_p_);
if (prtlev_>1) shp.printHeader(std::cout);
......@@ -149,6 +152,7 @@ DataTable SHPReader::FillDT(SphProjectionBase* proj)
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();
......@@ -189,22 +193,44 @@ DataTable SHPReader::FillDT(SphProjectionBase* proj)
<<" -> "<<shp.ShapePartType2Name(ptyp)<<std::endl;
continue;
}
rec[5]=iPart;
int npts=0;
double x,y,xlast,ylast,xp,yp,xplast,yplast;
double delx=0, dely=0, delxp=0, delyp=0;
bool fgokproj=false;
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.;
x=shp.getX(j);
y=shp.getY(j);
rec[6]=x;
rec[7]=y;
rec[8]=shp.getZ(j);;
for(size_t ii=9; ii<20; ii++) rec[ii]=-99999999.;
fgokproj=false;
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;
fgokproj=proj->SphToPlane(ll,xp,yp);
if (fgokproj) {
rec[9]=xp; rec[10]=yp;
}
}
if (npts>0) {
rec[11]=x-xlast; rec[12]=y-ylast;
if (fgokproj) {
rec[13]=xp-xplast; rec[14]=yp-yplast;
}
}
else {
rec[11]=rec[12]=rec[13]=rec[14]=0.;
}
dt.AddLine(rec);
xlast=x; ylast=y;
if (fgokproj) {
xplast=xp; yplast=yp;
}
npts++;
}
}
}
......
......@@ -157,7 +157,7 @@ public:
//! return the Y-ccordinate for vertex k
inline double getY(int k) { return psShape_->padfY[k]; }
//! return the Z-ccordinate for vertex k
inline double getZ(int k) { return psShape_->padfY[k]; }
inline double getZ(int k) { return psShape_->padfZ[k]; }
//! prints the file header information (number of shapes, type and bounds)
std::ostream& printHeader(std::ostream& os);
......
......@@ -18,15 +18,21 @@
\ingroup PrgUtil
\file rdshp.cc
\brief \b rdshp: Read content of shapefile (geographic objects) and create GeoShapeCollection object
Scan HDF5 files and prints information on each Datast in file.
Uses the HDF5IOServer module.
Can also create a DataTable, with projected coordinates, using the specified spherical projection.
Uses the ShapeLib module.
\verbatim
sh> rdshp -h
SOPHYA Version 2.5 Revision 60 (V_Avr2022) -- May 4 2022 10:00:00 gcc Apple LLVM 13.0.0 (clang-1300.0.29.30)
SOPHYA Version 2.5 Revision 60 (V_Avr2022) -- May 12 2022 23:52:06 gcc 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)
Usage: rdshp [flags] inPath outPPFName [attname1 attname2 ...]
flags = -dbf -prt=lev -dt -nm=nameatt -id=idatt
flags = -proj=projType -dbf -prt=lev -dt -nm=nameatt -id=idatt
-proj=projType : projType=mollweide,sinus,lambertcyl,mercator
-dbf : open also attribute file (.dbf) - default: NO
-prt=lev : define print level (default=0)
-dt : fill and save DataTable instead of GeoShapeCollection object
-nm=nameatt : attribute for ShapeObject name to be extracted from attributes
-id=idatt : attribute for ShapeObject id to be extracted from attributes
\endverbatim
*/
......@@ -41,7 +47,13 @@ int main(int narg, char *arg[])
{
if ((narg < 2) || (strcmp(arg[1],"-h") == 0) ) {
cout << " Usage: rdshp [flags] inPath outPPFName [attname1 attname2 ...] \n"
<< " flags = -dbf -prt=lev -dt -nm=nameatt -id=idatt \n"<<endl;
<< " flags = -proj=projType -dbf -prt=lev -dt -nm=nameatt -id=idatt \n"
<< " -proj=projType : projType=mollweide,sinus,lambertcyl,mercator \n"
<< " -dbf : open also attribute file (.dbf) - default: NO \n"
<< " -prt=lev : define print level (default=0) \n"
<< " -dt : fill and save DataTable instead of GeoShapeCollection object \n"
<< " -nm=nameatt : attribute for ShapeObject name to be extracted from attributes \n"
<< " -id=idatt : attribute for ShapeObject id to be extracted from attributes \n"<<endl;
return(0);
}
......@@ -53,6 +65,7 @@ int main(int narg, char *arg[])
string outppf;
int aoff=0;
string nameatt, idatt;
string prjtype="mollweide";
for (int k=1; k<narg; k++) {
string aopt = arg[k];
if (aopt=="-dbf") {
......@@ -70,6 +83,9 @@ int main(int narg, char *arg[])
else if (aopt.substr(0,4)=="-id=") {
idatt=aopt.c_str()+4; aoff++;
}
else if (aopt.substr(0,6)=="-proj=") {
prjtype=aopt.c_str()+6; aoff++;
}
else break;
}
if (narg<3+aoff) {
......@@ -94,10 +110,13 @@ int main(int narg, char *arg[])
GeoShapeCollection gsc;
POutPersist po(outppf);
if (fgdt) {
MollweideProjection mollproj(Angle(0.));
cout << " --- Creating DataTable from shapefile..."<<endl;
dt=shprd.FillDT(&mollproj);
SphProjectionBase* prj=SphProjectionBase::getProjector(prjtype);
// MollweideProjection mollproj(Angle(0.));
cout << " --- Creating DataTable from shapefile - Projection type="<<prjtype<<endl;
dt=shprd.FillDT(prj);
po<<dt;
dt.SetShowMinMaxFlag(true);
cout<<dt;
}
else {
......
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