1. FAO-SDLC Avatar FAO-SDLC
  2. WaPOR ET Look
  3. pywapor
  4. Issues

Georeferencing issues with output

Issue #6 resolved
Claire Michailovsky created an issue

The output file et_look_out.nc does not have a coordinate reference system.

gdal_info on the output file is incorrect.

Opening in QGIS, one layer of the output is North South flipped.

‌

Comments (15)

  2. bert.coerver
    • changed status to open

    Same for "se_root_in.nc", "et_look_in.nc" and "et_look_out.nc".

    "se_root_out.nc" does have CRS.

  3. Claire Michailovsky reporter

    Update: It appears that when run from Colab, the .nc file opens properly in QGIS.

  4. bert.coerver

    ok thanks, thats useful to know. I’m checking and it looks like it's maybe a bug in (rio)xarray and so Colab perhaps uses an older version.

  5. bert.coerver

    Hi Claire,

    I hadnt checked actually opening the files in QGIS myself, only checked gdalinfo on the file and found the same as you did.

    However, for netcdfs with subdatasets, you need to check the gdalinfo of the subdataset, not of the overarching netcdf file. E.g. when I run:

    gdalinfo 'NETCDF:"/Users/hmcoerver/Local/bich_pywapor_test/et_look_out.nc":int_mm'

    I get the info shown below. The file also opens correctly in QGIS:

    So I can’t reproduce this issue, you also mentioned that files coming from Colab are working. So can you maybe double check if this really an issue? Or provide me with more info to reproduce.

    ‌

    GDALINFO output:

    Driver: netCDF/Network Common Data Format
Files: /Users/hmcoerver/Local/bich_pywapor_test/et_look_out.nc
Size is 1709, 1368
Coordinate System is:
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["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (30.199910779347409,29.700531598611771)
Pixel Size = (0.000585351719396,-0.000585351719396)
Metadata:
  int_mm#calculated_with={p_24,vc,lai}
  int_mm#et_look_module=pywapor.et_look_v2_v3.evapotranspiration
  int_mm#grid_mapping=spatial_ref
  int_mm#long_name=Interception
  int_mm#scale_factor=1e-08
  int_mm#source_data={ndvi,p_24}
  int_mm#units=mm day-1
  int_mm#_FillValue=-9999
  NETCDF_DIM_EXTRA={time_bins}
  NETCDF_DIM_time_bins_DEF={3,6}
  NETCDF_DIM_time_bins_VALUES={0,1,2}
  spatial_ref#crs_wkt=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  spatial_ref#geographic_crs_name=WGS 84
  spatial_ref#GeoTransform=30.19991077934741 0.0005853517193961209 0.0 29.70053159861177 0.0 -0.0005853517193961205
  spatial_ref#grid_mapping_name=latitude_longitude
  spatial_ref#horizontal_datum_name=World Geodetic System 1984
  spatial_ref#inverse_flattening=298.257223563
  spatial_ref#longitude_of_prime_meridian=0
  spatial_ref#prime_meridian_name=Greenwich
  spatial_ref#reference_ellipsoid_name=WGS 84
  spatial_ref#semi_major_axis=6378137
  spatial_ref#semi_minor_axis=6356752.314245179
  spatial_ref#spatial_ref=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]
  time_bins#calendar=proleptic_gregorian
  time_bins#units=days since 2021-07-01 00:00:00
  time_bins#_FillValue=nan
  x#_FillValue=nan
  y#_FillValue=nan
Corner Coordinates:
Upper Left  (  30.1999108,  29.7005316) ( 30d11'59.68"E, 29d42' 1.91"N)
Lower Left  (  30.1999108,  28.8997704) ( 30d11'59.68"E, 28d53'59.17"N)
Upper Right (  31.2002769,  29.7005316) ( 31d12' 1.00"E, 29d42' 1.91"N)
Lower Right (  31.2002769,  28.8997704) ( 31d12' 1.00"E, 28d53'59.17"N)
Center      (  30.7000938,  29.3001510) ( 30d42' 0.34"E, 29d18' 0.54"N)
Band 1 Block=500x500 Type=Int64, ColorInterp=Undefined
  NoData Value=-9999
  Unit Type: mm day-1
  Offset: 0,   Scale:1e-08
  Metadata:
    calculated_with={p_24,vc,lai}
    et_look_module=pywapor.et_look_v2_v3.evapotranspiration
    grid_mapping=spatial_ref
    long_name=Interception
    NETCDF_DIM_time_bins=0
    NETCDF_VARNAME=int_mm
    scale_factor=1e-08
    source_data={ndvi,p_24}
    units=mm day-1
    _FillValue=-9999
Band 2 Block=500x500 Type=Int64, ColorInterp=Undefined
  NoData Value=-9999
  Unit Type: mm day-1
  Offset: 0,   Scale:1e-08
  Metadata:
    calculated_with={p_24,vc,lai}
    et_look_module=pywapor.et_look_v2_v3.evapotranspiration
    grid_mapping=spatial_ref
    long_name=Interception
    NETCDF_DIM_time_bins=1
    NETCDF_VARNAME=int_mm
    scale_factor=1e-08
    source_data={ndvi,p_24}
    units=mm day-1
    _FillValue=-9999
Band 3 Block=500x500 Type=Int64, ColorInterp=Undefined
  NoData Value=-9999
  Unit Type: mm day-1
  Offset: 0,   Scale:1e-08
  Metadata:
    calculated_with={p_24,vc,lai}
    et_look_module=pywapor.et_look_v2_v3.evapotranspiration
    grid_mapping=spatial_ref
    long_name=Interception
    NETCDF_DIM_time_bins=2
    NETCDF_VARNAME=int_mm
    scale_factor=1e-08
    source_data={ndvi,p_24}
    units=mm day-1
    _FillValue=-9999

    ‌

  7. Claire Michailovsky reporter

    Export to GTiff with gdal_translate works fine for individual variables/bands.

    Could be a QGIS issue - unsure… this is what I get (Bing maps in background).

    ‌

  8. bert.coerver
    1. So when check gdalinfo on the subdataset you do see the CRS etc?
    2. When you translate from netcdf to geotiff the files open correctly in QGIS?

    ‌

  10. bert.coerver

    ok if gdal (and xarray as well) handles the data correctly and shows the CRS, then it looks like this is a problem with QGIS. You could still try to compare the file you got from Colab with the one that doesnt work and see if you can find any difference.

  11. Claire Michailovsky reporter

    Sorry Bert, I had used the wrong output for gdal_translate. The one for which I put the image above in fact does not work well with gdal_translate, and the gdalinfo is also wrong. See below.

    ‌

  13. Bich Tran

    As Bert said, it’s a bug in rioxarray. I produced output using rioxarray version 0.15.0 before. After updating to rioxarray 0.15.5, the output is ok (Tested on both Windows and Linux).

    Here’s the gdalinfo of the outputs using old rioxarray version and new rioxarray version. The old version dataset does not have CRS information:

    ‌

    Suggestion: add requirement of xarray version >=0.15.5 in setup.py, so that for people with old version of xarray already installed, conda will have to upgrade it when installing pywapor.

    ‌

  16. Log in to comment
Assignee
–
Type
bug
Priority
major
Status
resolved
Votes
2
Watchers
3
Jira: the preferred issue tracker for Bitbucket. Join the team!