(help request) using ArchGDAL.gdalwarp() for resampling a raster dataset

Specific Domains Geo
,
1

Hi all, I’m struggling to figure out how to use ArchGDAL.gdalwarp (currently only unsafe_gdalwarp is documented). ArchGDAL.gdalwarp’s first parameter being a function is a bit confusing as well.

Here’s the GDAL CLI version of what I’m trying to implement in ArchGDAL (params I already have in Julia are represented by the all caps strings):

gdalwarp -t_srs WKT_STRING -te LEFT_EXTENT LOWER_EXTENT RIGHT_EXTENT \
    UPPER_EXTENT -te_srs WKT_STRING -tr X_RES Y_RES -r cubic \
    GDAL_DATASET output_dataset

I’m also hoping to return an ArchGDAL.Dataset if possible (without writing anything to disk).

Thanks for any help!

2

Okay! I think I figured it out…?

ds = ArchGDAL.read(path) do source
    ArchGDAL.gdalwarp([source], ["-t_srs", "$(wkt)",
                               "-te", "$(west_bound)", "$(south_bound)",
                               "$(east_bound)", "$(north_bound)",
                               "-te_srs", "$(wkt)",
                               "-tr", "$(step(x_dim))", "$(step(y_dim))",
                               "-r", "cubic"]) do warped
        band = ArchGDAL.getband(warped, 1)

        ArchGDAL.read(band)
    end
end
2 Likes
3

We may as well wrap this in GeoData.jl, that could be a one liner

4

I can do a PR if you’d like. Feel free to create an issue in GeoData.jl and tag me! Can probably get to that next week. Also feel free to copy/paste the above if you want to implement it.

5

Done: https://github.com/rafaqz/GeoData.jl/issues/68

2 Likes
6

I have a similar issue trying to resample with gdalwarp(). I want to write a function to save both the warped dataset and the Matrix of the band for processing. I need the dataset to use the getgeotransform(dataset) and getproj(dataset) methods to retrieve the x/y from the row/col of the pixel.

From my understanding, it seems like gdalwarp() calls the unsafe_gdalwarp() method and then destroys the warped dataset. This is what I tried, but the returned dataset is NULL

function reproject_dem(param, output_bounds)
    dst_cs = AG.importEPSG(point2epsg(param.lon, param.lat)) #  UTM projected coordinate system
    AG.read(param.dem_file) do source 
        AG.gdalwarp([source], ["-t_srs", "$(AG.toWKT(dst_cs))",
                        "-te",string.(output_bounds)...,
                        "-te_srs", "$(AG.toWKT(dst_cs))",
                        "-tr","$(param.pixel_spacing_in_m)","$(param.pixel_spacing_in_m)",
                        "-r", "cubic"]) do warped
            warped,
            AG.read(AG.getband(warped, 1))
        end
    end
end

Thanks for any help!

7

Rasters.jl wraps AG.gdalwarp for easy warping and resampling. The resample method is pretty straightforward compared to this.

8

And with GMT you can do things like

using GMT

julia> G = GMT.peaks();

julia> grdinfo(G)
11-element Vector{String}:
 ": Title: "
 ": Command: "
 ": Remark: "
 ": Gridline node registration used [Cartesian grid]"
 ": Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7"
 ": x_min: -3 x_max: 3 x_inc: 0.125 name: x n_columns: 49"
 ": y_min: -3 y_max: 3 y_inc: 0.125 name: y n_rows: 49"
 ": v_min: -6.54166173935 v_max: 8.08391857147 name: z"
 ": scale_factor: 1 add_offset: 0"
 ": format: classic"
 ": Default CPT: "

# Reinterpolate to 0.1
julia> G2 = gdalwarp(G, I=0.1);

julia> grdinfo(G2)
11-element Vector{String}:
 ": Title: "
 ": Command: "
 ": Remark: "
 ": Gridline node registration used [Cartesian grid]"
 ": Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7"
 ": x_min: -3.0125 x_max: 2.9875 x_inc: 0.1 name: x n_columns: 61"
 ": y_min: -2.9875 y_max: 3.0125 y_inc: 0.1 name: y n_rows: 61"
 ": v_min: -6.54166173935 v_max: 8.08391857147 name: z"
 ": scale_factor: 1 add_offset: 0"
 ": format: classic"
 ": Default CPT: "
9

Thank both of you for the replies!

I ended up using Rasters.warp() because I have to limit the output with a set of bounds defined in the destination CRS, and it seems not possible with Rasters.resample()

10

If you have time please make a github issue spelling out what should happen and why it doesn’t and I’ll fix it. resample should be able to handle this.

1 Like