crs not parsed from NetCDF

[alex-s-gardner]

When reading in a netcdf it seems that the crs is not parsed from the mapping_data... I haven't had this issue with this file and other geospatial software so I'm wondering if there is a bug somewhere.

In this example Rasters.crs() returns nothing

paths2granule = "https://hyp3-jhk-sandbox-contentbucket-wxvu5zy6v7d6.s3.us-west-2.amazonaws.com/100f0e79-d569-42b9-8a20-3c3384b34dee/S1B_IW_SLC__1SSH_20181227T194309_20181227T194320_014232_01A75C_12D4_X_S1A_IW_SLC__1SSH_20190102T194350_20190102T194401_025303_02CCA3_D38F_G0120V02_P099.nc"

local_path = joinpath(tempdir(), splitpath(paths2granule)[end])

isfile(local_path) || Downloads.download(paths2granule, local_path)

rs = RasterStack(local_path)

Rasters.crs(rs[:v])

[asinghvi17]

Rasters does not support parsing CF compliant CRS yet. It's something we were looking to do even when I was at JPL but never got around to.

You can force reading via source = Rasters.GDALsource() to use GDAL which will parse the CRS correctly...

[alex-s-gardner]

I've got the memory of a fish!

Thanks for the GDAL tip!

I see this has been open for a number of years now: https://github.com/rafaqz/Rasters.jl/issues/31 but I'll leave this open as it might be helpful for others to see the GDAL suggestion.

[asinghvi17]

We should at least detect WKT...

[asinghvi17]

 mapping
    Attributes:
     grid_mapping_name    = polar_stereographic
     straight_vertical_longitude_from_pole = -45.0
     false_easting        = 0.0
     false_northing       = 0.0
     latitude_of_projection_origin = 90.0
     latitude_of_origin   = 70.0
     semi_major_axis      = 6.378137e6
     scale_factor_at_projection_origin = 1
     inverse_flattening   = 298.257223563
     spatial_ref          = PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",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"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH],AUTHORITY["EPSG","3413"]]
     crs_wkt              = PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",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"]],PROJECTION["Polar_Stereographic"],PARAMETER["latitude_of_origin",70],PARAMETER["central_meridian",-45],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",SOUTH],AXIS["Northing",SOUTH],AUTHORITY["EPSG","3413"]]
     proj4text            = +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
     spatial_epsg         = 3413
     GeoTransform         = -60787.5 120.0 0 -1340932.5 0 -120.0

[asinghvi17]

it was https://github.com/rafaqz/Rasters.jl/pull/803/ where I edited functionality that is associated with this, and https://github.com/rafaqz/Rasters.jl/tree/netcdf_crs

has some translation tools to go from CF to a projstring.

[asinghvi17]

given that this has spatial_ref, crs_wkt, proj4text and spatial_epsg fields we ought to at least check for those and assign potentially

at least crs_wkt, proj4text and spatial_epsg since we know what types those are

[asinghvi17]

https://github.com/rafaqz/Rasters.jl/blob/0a32706ca39bcaf6dc895fcbabac4af08991a2d1/src/sources/commondatamodel.jl#L296-L341 is where we create lookups for cf compliant stuff

[asinghvi17]

PR is up, @alex-s-gardner can you try that?

[alex-s-gardner]

Something broke between v0.13.0 and v0.14.0:

v0.13.0 loads the netcdf file fine, v0.14.0 throws this error so I can't test your new branch:

ERROR: Argument is an abstract type and does not have a definite size.
Stacktrace:
  [1] sizeof(x::Type)
    @ Base ./essentials.jl:780
  [2] _sizeof(A::Rasters.ModifiedDiskArray{false, Union{…}, 0, NCDatasets.Variable{…}, Rasters.Mod{…}})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/utils.jl:395
  [3] _checkobjmem(f::Function, obj::Rasters.ModifiedDiskArray{false, Union{…}, 0, NCDatasets.Variable{…}, Rasters.Mod{…}})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/utils.jl:391
  [4] _checkobjmem(obj::Rasters.ModifiedDiskArray{false, Union{…}, 0, NCDatasets.Variable{…}, Rasters.Mod{…}})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/utils.jl:389
  [5] #134
    @ ~/.julia/packages/Rasters/TYVYz/src/stack.jl:561 [inlined]
  [6] (::Base.var"#4#5"{…})(a::Tuple{…})
    @ Base ./generator.jl:37
  [7] iterate
    @ ./generator.jl:48 [inlined]
  [8] collect(itr::Base.Generator{Base.Iterators.Zip{Tuple{…}}, Base.var"#4#5"{Rasters.var"#134#139"{…}}})
    @ Base ./array.jl:780
  [9] map(::Function, ::Vector{…}, ::Vector{…}, ::Vector{…})
    @ Base ./abstractarray.jl:3495
 [10] (::Rasters.var"#130#135"{…})(ds::NCDatasets.NCDataset{…})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/stack.jl:559
 [11] #_open#1039
    @ ~/.julia/packages/Rasters/TYVYz/src/sources/commondatamodel.jl:68 [inlined]
 [12] _open
    @ ~/.julia/packages/Rasters/TYVYz/src/sources/commondatamodel.jl:55 [inlined]
 [13] #6
    @ ~/.julia/packages/Rasters/TYVYz/ext/RastersNCDatasetsExt/ncdatasets_source.jl:66 [inlined]
 [14] NCDatasets.NCDataset(::RastersNCDatasetsExt.var"#6#7"{…}, ::String, ::Vararg{…}; kwargs::@Kwargs{})
    @ NCDatasets ~/.julia/packages/NCDatasets/ydsnA/src/dataset.jl:252
 [15] NCDataset
    @ ~/.julia/packages/NCDatasets/ydsnA/src/dataset.jl:249 [inlined]
 [16] #_open#5
    @ ~/.julia/packages/Rasters/TYVYz/ext/RastersNCDatasetsExt/ncdatasets_source.jl:65 [inlined]
 [17] _open
    @ ~/.julia/packages/Rasters/TYVYz/ext/RastersNCDatasetsExt/ncdatasets_source.jl:62 [inlined]
 [18] #_open#21
    @ ~/.julia/packages/Rasters/TYVYz/src/sources/sources.jl:86 [inlined]
 [19] _open
    @ ~/.julia/packages/Rasters/TYVYz/src/sources/sources.jl:85 [inlined]
 [20] _layer_stack(filename::String; source::Rasters.NCDsource, dims::Rasters.NoKW, refdims::Tuple{}, name::Rasters.NoKW, group::Rasters.NoKW, metadata::Rasters.NoKW, layermetadata::Rasters.NoKW, layerdims::Rasters.NoKW, missingval::Rasters.NoKW, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW, coerce::Rasters.NoKW, scaled::Bool, checkmem::Bool, lazy::Bool, kw::@Kwargs{})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/stack.jl:526
 [21] RasterStack(filename::String; lazy::Bool, dropband::Bool, raw::Bool, source::Rasters.NoKW, missingval::Rasters.NoKW, name::Rasters.NoKW, group::Rasters.NoKW, scaled::Rasters.NoKW, coerce::Rasters.NoKW, verbose::Bool, replace_missing::Rasters.NoKW, kw::@Kwargs{})
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/stack.jl:435
 [22] RasterStack(filename::String)
    @ Rasters ~/.julia/packages/Rasters/TYVYz/src/stack.jl:399
 [23] top-level scope
    @ ~/Documents/GitHub/Altim.jl/src/junk1.jl:12
Some type information was truncated. Use `show(err)` to see complete types.

[asinghvi17]

Raf is on it

[asinghvi17]

But you can load a single layer using Raster for now

[rafaqz]

Yeah it's the layer that just has a String is broken. Will push a fix tonight