diff --git a/Taudem5PCVS2010/AreaD8/AreaD8.vcxproj b/Taudem5PCVS2010/AreaD8/AreaD8.vcxproj index 379e4094..6e1b15f2 100644 --- a/Taudem5PCVS2010/AreaD8/AreaD8.vcxproj +++ b/Taudem5PCVS2010/AreaD8/AreaD8.vcxproj @@ -127,7 +127,7 @@ ProgramDatabase - libxml2.lib;gdal_i.lib;expat.lib;msmpi.lib;%(AdditionalDependencies) + gdal_i.lib;msmpi.lib;%(AdditionalDependencies) C:\GDAL\lib;C:\Program Files\Microsoft HPC Pack 2012\Lib\amd64;%(AdditionalLibraryDirectories) true Console diff --git a/pyfiles/TauDEM Tools 10 2.tbx b/pyfiles/TauDEM Tools 10 2.tbx new file mode 100644 index 00000000..94cea6b8 Binary files /dev/null and b/pyfiles/TauDEM Tools 10 2.tbx differ diff --git a/src/ReadOutlets.cpp b/src/ReadOutlets.cpp index 9c15d845..6637a01c 100644 --- a/src/ReadOutlets.cpp +++ b/src/ReadOutlets.cpp @@ -38,17 +38,14 @@ email: dtarb@usu.edu // This software is distributed from http://hydrology.usu.edu/taudem/ -// 1/25/14. Modified to use shapelib by Chris George +// 5/2/16. Modified to use OGR by Nazmus Sazib #include #include #include "commonLib.h" #include "ogr_api.h" - - - -int readoutlets(char *outletsds,char *lyrname, int uselayername,int outletslyr,OGRSpatialReferenceH hSRSRaster,int *noutlets, double*& x, double*& y) +int readoutlets(char *outletsds,char *lyrname, int uselayername,int outletslyr,OGRSpatialReferenceH hSRSRaster,int *noutlets, double*& x, double*& y,int*& id) { // initializing datasoruce,layer,feature, geomtery, spatial reference @@ -60,107 +57,15 @@ int readoutlets(char *outletsds,char *lyrname, int uselayername,int outletslyr,O OGRFeatureH hFeature1; OGRGeometryH geometry, line; OGRSpatialReferenceH hRSOutlet; - // regiser all ogr driver related to OGR + // register all ogr driver related to OGR OGRRegisterAll(); - // open datasource - hDS1 = OGROpen(outletsds, FALSE,NULL); - if( hDS1 == NULL ) - { - printf( "Error Opening datasource %s.\n",outletsds ); - exit( 1 ); - } - - //get layer from layer name - if(uselayername==1) { hLayer1 = OGR_DS_GetLayerByName(hDS1,lyrname);} - //get layerinfo from layer number - else { hLayer1 = OGR_DS_GetLayer(hDS1,outletslyr);} // get layerinfo from layername - - if(hLayer1 == NULL)getlayerfail(hDS1,outletsds,outletslyr); - OGRwkbGeometryType gtype; - gtype=OGR_L_GetGeomType(hLayer1); - if(gtype != wkbPoint)getlayerfail(hDS1,outletsds,outletslyr); - // get spatial reference - hRSOutlet = OGR_L_GetSpatialRef(hLayer1); - - //if there is spatial reference then write warnings - - if(hRSOutlet!=NULL) { - - const char* epsgAuthorityIdRaster; - const char* epsgAuthorityIdOutlet; - int pj_raster=OSRIsProjected(hSRSRaster); // find if projected or not - int pj_outlet=OSRIsProjected(hRSOutlet); - OSRAutoIdentifyEPSG(hSRSRaster); //identify EPSG code - OSRAutoIdentifyEPSG(hRSOutlet); - const char *sprs; - if(pj_raster==0) {sprs="GEOGCS";} else { sprs="PROJCS"; } - if (pj_raster==pj_outlet){ - epsgAuthorityIdRaster=OSRGetAuthorityCode(hSRSRaster,sprs);// get EPSG code. TODO. Make sure these functions do not fail if there is no EPSG code - epsgAuthorityIdOutlet=OSRGetAuthorityCode(hRSOutlet,sprs); - - if(atoi(epsgAuthorityIdRaster)!=atoi( epsgAuthorityIdOutlet)){ - printf( "Warning: EPSG code of Outlet shapefile and Raster data are different.\n" ); - // TODO - Print the WKT and EPSG code of each. If no spatial reference information, print unknown - // TODO - Test how this works if spatial reference information is incomplete, and create at least one of the unit test functions with a shapefile without a .prj file, and one of the unit test functions a raster without a projection (eg an ASCII file) - } - } - - else { - printf( "Warning: Spatial References of Outlet shapefile and Raster data are different.\n" ); - // TODO - Print the WKT of each. The general idea is that if these match, do not print anything. - // If these do not match give the user a warning. Only give an error if the program can not proceed, such as would be the case if rows and columns did not match. - } - - } - - - long countPts=0; - // count number of feature - countPts=OGR_L_GetFeatureCount(hLayer1,0); - // get schema i.e geometry, properties (e.g. ID) - hFDefn1 = OGR_L_GetLayerDefn(hLayer1); - x = new double[countPts]; - y = new double[countPts]; - int iField; - int nxy=0; - // loop through each feature and get the latitude and longitude for each feature - OGR_L_ResetReading(hLayer1); - while( (hFeature1 = OGR_L_GetNextFeature(hLayer1)) != NULL ){ - //for( int j=0; j= 0) id = new int[countPts]; // loop through each feature and get lat,lon and id information @@ -226,20 +154,30 @@ int readoutlets(char *outletsds,char *lyrname,int uselayername,int outletslyr,OG x[nxy] = OGR_G_GetX(geometry, 0); y[nxy] = OGR_G_GetY(geometry, 0); int idfld =OGR_F_GetFieldIndex(hFeature1,"id"); - if (idfld >= 0) { - hFieldDefn = OGR_FD_GetFieldDefn( hFDefn1,idfld); // get field definiton based on index - if( OGR_Fld_GetType(hFieldDefn) == OFTInteger ) { + hFieldDefn1 = OGR_FD_GetFieldDefn( hFDefn1,idfld); // get field definiton based on index + if( OGR_Fld_GetType(hFieldDefn1) == OFTInteger ) { id[nxy] =OGR_F_GetFieldAsInteger( hFeature1, idfld );} // get id value } - - nxy++; // count number of outlets point - OGR_F_Destroy( hFeature1 ); // destroy feature - } + else { + id[nxy]=1;// if there is no id field + } + nxy++; // count number of outlets point + OGR_F_Destroy( hFeature1 ); // destroy feature + } *noutlets=nxy; - OGR_DS_Destroy( hDS1); // destroy data source + OGR_DS_Destroy( hDS1); // destroy data source return 0; } + + +int readoutlets(char *outletsds,char *lyrname,int uselayername,int outletslyr,OGRSpatialReferenceH hSRSRaster, int *noutlets, double*& x, double*& y) + +{ + int *id; + int retval=readoutlets(outletsds,lyrname,uselayername,outletslyr,hSRSRaster, noutlets,x,y,id); + return retval; // Do this so that if ever we use the return from readoutlets with id it gets captured in the wrapper function +}