-
Notifications
You must be signed in to change notification settings - Fork 116
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #64 from dtarb/63-EPSG-Code
63 epsg code
- Loading branch information
Showing
3 changed files
with
79 additions
and
141 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
|
|
@@ -38,17 +38,14 @@ email: [email protected] | |
|
|
||
| // 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 <stdio.h> | ||
| #include <string.h> | ||
| #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<countPts; j++) { //does not work for geojson file | ||
| geometry = OGR_F_GetGeometryRef(hFeature1); // get geometry type | ||
| x[nxy] = OGR_G_GetX(geometry, 0); | ||
| y[nxy] = OGR_G_GetY(geometry, 0); | ||
| OGR_F_Destroy( hFeature1); // destroy feature | ||
| nxy++; | ||
| } | ||
| *noutlets=nxy; // total number of outlets point | ||
|
|
||
| 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) | ||
|
|
||
| { | ||
|
|
||
| // initializing datasoruce,layer,feature, geomtery, spatial reference | ||
| OGRSFDriverH driver; | ||
| OGRDataSourceH hDS1; | ||
| OGRLayerH hLayer1; | ||
| OGRFeatureDefnH hFDefn1; | ||
| OGRFieldDefnH hFieldDefn1; | ||
| OGRFeatureH hFeature1; | ||
| OGRGeometryH geometry, line; | ||
| OGRSpatialReferenceH hRSOutlet; | ||
| OGRFieldDefnH hFieldDefn; | ||
| OGRRegisterAll(); | ||
| // open data soruce | ||
|
|
||
| // open data source | ||
| hDS1 = OGROpen(outletsds, FALSE, NULL ); | ||
| if( hDS1 == NULL ) | ||
| { | ||
| printf( "Error Opening in Shapefile .\n" ); | ||
| //exit( 1 ); | ||
| printf( "Error Opening OGR Data Source .\n" ); | ||
| return 1; | ||
| } | ||
|
|
||
| //get layer from layer name | ||
|
|
@@ -171,50 +76,73 @@ int readoutlets(char *outletsds,char *lyrname,int uselayername,int outletslyr,OG | |
| if(hLayer1 == NULL)getlayerfail(hDS1,outletsds,outletslyr); | ||
| OGRwkbGeometryType gtype; | ||
| gtype=OGR_L_GetGeomType(hLayer1); | ||
|
|
||
| // Test that the type is a point | ||
| if(gtype != wkbPoint)getlayerfail(hDS1,outletsds,outletslyr); | ||
| //OGR_L_ResetReading(hLayer1); | ||
|
|
||
| const char* RasterProjectionName; | ||
| const char* sprs; | ||
| const char* sprso; | ||
| const char* OutletProjectionName; | ||
| int pj_raster,pj_outlet; | ||
|
|
||
| // Spatial reference of outlet | ||
| hRSOutlet = OGR_L_GetSpatialRef(hLayer1); | ||
| //if there is spatial reference then write warnings | ||
| if(hSRSRaster!=NULL){ | ||
| pj_raster=OSRIsProjected(hSRSRaster); // find if projected or not | ||
| if(pj_raster==0) {sprs="GEOGCS";} else { sprs="PROJCS"; } | ||
| RasterProjectionName = OSRGetAttrValue(hSRSRaster,sprs,0); // get projection name | ||
| } | ||
| 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) | ||
| } | ||
| pj_outlet=OSRIsProjected(hRSOutlet); | ||
| if(pj_outlet==0) {sprso="GEOGCS";} else { sprso="PROJCS"; } | ||
| OutletProjectionName = OSRGetAttrValue(hRSOutlet,sprso,0); | ||
| } | ||
|
|
||
| //Write warnings where projections may not match | ||
| if(hRSOutlet!=NULL && hSRSRaster!=NULL){ | ||
|
|
||
| if (pj_raster==pj_outlet){ | ||
|
|
||
| int rc= strcmp(RasterProjectionName,OutletProjectionName); // compare string | ||
| if(rc!=0){ | ||
| printf( "Warning: Projection of Outlet feature and Raster data may be different.\n" ); | ||
| printf("Projection of Raster datasource %s.\n",RasterProjectionName); | ||
| printf("Projection of Outlet feature %s.\n",OutletProjectionName); | ||
| } | ||
| } | ||
|
|
||
| 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. | ||
| else { | ||
| printf( "Warning: Spatial References of Outlet feature and Raster data are different.\n" ); | ||
| printf("Projection of Raster datasource %s.\n",RasterProjectionName); | ||
| printf("Projection of Outlet feature %s.\n",OutletProjectionName); | ||
| } | ||
| } | ||
|
|
||
| else if(hSRSRaster==NULL && hRSOutlet!=NULL) { | ||
| printf( "Warning: Spatial Reference of Raster is missing.\n" ); | ||
| printf("Projection of Outlet feature %s.\n",OutletProjectionName); | ||
|
|
||
| } | ||
| else if(hSRSRaster!=NULL && hRSOutlet==NULL) { | ||
| printf( "Warning: Spatial Reference of Outlet feature is missing.\n" ); | ||
| printf("Projection of Raster datasource %s.\n",RasterProjectionName); | ||
| } | ||
| else { | ||
| printf( "Warning: Spatial References of Outlet feature and Raster data are missing.\n" ); | ||
| } | ||
|
|
||
|
|
||
|
|
||
| long countPts=0; | ||
| countPts=OGR_L_GetFeatureCount(hLayer1,0); // get feature count | ||
| hFDefn1 = OGR_L_GetLayerDefn(hLayer1); // get schema i.e geometry, properties (e.g. ID) | ||
| // 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; | ||
|
|
||
| //hFeature1 = OGR_L_GetNextFeature(hLayer1); | ||
| hFDefn1 = OGR_L_GetLayerDefn(hLayer1); | ||
| //hFeature1=OGR_L_GetFeature(hLayer1,0); // read first feature to get all field info | ||
| //int idfld =OGR_F_GetFieldIndex(hFeature1,"Id"); // get index for the 'id' field | ||
| //if (idfld >= 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 | ||
| } | ||
|
|
||