search

unige-geohealth / accessmod

accessmod 5 : anisotropic accessibility analysis.
GNU Lesser General Public License v3.0
41 stars 14 forks source link

CRS of merged landcover layer is not recognized when exported #433

Closed RikLubbers closed 2 weeks ago

RikLubbers commented 6 months ago

Current Behavior

When exporting the merged landcover raster from AccessMod and then importing it into an R or GIS environment, the coordinate reference system (CRS) is not recognized. This issue can lead to problems in subsequent layering and spatial analyses due to the absence of correct spatial referencing.

Expected Behavior

Upon export from AccessMod and import into an R or GIS environment, the merged landcover raster should retain its CRS information correctly.

Possible Solution

This problem might be due to the use of deprecated PROJ.4 strings for CRS definitions (I observe that this string is provided in R). A potential solution could be to update AccessMod's handling of CRS.

Steps to Reproduce

  1. Export a merged landcover raster from AccessMod using the standard procedure.
  2. Import the exported raster file into an R environment or another GIS software.
  3. Notice that the CRS of the imported raster is not recognized or is incorrectly defined.

Detailed Description

The inability to recognize the CRS upon import into R or GIS environments indicates a possible issue with how AccessMod handles CRS information during the export of spatial data, potentially due to a reliance on outdated spatial reference definitions like PROJ.4 strings. Given that accurate spatial referencing is vital for geospatial analysis, particularly in fields where spatial alignment and precision are critical, addressing this issue is imperative. By updating the CRS handling mechanism to align with current standards and ensuring the correct export of spatial metadata, the utility of AccessMod for subsequent spatial analyses could be significantly enhanced.

fxi commented 5 months ago

Hi,

Sorry for the delay.

Could you provide the version of AccessMod? The version number should be something like AccessMod Desktop 5.8.2.

There are some outdated spatial components in AccessMod, but that would be the first time we encounter a blocking situation, which would be an urgent issue to solve indeed.

nicolasray commented 5 months ago

I can confirm that Rik was working with us when this issue arised, and that he used the latest version of AccessMod.

fxi commented 5 months ago

I just exported a merged landcover raster and the proj / wkt2 seems to be parsable, except the label. It's recognized in -gdal : ok -qgis : ok -R via terra : ok

╰$ gdalsrsinfo raster_land_cover_merged_crossnew.img 

PROJ.4 : +proj=utm +zone=35 +south +datum=WGS84 +units=m +no_defs

OGC WKT2:2019 :
PROJCRS["unknown",
    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["UTM zone 35S",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",27,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",10000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["meters",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["meters",1]],
    ID["EPSG",32735]]

Am I missing something ? Can you provide a landcover merged layer that's not recognized for examination ? If the file is large, please provide a https://www.swisstransfer.com/ link or similar.

RikLubbers commented 5 months ago

Sorry for the late reply, I was on holiday.

This is a .img file of a MLC directly exported from accessmod with PROJCRS[\"Arc 1960 / UTM zone 36N\",\n for all inputted files, prepared using the inAccessMod r package.

link

I load this .img file using terra::raster and when verifying the crs using terra::crs I get this output. My understanding is that subsequently problems are created further in the analysis due to the PROJCRS[\"unknown\",\n and DATUM[\"unknown\",\n part.

[1] "BOUNDCRS[\n    SOURCECRS[\n        PROJCRS[\"unknown\",\n            BASEGEOGCRS[\"GCS_clark80\",\n                DATUM[\"unknown\",\n                    ELLIPSOID[\"Clarke 1880\",6378249.145,293.465,\n                        LENGTHUNIT[\"metre\",1],\n                        ID[\"EPSG\",7034]]],\n                PRIMEM[\"Greenwich\",0,\n                    ANGLEUNIT[\"Degree\",0.0174532925199433]]],\n            CONVERSION[\"UTM zone 36N\",\n                METHOD[\"Transverse Mercator\",\n                    ID[\"EPSG\",9807]],\n                PARAMETER[\"Latitude of natural origin\",0,\n                    ANGLEUNIT[\"Degree\",0.0174532925199433],\n                    ID[\"EPSG\",8801]],\n                PARAMETER[\"Longitude of natural origin\",33,\n
  ANGLEUNIT[\"Degree\",0.0174532925199433],\n                    ID[\"EPSG\",8802]],\n                PARAMETER[\"Scale factor at natural origin\",0.9996,\n                    SCALEUNIT[\"unity\",1],\n                    ID[\"EPSG\",8805]],\n                PARAMETER[\"False easting\",500000,\n                    LENGTHUNIT[\"metre\",1],\n                    ID[\"EPSG\",8806]],\n
    PARAMETER[\"False northing\",0,\n                    LENGTHUNIT[\"metre\",1],\n                    ID[\"EPSG\",8807]],\n                ID[\"EPSG\",16036]],\n            CS[Cartesian,2],\n     
           AXIS[\"(E)\",east,\n                    ORDER[1],\n                    LENGTHUNIT[\"metre\",1,\n                        ID[\"EPSG\",9001]]],\n                AXIS[\"(N)\",north,\n       
             ORDER[2],\n                    LENGTHUNIT[\"metre\",1,\n                        ID[\"EPSG\",9001]]]]],\n    TARGETCRS[\n        GEOGCRS[\"WGS 84\",\n            ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n                MEMBER[\"World Geodetic System 1984 (Transit)\"],\n                MEMBER[\"World Geodetic System 1984 (G730)\"],\n                MEMBER[\"World Geodetic System 1984 (G873)\"],\n                MEMBER[\"World Geodetic System 1984 (G1150)\"],\n                MEMBER[\"World Geodetic System 1984 (G1674)\"],\n                MEMBER[\"World Geodetic System 1984 (G1762)\"],\n                MEMBER[\"World Geodetic System 1984 (G2139)\"],\n                ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                    LENGTHUNIT[\"metre\",1]],\n                ENSEMBLEACCURACY[2.0]],\n            PRIMEM[\"Greenwich\",0,\n                ANGLEUNIT[\"degree\",0.0174532925199433]],\n            CS[ellipsoidal,2],\n                AXIS[\"geodetic latitude (Lat)\",north,\n                    ORDER[1],\n                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n                AXIS[\"geodetic longitude (Lon)\",east,\n
     ORDER[2],\n                    ANGLEUNIT[\"degree\",0.0174532925199433]],\n            USAGE[\n                SCOPE[\"Horizontal component of 3D system.\"],\n                AREA[\"World.\"],\n                BBOX[-90,-180,90,180]],\n            ID[\"EPSG\",4326]]],\n    ABRIDGEDTRANSFORMATION[\"Transformation to WGS84\",\n        METHOD[\"Position Vector transformation (geog2D domain)\",\n            ID[\"EPSG\",9606]],\n        PARAMETER[\"X-axis translation\",0,\n            ID[\"EPSG\",8605]],\n        PARAMETER[\"Y-axis translation\",0,\n            ID[\"EPSG\",8606]],\n   
     PARAMETER[\"Z-axis translation\",0,\n            ID[\"EPSG\",8607]],\n        PARAMETER[\"X-axis rotation\",0,\n            ID[\"EPSG\",8608]],\n        PARAMETER[\"Y-axis rotation\",0,\n     
       ID[\"EPSG\",8609]],\n        PARAMETER[\"Z-axis rotation\",0,\n            ID[\"EPSG\",8610]],\n        PARAMETER[\"Scale difference\",1,\n            ID[\"EPSG\",8611]]]]"
fxi commented 5 months ago

I confirm that this file has no datum set in the SRS.

SRS info of the provided file :

-> % gdalsrsinfo  raster_land_cover_merged_MLC0_major_roads_10m.img

PROJ.4 : +proj=utm +zone=36 +a=6378249.145 +rf=293.465 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

OGC WKT2:2019 :
BOUNDCRS[
    SOURCECRS[
        PROJCRS["unknown",
            BASEGEOGCRS["GCS_clark80",
                DATUM["unknown",
                    ELLIPSOID["Clarke 1880",6378249.145,293.465,
                        LENGTHUNIT["metre",1],
                        ID["EPSG",7034]]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["Degree",0.0174532925199433]]],
            CONVERSION["UTM zone 36N",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["Latitude of natural origin",0,
                    ANGLEUNIT["Degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",33,
                    ANGLEUNIT["Degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9996,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",500000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]],
                ID["EPSG",16036]],
            CS[Cartesian,2],
                AXIS["(E)",east,
                    ORDER[1],
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]],
                AXIS["(N)",north,
                    ORDER[2],
                    LENGTHUNIT["metre",1,
                        ID["EPSG",9001]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            ENSEMBLE["World Geodetic System 1984 ensemble",
                MEMBER["World Geodetic System 1984 (Transit)"],
                MEMBER["World Geodetic System 1984 (G730)"],
                MEMBER["World Geodetic System 1984 (G873)"],
                MEMBER["World Geodetic System 1984 (G1150)"],
                MEMBER["World Geodetic System 1984 (G1674)"],
                MEMBER["World Geodetic System 1984 (G1762)"],
                MEMBER["World Geodetic System 1984 (G2139)"],
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ENSEMBLEACCURACY[2.0]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["geodetic latitude (Lat)",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["geodetic longitude (Lon)",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            USAGE[
                SCOPE["Horizontal component of 3D system."],
                AREA["World."],
                BBOX[-90,-180,90,180]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]
RikLubbers commented 5 months ago

Version that was used is Version[5.8.2-alpha.2]

Checking the CRS of the DEM and landcover that are exported from inAccessMod to AccessMod both have the following

PROJCRS["Arc 1960 / UTM zone 36N",
    BASEGEOGCRS["Arc 1960",
        DATUM["Arc 1960",
            ELLIPSOID["Clarke 1880 (RGS)",6378249.145,293.465,
 LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4210]],
    CONVERSION["UTM zone 36N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",33,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            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["Engineering survey, topographic mapping."],
        AREA["Kenya - north of equator and west of 36°E; Uganda - north of equator and east of 30°E."],
        BBOX[0,29.99,4.63,36]],
    ID["EPSG",21096]]

[edit : formating]

fxi commented 5 months ago

Ok. I don't have much time today, but it looks like some parameters are skipped after the export, and hopefully only at that stage, hopefully in an edge case... Comparing the WKT from the provided file (b) and the WKT you reported (a). https://colab.research.google.com/drive/16rpT0Ur85qsKrzraXLdGXjOuCELwljwB?usp=sharing

image
RikLubbers commented 5 months ago

Thanks for your input and I hope indeed it is an edge case. Luckily I found the mistake when verifying my data for further analysis, so it was fixed.

fxi commented 5 months ago

I've done a quick check on 13 projects, and all datums seem to be preserved and reported. The SRS seems to be correctly recognized in GIS tools. The issue reported, "CRS of merged landcover layer is not recognized when exported," doesn't seem to be reproducible with the current test data.

image

However, as mentioned, proj4 strings are no longer considered sufficient and could lead to mismatch in datums after export. Still unclear for me how/why. The whole import system in AccessMod is based on that, as it was the norm for a long time, and the export uses mostly default parameters.

https://github.com/unige-geohealth/accessmod/blob/3777ea4ac733e8ac55f4b8389813eb419718476d/tools/R/amExport.R#L231

AccessMod uses R and GRASS to import layers. AccessMod should not be directly involved, as it uses lower-level tools to do that. The underlying packages use standard Proj and GDAL libraries, which AccessMod has kept relatively up-to-date with GDAL 3.6.3 and Proj 9.2.0 (~March 2023).

https://github.com/unige-geohealth/accessmod/blob/3777ea4ac733e8ac55f4b8389813eb419718476d/docker/alpine_base/Dockerfile#L4

This issue could lead to a change in how GRASS projects (location & mapset) are set up, and this part could probably be improved, e.g., using the wkt (ref) option instead of proj4 in:

https://github.com/unige-geohealth/accessmod/blob/3777ea4ac733e8ac55f4b8389813eb419718476d/tools/R/amProjectImportExport.R#L297

More reading:

nicolasray commented 5 months ago

Re-opening the issue as it was not solved, and we'll follow up soon on this.

fxi commented 3 weeks ago

I've built a script/pipeline to fetch and build a project using inAccessMod and the landcover script for a country in one go. I then imported and exported the landcover data to test for projection issues, such as datum loss.

I used the latest version that will be available in 5.8.3-beta.0, (later this week) which includes the refactoring of all import/export spatial operations in AccessMod, the switch from proj.4 to WKT, and updates to the latest versions of all major dependencies.

It appears that the issue is resolved, according to this script:

Result

CRS Information for /content/drive/MyDrive/accessmod/ldc_merged_imported.img:
EPSG: 32632
Proj4: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

CRS Information for /content/drive/MyDrive/accessmod/ldc_merged_exported.tif:
EPSG: 32632
Proj4: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

Are the CRS the same?
Yes

Finished comparing CRS.
Python code ```python from osgeo import gdal from osgeo import osr def get_crs_info(file_path): ds = gdal.Open(file_path) if ds is None: return None, None, None wkt = ds.GetProjection() srs = osr.SpatialReference() srs.ImportFromWkt(wkt) srs.AutoIdentifyEPSG() epsg = srs.GetAuthorityCode(None) proj4 = srs.ExportToProj4() return wkt, epsg, proj4 # File paths file1 = '/content/drive/MyDrive/accessmod/ldc_merged_imported.img' file2 = '/content/drive/MyDrive/accessmod/ldc_merged_exported.tif' # Get CRS information for both files wkt1, epsg1, proj4_1 = get_crs_info(file1) wkt2, epsg2, proj4_2 = get_crs_info(file2) # Print results print(f"CRS Information for {file1}:") print(f"EPSG: {epsg1}") print(f"Proj4: {proj4_1}") print(f"\nCRS Information for {file2}:") print(f"EPSG: {epsg2}") print(f"Proj4: {proj4_2}") # Compare CRS if wkt1 and wkt2: srs1 = osr.SpatialReference() srs2 = osr.SpatialReference() srs1.ImportFromWkt(wkt1) srs2.ImportFromWkt(wkt2) print("\nAre the CRS the same?") print("Yes" if srs1.IsSame(srs2) == 1 else "No") else: print("\nUnable to compare CRS due to error in reading one or both files.") print("\nFinished comparing CRS.") ```
RikLubbers commented 2 weeks ago

I'm pleased to report that the fix is now working. After reprocessing all the raw files (DEM, population raster) again using the inAccessMod package in R and uploading them into a new AccessMod project, the newly created MLC retains the assigned CRS when exporting.

Old files still lack the CRS information, but that's expected since they were processed before the fix was applied. It appears the problem is now solved for new projects. Thank you for addressing this issue!