Add resample wrapper and methods to convert between ArchGDAL.Dataset and GeoArray #73
Conversation
GeoArray to ArchGDAL.dataset and vice versa
| function unsafe_ArchGDALdataset(A::GeoArray) | ||
| width = size(A, Lon) | ||
| height = size(A, Lat) | ||
| nbands = size(A, Band) |
There was a problem hiding this comment.
We can't actually be sure that a GeoArray has any of these dimensions - especially Band but the others as well as it could have been sliced. You can use e.g. if hasdim(A, Lon) to check if A has a Lat dim.
You can make a Dataset from 2d without a Band - check the write methods.
The GeoArray may also have other dimensions (that you can get with otherdims(A, (Lat, Lon, Band))) I was wondering if it would be cool to just loop over the time/vertical etc dimension somewhere and slice it - then convert each slice with GDAL and cat them together afterwards. But don't worry about that now :)
There was a problem hiding this comment.
Ahh I see -- is there a way to identify the X dimension, the Y dimension, and the number of bands without relying on dim names?
There was a problem hiding this comment.
They should be Lat/Lon/Band for GeoArray - but a GeoArray might be just a vector after indexing or something and not have them all anymore. And other formats like netcdf wont have Band - they might just be 2d lat/lon arrays.
But maybe don't worry about this for now I'll fix it up later - I'm distracting you from the main use case with the edge cases.
also update GeoArray() function to use a Symbol for the name arg
|
|
||
| ## Keyword Arguments | ||
| - `crs`: A `GeoFormatTypes.GeoFormat` specifying an output crs (`A` with be reprojected to `crs` in addition to being resampled). Defaults to `crs(A)` | ||
| - `method`: A `String` specifying the method to use for resampling. Defaults to `"near"` (nearest neighbor resampling). See [resampling method](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r) in the gdalwarp docs for a complete list of possible values. |
There was a problem hiding this comment.
@rafaqz this link does not render when I run ?resample -- but maybe it will for the Documenter.jl-generated docs? Just wanted to flag this.
Codecov Report
@@ Coverage Diff @@
## master #73 +/- ##
==========================================
+ Coverage 76.70% 79.49% +2.78%
==========================================
Files 16 17 +1
Lines 790 878 +88
==========================================
+ Hits 606 698 +92
+ Misses 184 180 -4
Continue to review full report at Codecov.
|
|
Ok this is looking great. Just needs a few tests! |
|
Also I was thinking I will add I have an immediate use case for your snap raster idea as well - so I might have a go at adding that too. Were you thinking we could just use the "-te" and "-te-srs" flags in gdalwarp to add new extent bounds? |
Yep, exactly! Though it might tricky if the snap raster bounds and input raster bounds don't match. Those edge cases were what prevented me from going for it right away. Though I suppose each could be padded with NoData so the extents do match? That feel a little hacky though. For tests, I can think of a couple. 1) check that the resolution of the output matches the resolution specified in the function call and 2) check that the extent bounds are the same for the input and output GeoArrays (allowing error up to +/- the step of each dim). |
|
Yeah I was thinking about that too. Its easy if the input raster is larger in every bound. Padding the edge is an option. I'm not sure what GDAL will do if you specify "te" outside of the bounds of the dataset you pass it? Does it just fill with NoData anyway? (I've never actually used it manually, lol) And those tests sound about right. Maybe we could also test something using the command line gdalwarp on a tiff to make sure they match? Not sure if that's overkill, but it would be conclusive. |
I imagine ArchGDAL tests for that, so maybe testing |
|
Oh right that makes more sense. |
Okay I'll add that test against ArchGDAL too.
I'm not 100% sure either. I remember this being an issue in GeoData when subsetting (I know that's a totally different animal though) for another workflow. I had to do something like: left_bound = max(left_extent_raster1, left_extent_raster2)
right_bound = min(right_extent_raster1, right_extent_raster2)
lower_bound = max(lower_extent_raster1, lower_extent_raster2)
upper_extent = min(upper_extent_raster1, upper_extent_raster2)Here's the complete code if you want to check it out. Something similar might be needed for setting |
|
Also I did use gdalwarp with a snap raster in that same workflow, the extent of the snap (and also mask) raster was just always within the bounds of the input raster so I didn't have to worry about cases of non-overlapping rasters. Here's that code. |
|
@rafaqz what are you using for data sources for testing? I'll need a .tif to test ArchGDAL output -- should I create a GeoArray, write it to test/data, then read it back in with ArchGDAL? |
|
I just use some random tif file from somewhere... It's in the gdal.jl test file. I need to source more and do that in an organised way. |
|
@rafaqz tests added! I skipped checking the output extents for now (which is a test I suggested above) -- I just test for the correct output resolution and that the data (array) for the |
|
Awesome! I'll add some more when I have a go at snap. It seems the download is broken somehow, but I can merge this when that passes. |
|
Sweet! Yeah it's something simple with |
|
🎉 |
|
Looks like gdal just fills the edge :) I'll get a PR sorted tomorrow. |
|
Nice! That makes things super easy! |
Closes #68. Still a WIP -- no formal tests or docs yet, but can confirm that it is working for resampling a single-band raster from 1km resolution in a UTM projection (left panel in figure below) to 0.002° in EPSG:4326 (right panel in figure below) 🎉. I think it should work for multi-band rasters. An edge case that would have errors would be if different bands have different NoData values, but I'm not sure if GeoData handles that to start.
I had some issues/trouble converting
ArchGDAL.Datasets andRasterDatasets toGDALarrays, so I just skipped the middle man since theArchGDAL.Datasetis being created in memory anyway. I added a method to convert back and forth betweenArchGDAL.DatasetandGeoData.GeoArray.@rafaqz I wanted to give you a chance to review what I have so far before proceeding further.