search

orbisgis / geoclimate

Geospatial processing toolbox for environmental and climate studies
GNU Lesser General Public License v3.0
59 stars 16 forks source link

Geoclimate SNAPSHOT-0.0.2 Nigthly- "The process has been stopped since the table *** has a no SRID" error #739

Closed sfaraut closed 2 years ago

sfaraut commented 2 years ago

Hello,

When I try to use all recent Geoclimate SNAPSHOT-0.0.2 Nightly versions with BDTOPO data downloaded from IGN, I have this error: The process has been stopped since the table BATI_INDIFFERENCIE has a no SRID.

C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate>"C:\Program Files\Java\jdk-11.0.13\bin\java.exe" -jar geoclimate-0.0.2-SNAPSHOT-2022-05-16.jar -f my_first_try_BDTOPO_Aigrefeuille-id_zone-converted-0.0.2.json -w BDTOPO_V2.2
[main] INFO org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - Starting to process insee id_zone 31480
[main] ERROR org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - The process has been stopped since the table BATI_INDIFFERENCIE has a no SRID
[main] ERROR org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - Cannot import preprocess.
The BDTOPO_V2.2 workflow has been successfully executed

I tried to convert DATA to other projections using "ogr2og" , and it makes no difference. Tried all RGF93 (Lambert93, CC43,...) and even "WGS 84 / World Mercator - EPSG:3395]"

Is only the projection of "sample_12174" to be recognized bhy Geoclimate?

For example with WGS 84 / World Mercator - EPSG:3395 :

Config file (my_first_try_BDTOPO_Aigrefeuille-id_zone-converted-0.0.2.json):

{
    "description": "Processing BD Topo v2 data",
    "output": {
        "folder": "C:\\Users\\faraut\\Documents\\DEVELOPPEMENT\\Geoclimate\\TEMP-RESULTS"    
    },
    "parameters": {        "rsu_indicators": {
            "indicatorUse": [
                "LCZ",
                "TEB",
                "UTRF"
            ],
            "svfSimplified": true,
            "estimateHeight": true
        },
        "grid_indicators": {
            "x_size": 100,
        "y_size": 100,
        "rowCol": false,
        "output" : "geojson",
        "indicators" :[
                 "BUILDING_FRACTION", 
                 "BUILDING_HEIGHT", 
                 "WATER_FRACTION",
                 "VEGETATION_FRACTION", 
                 "ROAD_FRACTION", 
                 "IMPERVIOUS_FRACTION", 
                 "LCZ_FRACTION"
             ]
        }
     },
    "input": { 

    "folder": {"path": "C:\\DATA\\IGN-DATA\\BDTOPO_2-2_EPSG3395_D031_2018-09-25-converted",
                   "locations": [ "31480" ],
      }
    }
}

Execution:

C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate>"C:\Program Files\Java\jdk-11.0.13\bin\java.exe" -jar geoclimate-0.0.2-SNAPSHOT-2022-05-16.jar -f my_first_try_BDTOPO_Aigrefeuille-id_zone-converted-0.0.2.json -w BDTOPO_V2.2
[main] INFO org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - Starting to process insee id_zone 31480
[main] ERROR org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - The process has been stopped since the table BATI_INDIFFERENCIE has a no SRID
[main] ERROR org.orbisgis.geoclimate.bdtopo_v2.BDTopo_V2 - Cannot import preprocess.
The BDTOPO_V2.2 workflow has been successfully executed

C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate>

Informations on shapefile C:\DATA\IGN-DATA\BDTOPO_2-2_EPSG3395_D031_2018-09-25-converted\BATI_INDIFFERENCIE.shp :

gdalsrsinfo -o wkt C:\DATA\IGN-DATA\BDTOPO_2-2_EPSG3395_D031_2018-09-25-converted\BATI_INDIFFERENCIE.shp

PROJCRS["WGS 84 / World Mercator",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["World Mercator",
        METHOD["Mercator (variant A)",
            ID["EPSG",9804]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["World - between 80┬░S and 84┬░N"],
        BBOX[-80,-180,84,180]],
    ID["EPSG",3395]]

Same thing with RGF93 / Lambert-93 - EPSG 2154, with in config file concerning new input folder:

    "input": { 

    "folder": {"path": "C:\\DATA\\IGN-DATA\\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted",
                   "locations": [ "31480" ],      }
     }

gdalsrsinfo -o wkt C:\DATA\IGN-DATA\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted\BATI_INDIFFERENCIE.shp

PROJCRS["RGF93 / Lambert-93",
    BASEGEOGCRS["RGF93",
        DATUM["Reseau Geodesique Francais 1993",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4171]],
    CONVERSION["Lambert-93",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",46.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",3,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",49,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",44,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",700000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",6600000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["France"],
        BBOX[41.15,-9.86,51.56,10.38]],
    ID["EPSG",2154]]

For reference, here is what I get with the "sample_12174" example which work (but with the road width bug corrected):

gdalsrsinfo -o wkt C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate\BD_TOPO_v2\sample_12174\ROUTE.shp

BOUNDCRS[
    SOURCECRS[
        PROJCRS["RGF93 / Lambert-93",
            BASEGEOGCRS["RGF93",
                DATUM["Reseau Geodesique Francais 1993",
                    ELLIPSOID["GRS 1980",6378137,298.257222101,
                        LENGTHUNIT["metre",1]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]],
                ID["EPSG",4171]],
            CONVERSION["Lambert-93",
                METHOD["Lambert Conic Conformal (2SP)",
                    ID["EPSG",9802]],
                PARAMETER["Latitude of false origin",46.5,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8821]],
                PARAMETER["Longitude of false origin",3,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8822]],
                PARAMETER["Latitude of 1st standard parallel",49,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8823]],
                PARAMETER["Latitude of 2nd standard parallel",44,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8824]],
                PARAMETER["Easting at false origin",700000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8826]],
                PARAMETER["Northing at false origin",6600000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8827]]],
            CS[Cartesian,2],
                AXIS["easting (X)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1]],
                AXIS["northing (Y)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1]],
            USAGE[
                SCOPE["unknown"],
                AREA["France"],
                BBOX[41.15,-9.86,51.56,10.38]],
            ID["EPSG",2154]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["RGF93 to WGS 84 (1)",
        VERSION["EPSG-Fra"],
        METHOD["Geocentric translations (geog2D domain)",
            ID["EPSG",9603]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        USAGE[
            SCOPE["Approximation at the +/- 1m level."],
            AREA["France"],
            BBOX[41.15,-9.86,51.56,10.38]],
        ID["EPSG",1671],
        REMARK["Parameter values from RGF93 to ETRS89 (1) (code 1591) assuming that ETRS89 is equivalent to WGS 84 within the accuracy of the transformation."]]]

These SRS information are more "complex" to understand, even its seems that main parameters are the same... What are the differences as they should be all of them some "RGF93 / Lambert-93" projection? And how this shapefile (or this projection) have been obtained?

I convert (by script) all data files using GDAL/OGR tools, on Windows. For forcing RGF93/Lamber-93 EPSG:2154" projection I used:

REM Convert to RGF93-Lambert93 - EPSG:2154 (NOT RECOGNIZED BY Geoclimate)
CD C:\DATA\IGN-DATA\BDTOPO_2-2_RGF93_D031_2018-09-25-converted
FOR %F IN (*.shp) DO ogr2ogr  -t_srs EPSG:2154 -a_srs EPSG:2154 "C:\DATA\IGN-DATA\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted\%F" "%F"

I don't understand why my "standard" projection (even WGS84) using "GDAL/OGR" don't work... Can't we use them do make batch process of data? Is it possible to know what projection can be recognized in Geoclimate by "CTS" library?

Or is it a bug to fix in SRID verification?

Thanks by advance,

Serge.

sfaraut commented 2 years ago

Same problems with SNAPSHOT-v0.0.2 Nightly downloaded this morning (18-05-2022)...

sfaraut commented 2 years ago

I gave before WTK's GDAL SRS informations (using gdalsrsinfo). Concerning data using RGF93-Lambert 93 projection, my BDTOPO data and "sample_12174" data are both seen as EPSG:2154 with "epsg" output type:

>gdalsrsinfo -o epsg  C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate\BD_TOPO_v2\sample_12174-orig\BATI_INDIFFERENCIE.shp

EPSG:2154

>gdalsrsinfo -o epsg  C:\\DATA\\IGN-DATA\\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted\BATI_INDIFFERENCIE.shp

EPSG:2154

Better, I could verify that their proj4 are also exactly the same:

>gdalsrsinfo -o proj4  C:\Users\faraut\Documents\DEVELOPPEMENT\Geoclimate\BD_TOPO_v2\sample_12174-orig\BATI_INDIFFERENCIE.shp

+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

>gdalsrsinfo -o proj4  C:\\DATA\\IGN-DATA\\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted\BATI_INDIFFERENCIE.shp

+proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

So it seems that Geoclimate (CTS?), somehow, is not able to parse correctly the "standard" .proj file of shapefiles provided by IGN or created by GDAL/OGR....

ebocher commented 2 years ago

GeoClimate uses H2GIS to process data and then CTS to read prj file. CTS parses only the WKT version 1 that describes the CRS information in the prj file. Could you please copy-paste, the input of one prj file (not the GDAL info).

ebocher commented 2 years ago

But this EPSG code is necessary to declare the SRID of the geometries in the database as PostGIS does

sfaraut commented 2 years ago

Sorry for the delay. My version of GDAL(gdalsrsinfo ) on Windows doesn't support WKT output format... Here is the raw .proj file content of BTDOPO V2 data converted to "EPSG21:54" (with GDAL/ogr2ogr):

C:\>type C:\\DATA\\IGN-DATA\\BDTOPO_2-2_EPSG2154_D031_2018-09-25-converted\BATI_INDIFFERENCIE.prj
PROJCS["RGF_1993_Lambert_93",GEOGCS["GCS_RGF_1993",DATUM["D_RGF_1993",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",700000.0],PARAMETER["False_Northing",6600000.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Standard_Parallel_1",49.0],PARAMETER["Standard_Parallel_2",44.0],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]
C:\>

This SRID is well identified by GDAL/gdalsrsinfo as as "EPSG:2154", just like GIS softwares (QGIS,...).

One way to solve this problem could be having a standard way/tool (in CLI mode like GDAL ?) to produice .proj files comptabible with Geoclimate/CTS. If not implementing a way to work as GDAL/gdalsrsinfo in CTS...

Tests must still have to be done on other system than Windows... or latest versions ( A couple of improvements... (PR #740))

sfaraut commented 2 years ago

To be complete, the original BDTOPO V2 files .proj files contain:

C:\>type C:\DATA\IGN-DATA\BDTOPO_2-2_TOUSTHEMES_SHP_LAMB93_D031_2018-09-25\BATI_INDIFFERENCIE.prj
PROJCS["RGF93_Lambert_93",GEOGCS["GCS_RGF_1993",DATUM["D_RGF_1993",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",700000.0],PARAMETER["False_Northing",6600000.0],PARAMETER["Central_Meridian",3.0],PARAMETER["Standard_Parallel_1",44.0],PARAMETER["Standard_Parallel_2",49.0],PARAMETER["Latitude_Of_Origin",46.5],UNIT["Meter",1.0]]
C:\>

The only difference with converted to "EPSG:2154" version is only "RGF93_Lambert_93" in place of "RGF_1993_Lambert_93". All other parameters are the same.

ebocher commented 2 years ago

Here is the raw .proj file content of BTDOPO V2 data converted to "EPSG21:54" (with GDAL/ogr2ogr):

As you can see the prj file doesn't contain any information about the epsg. GDal uses PROJ4 lib that offers a mecanism to retrieve the best epsg code according the CRS parameters. QGIS is based on GDAL so yes it detects the EPSG code.

One way to solve this problem could be having a standard way/tool (in CLI mode like GDAL ?) to produice .proj files comptabible with Geoclimate/CTS. If not implementing a way to work as GDAL/gdalsrsinfo in CTS...

It's a big job. CTS is a free software so a PR is welcome to solve this problem.

Tests must still have to be done on other system than Windows... or latest versions ( A couple of improvements... (PR https://github.com/orbisgis/geoclimate/pull/740))

Same comments as before. Any PR that can improved GeoClimate and help the community are welcome.

The only difference with converted to "EPSG:2154" version is only "RGF93_Lambert_93" in place of "RGF_1993_Lambert_93". All other parameters are the same.

A the end of prj formatting yes.

ebocher commented 2 years ago

The SRID error with non standardized prj has been fixed on H2GIS side see https://github.com/orbisgis/h2gis/pull/1312