qgis / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
https://qgis.org
GNU General Public License v2.0
10.52k stars 2.99k forks source link

netcdf file (.nc4) added with latitude and longitude swapped #33082

Closed rusanpe closed 4 years ago

rusanpe commented 4 years ago

Dear Friends,

I have been looking trough Google about this problem but I didn't found any solution. When I add a nc4 file, Qgis assign it the WGS84 crs (EPSG:4326). For some reason, the file has the latitude coordinates first, and Qgis swap the coordinates, so the image appears over the antartida (-6,-74) instead of over Perú, (-74,-6) There is no any way to change it. When you create a WFS connection you can mark "reverse axis orientation" for example. Something like this would fix this problem.

Thanks a lot

gioman commented 4 years ago

@rusanpe QGIS version? platform? sample data?

rusanpe commented 4 years ago

@gioman, sorry, this was my first contribution. I've try to add nc4 files in 3.6, 3.8 and 3.10 versión. Even 2.x version, with plugin "netcdf browser" add the file with the latitude and longitude swapped.

I've try to upload a file as example but I didn't found any option to do it

Thanks

gioman commented 4 years ago

I've try to upload a file as example but I didn't found any option to do it

@rusanpe you can attach files here, zip it before trying do so.

rusanpe commented 4 years ago

3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.zip Sorry, this is the sample file

gioman commented 4 years ago

I've try to add nc4 files in 3.6, 3.8 and 3.10 versión

@rusanpe on what platform?

rusanpe commented 4 years ago

Windows 10 with the 3.x version, Linux (2.x version with osgeo-live version 11 and 3.4 version in osgeo-live version 13)

gioman commented 4 years ago

For some reason, the file has the latitude coordinates first, and Qgis swap the coordinates, so the imag

gdalinfo 3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4 
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: netCDF/Network Common Data Format
Files: 3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4
Size is 142, 170
Coordinate System is `'
**Origin = (-12.199999617346634,-66.299998469606663)**
Pixel Size = (0.099999997632723,-0.100000009028869)
Metadata:
  lat#axis=Y
  lat#bounds=lat_bnds
  lat#DimensionNames=lat
  lat#fullnamepath=/Grid/lat
  lat#LongName=Latitude at the center of
                        0.10 degree grid intervals of latitude
                        from -90 to 90.
  lat#origname=lat
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#axis=X
  lon#bounds=lon_bnds
  lon#DimensionNames=lon
  lon#fullnamepath=/Grid/lon
  lon#LongName=Longitude at the center of
                        0.10 degree grid intervals of longitude                                                 
                        from -180 to 180.                                                                       
  lon#origname=lon
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#DODS_EXTRA.Unlimited_Dimension=Grid_time
  NC_GLOBAL#FileHeader=DOI=10.5067/GPM/IMERG/3B-HH/06;
DOIauthority=http://dx.doi.org/;
DOIshortName=3IMERGHH;
AlgorithmID=3IMERGHH;
AlgorithmVersion=3IMERGH_6.3;
FileName=3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5;
SatelliteName=MULTI;
InstrumentName=MERGED;
GenerationDateTime=2019-06-25T16:04:14.000Z;
StartGranuleDateTime=2008-03-26T00:00:00.000Z;
StopGranuleDateTime=2008-03-26T00:29:59.999Z;
GranuleNumber=;
NumberOfSwaths=0;
NumberOfGrids=1;
GranuleStart=;
TimeInterval=HALF_HOUR;
ProcessingSystem=PPS;
ProductVersion=V06B;
EmptyGranule=NOT_EMPTY;
MissingData=;

  NC_GLOBAL#FileInfo=DataFormatVersion=6a;
TKCodeBuildVersion=0;
MetadataVersion=6a;
FormatPackage=HDF5-1.8.9;
BlueprintFilename=GPM.V6.3IMERGHH.blueprint.xml;
BlueprintVersion=BV_62;
TKIOVersion=3.93;
MetadataStyle=PVL;
EndianType=LITTLE_ENDIAN;

  NC_GLOBAL#Grid.fullnamepath=/Grid
  NC_GLOBAL#Grid.GridHeader=BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.1;
LongitudeResolution=0.1;
NorthBoundingCoordinate=90;
SouthBoundingCoordinate=-90;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;

  NC_GLOBAL#history=2019-09-30 09:18:18 GMT Hyrax-1.15.1 https://gpm1.gesdisc.eosdis.nasa.gov:443/opendap/GPM_L3/GPM_3IMERGHH.06/2008/086/3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4?precipitationCal[0:0][967:1136][778:919],time,lon[967:1136],lat[778:919]
  NETCDF_DIM_EXTRA={time}
  NETCDF_DIM_time_DEF={1,4}
  NETCDF_DIM_time_VALUES=1206489600
  precipitationCal#CodeMissingValue=-9999.9
  precipitationCal#coordinates=time lon lat
  precipitationCal#DimensionNames=time,lon,lat
  precipitationCal#fullnamepath=/Grid/precipitationCal
  precipitationCal#origname=precipitationCal
  precipitationCal#units=mm/hr
  precipitationCal#_FillValue=-9999.9004
  time#axis=T
  time#bounds=time_bnds
  time#calendar=julian
  time#DimensionNames=time
  time#fullnamepath=/Grid/time
  time#LongName=Representative time of data in 
                        seconds since 1970-01-01 00:00:00 UTC.
  time#origname=time
  time#standard_name=time
  time#units=seconds since 1970-01-01 00:00:00 UTC
Subdatasets:
  SUBDATASET_1_NAME=NETCDF:"3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4":precipitationCal
  SUBDATASET_1_DESC=[1x170x142] precipitationCal (32-bit floating-point)
Geolocation:
  LINE_OFFSET=0
  LINE_STEP=1
  PIXEL_OFFSET=0
  PIXEL_STEP=1
  SRS=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"]],AUTHORITY["EPSG","4326"]]
  X_BAND=1
  X_DATASET=NETCDF:"3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4":lon
  Y_BAND=1
  Y_DATASET=NETCDF:"3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4":lat
Corner Coordinates:
Upper Left  ( -12.1999996, -66.2999985) 
Lower Left  ( -12.1999996, -83.3000000) 
Upper Right (   2.0000000, -66.2999985) 
Lower Right (   2.0000000, -83.3000000) 
Center      (  -5.0999998, -74.7999992) 
Band 1 Block=142x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999.900390625
  Unit Type: mm/hr
  Metadata:
    CodeMissingValue=-9999.9
    coordinates=time lon lat
    DimensionNames=time,lon,lat
    fullnamepath=/Grid/precipitationCal
    NETCDF_DIM_time=1206489600
    NETCDF_VARNAME=precipitationCal
    origname=precipitationCal
    units=mm/hr
    _FillValue=-9999.9004
rusanpe commented 4 years ago

@gioman I've already seen this metadata. ArcGIS, through the "Make a raster layer from a netCDF file" command creates a correct raster image. Qgis don't need any command, but it's not be able to project the same file correctly.

roya0045 commented 4 years ago

QGIS did not swap the image if the coordinates are using inverse notation rather than the standard one. I guess that relying on the label might be risky rather than the order.

saberraz commented 4 years ago

@rusanpe have you tried to add it as a Mesh layer: https://docs.qgis.org/3.4/en/docs/user_manual/working_with_mesh/mesh_properties.html#loading-a-mesh-layer

rusanpe commented 4 years ago

@saberraz, thanks, the result is the same.

saberraz commented 4 years ago

@rusanpe Could you try this grib which I have converted from your nc4 file?

test.grb.zip

rusanpe commented 4 years ago

@saberraz It Works! Have you created a new specific SRC projection for this file?

saberraz commented 4 years ago

TLDR You need to run this command: cdo -v -f grb -copy 3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4 test.grb

(CDO can be installed in ubuntu by sudo apt-get install cdo For Windows and other platform see: https://www.isimip.org/protocol/isimip2b-files/cdo-help/).

Longer version:

The problem with your file is that precipitation variable is saved in this order: time, lon, lat. You can run ncdump -h 3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4

dimensions:
        time = 1 ;
        lon = 170 ;
        lat = 142 ;
variables:
        float precipitationCal(time, lon, lat) ;
                precipitationCal:DimensionNames = "time,lon,lat" ;
                precipitationCal:Units = "mm/hr" ;
                precipitationCal:units = "mm/hr" ;
                precipitationCal:coordinates = "time lon lat" ;
                precipitationCal:_FillValue = -9999.9f ;
                precipitationCal:CodeMissingValue = "-9999.9" ;
                precipitationCal:origname = "precipitationCal" ;
                precipitationCal:fullnamepath = "/Grid/precipitationCal" ;
        float lat(lat) ;
                lat:DimensionNames = "lat" ;
                lat:Units = "degrees_north" ;
                lat:units = "degrees_north" ;
                lat:standard_name = "latitude" ;
                lat:LongName = "Latitude at the center of\n",
                        "\t\t\t0.10 degree grid intervals of latitude\n",
                        "\t\t\tfrom -90 to 90." ;
                lat:bounds = "lat_bnds" ;
                lat:axis = "Y" ;
                lat:origname = "lat" ;
                lat:fullnamepath = "/Grid/lat" ;
        float lon(lon) ;
                lon:DimensionNames = "lon" ;
                lon:Units = "degrees_east" ;
                lon:units = "degrees_east" ;
                lon:standard_name = "longitude" ;
                lon:LongName = "Longitude at the center of\n",
                        "\t\t\t0.10 degree grid intervals of longitude \n",
                        "\t\t\tfrom -180 to 180." ;
                lon:bounds = "lon_bnds" ;
                lon:axis = "X" ;
                lon:origname = "lon" ;
                lon:fullnamepath = "/Grid/lon" ;
        int time(time) ;
                time:DimensionNames = "time" ;
                time:Units = "seconds since 1970-01-01 00:00:00 UTC" ;
                time:units = "seconds since 1970-01-01 00:00:00 UTC" ;
                time:standard_name = "time" ;
                time:LongName = "Representative time of data in \n",
                        "\t\t\tseconds since 1970-01-01 00:00:00 UTC." ;
                time:bounds = "time_bnds" ;
                time:axis = "T" ;
                time:calendar = "julian" ;
                time:origname = "time" ;
                time:fullnamepath = "/Grid/time" ;

// global attributes:
                :FileHeader = "DOI=10.5067/GPM/IMERG/3B-HH/06;\n",
                        "DOIauthority=http://dx.doi.org/;\n",
                        "DOIshortName=3IMERGHH;\n",
                        "AlgorithmID=3IMERGHH;\n",
                        "AlgorithmVersion=3IMERGH_6.3;\n",
                        "FileName=3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5;\n",
                        "SatelliteName=MULTI;\n",
                        "InstrumentName=MERGED;\n",
                        "GenerationDateTime=2019-06-25T16:04:14.000Z;\n",
                        "StartGranuleDateTime=2008-03-26T00:00:00.000Z;\n",
                        "StopGranuleDateTime=2008-03-26T00:29:59.999Z;\n",
                        "GranuleNumber=;\n",
                        "NumberOfSwaths=0;\n",
                        "NumberOfGrids=1;\n",
                        "GranuleStart=;\n",
                        "TimeInterval=HALF_HOUR;\n",
                        "ProcessingSystem=PPS;\n",
                        "ProductVersion=V06B;\n",
                        "EmptyGranule=NOT_EMPTY;\n",
                        "MissingData=;\n",
                        "" ;
                :FileInfo = "DataFormatVersion=6a;\n",
                        "TKCodeBuildVersion=0;\n",
                        "MetadataVersion=6a;\n",
                        "FormatPackage=HDF5-1.8.9;\n",
                        "BlueprintFilename=GPM.V6.3IMERGHH.blueprint.xml;\n",
                        "BlueprintVersion=BV_62;\n",
                        "TKIOVersion=3.93;\n",
                        "MetadataStyle=PVL;\n",
                        "EndianType=LITTLE_ENDIAN;\n",
                        "" ;
                :Grid.GridHeader = "BinMethod=ARITHMETIC_MEAN;\n",
                        "Registration=CENTER;\n",
                        "LatitudeResolution=0.1;\n",
                        "LongitudeResolution=0.1;\n",
                        "NorthBoundingCoordinate=90;\n",
                        "SouthBoundingCoordinate=-90;\n",
                        "EastBoundingCoordinate=180;\n",
                        "WestBoundingCoordinate=-180;\n",
                        "Origin=SOUTHWEST;\n",
                        "" ;
                :Grid.fullnamepath = "/Grid" ;
                :DODS_EXTRA.Unlimited_Dimension = "Grid_time" ;
                :history = "2019-09-30 09:18:18 GMT Hyrax-1.15.1 https://gpm1.gesdisc.eosdis.nasa.gov:443/opendap/GPM_L3/GPM_3IMERGHH.06/2008/086/3B-HHR.MS.MRG.3IMERG.20080326-S000000-E002959.0000.V06B.HDF5.nc4?precipitationCal[0:0][967:1136][778:919],time,lon[967:1136],lat[778:919]" ;
}

There are some tricks to re-order the variables and have it in lon, lat, time. But I just converted it to GRIB and it did the job.

I suggest file a ticket here: https://github.com/lutraconsulting/MDAL/issues

For GRIB and Netcdf files (especially when they have time dimension), it is better to load them as mesh (which uses MDAL...see my earlier comment about how to work with mesh).

I close this one, but feel free to reference to this ticket in your other issue report.