Need help adding a CRS to a Raster

[biskwikman]

Hello. I've got a Raster from a Matrix, which i read from a binary file. It should represent the entire Earth. It's 86400 by 43200 pixels (x and y). It is unprojected data. I just can't figure out how to get Rasters to attach a CRS to it. I know there is a crs keyword in the Raster function, but I don't really understand how to use it. I tried reading the documentation for GeoFormatTypes.jl but it's pretty sparse, and I didn't come away understanding how to use it with Rasters for the CRS. Can someone help me?

Right now I have this:

Raster(myarray; dims=(X, Y), missingval=32767)

[lazarusA]

maybe this: https://rafaqz.github.io/Rasters.jl/stable/tutorials/resample#A-Raster-from-scratch

[asinghvi17]

Yes, crs=EPSG(4326) needs to be a keyword argument to the Raster constructor. But that is not all - you also need to assign the correct X and Y dimensions such that this is actually in the correct long/lat space.

julia> using Rasters, Rasters.Lookups

julia> myarray = rand(864, 432);

julia> sx, sy = size(myarray)
(864, 432)

julia> xdim = X(Sampled(LinRange(-180, 180, sx+1)[1:end-1]; sampling = Intervals(Start())))
X Sampled{Float64} AutoOrder AutoSpan Intervals{Start}
wrapping: 864-element LinRange{Float64, Int64}:
 -180.0, -179.583, -179.167, -178.75, -178.333, -177.917, -177.5, -177.083, -176.667, -176.25, -175.833, …, 175.0, 175.417, 175.833, 176.25, 176.667, 177.083, 177.5, 177.917, 178.333, 178.75, 179.167, 179.583

julia> ydim = Y(Sampled(LinRange(-90, 90, sy+1)[1:end-1]; sampling = Intervals(Start())))
Y Sampled{Float64} AutoOrder AutoSpan Intervals{Start}
wrapping: 432-element LinRange{Float64, Int64}:
 -90.0, -89.5833, -89.1667, -88.75, -88.3333, -87.9167, -87.5, -87.0833, -86.6667, -86.25, -85.8333, -85.4167, …, 85.0, 85.4167, 85.8333, 86.25, 86.6667, 87.0833, 87.5, 87.9167, 88.3333, 88.75, 89.1667, 89.5833

julia> ras = Raster(myarray; dims = (xdim, ydim), missingval = 32767, crs = EPSG(4326))
┌ 864×432 Raster{Float64, 2} ┐
├────────────────────────────┴──────────────────────────────────────────────── dims ┐
  ↓ X Projected{Float64} LinRange{Float64}(-180.0, 179.58333333333331, 864) ForwardOrdered Regular Intervals{Start},
  → Y Projected{Float64} LinRange{Float64}(-90.0, 89.58333333333334, 432) ForwardOrdered Regular Intervals{Start}
├─────────────────────────────────────────────────────────────────────────── raster ┤
  missingval: 32767.0
  extent: Extent(X = (-180.0, 179.99999999999997), Y = (-90.0, 90.00000000000001))
  crs: EPSG:4326
└───────────────────────────────────────────────────────────────────────────────────┘
    ↓ →    -90.0         -89.5833     …  88.75       89.1667    89.5833
 -180.0      0.64777       0.861928       0.507233    0.656124   0.359831
 -179.583    0.00720807    0.865115       0.225164    0.35233    0.485258
 -179.167    0.0983575     0.657511       0.0055762   0.435763   0.931783
    ⋮                                 ⋱   ⋮
  178.75     0.584729      0.943945       0.791519    0.60058    0.754317
  179.167    0.510269      0.0838822      0.210341    0.741584   0.858012
  179.583    0.335859      0.641238       0.806762    0.877428   0.0742643

This is the "correct" way to construct a raster whose cells are treated as areas (intervals) rather than points, which spans the globe. @rafaqz this is really clumsy though - we should make this easier to do!

[rafaqz]

There already are much cleaner ways?

You can pass the sampling keyword to the Dimensions constructor! Yes that's kind of magic but if does exactly what you did above.

E.g. just Y(1:100; sampling=Intervals()))

Its not documented but Rasters.create let's you pass the sampling keyword for all dimensions.

[biskwikman]

Thanks for your advice everyone! The response by lazarusA made the most sense to me so that's what I used and it worked!