sits
sits copied to clipboard
[Probably BUG?] overlap warning when merge two adjacent classified maps
Sits work with tiles and generate products separated For example if i apply sits classify in a cube with two or more tiles, the number of classificated images in tiles is the same of input tiles.. I tried to merge outputs in one images using several packages as merge and mosaic, after open output raster in R. When it executed, a message warning appears, implying that the images have overlapping pixels or have an aligment problem
Steps to Reproduce
- use sits classify to generate a classified imagem for two or more adjacent tiles
- save this tiles somewhere in tif format
- open these tiles as rast or raster object in R
- use merge or mosaic function (terra or raster package) to join adjacent tiles
I think it may be a problem in my processing, but I think it's important to observe if this occurs in any object generated by sits
Hi @ils15!
This problem occurs because you are trying to mosaic non-adjacent tiles. See the figure below:
To mosaic these tiles, use the gdalwarp
function.
Here is an example:
library(gdalUtilities)
files_classified <- list.files("/path/to/dir", pattern = "_class_v1.tif$", full.names = TRUE)
out_file <- "/path/to/dir/my_mosaic.tif"
gdalUtilities::gdalwarp(
srcfile = files_classified,
dstfile = out_file,
co = c("COMPRESS=LZW", "PREDICTOR=2", "BIGTIFF=YES"),
ot = "Byte",
t_srs = "+proj=aea +lat_0=-12 +lon_0=-54 +lat_1=-2 +lat_2=-22 +x_0=5000000 +y_0=10000000 +ellps=GRS80 +units=m +no_defs "
)
Sorry, wrong print
i made a lot of tests and send wrong mosaic
lets check it again
same error with adjacent mosaics
its four adjacent areas
what ive tried
merge and/or mosaic four areas
merge and/or mosaic in pairs
i dont know what is happenend, but i saw BDC uses a albers equal area projection with some edited parameters
i think something is wrong with terra/raster packages interpretation os this srs/crs
But as i said, it can be my processing metodology mistake here
Look at these tiles please
Only way to merge it is using gdal?
Yes, try with gdal and see if it works. The warp admits different alignments.
@OldLipe assuming i use GDALwarp, mosaic will make a resample in cells to make image align right? Because with albers equal area i cannot merge two or more tiles without resample? correct? For example, I assume the risk of changing a pixel when I join two or more tiles (ajacent or not). What I understand is, in the case of using albers equal area and assuming the count is in meters and not in grid (lat long), you made adjacent tiles non paralell for sure. For this reason a resample is required every time we need to merge two or more adjacent areas. To make extent regular again, correct? I confess that I haven't read all the sits documentation (maybe it's detailed there), but it's something that can cause some small problems when treating the image by block. For a classificated image, impact are very small. But for a float image, resample can produce effects where images are joined (edges). Perhaps this explains the care with which the sits and BDC procedures are always performed tile by tile and not in blocks, sure?
Hi @ils15, In this case the images will not be resampled. There is a small difference of the bounding box between the images that is generated by floating point error. The problem is solved by numerical approximation of the spatial resolution, something that the gdal algorithm does very well. This has nothing to do with AEA projection. Regarding the classification procedure, I can assure you that the block strategy is solely for computational reasons. Consider the following scenario: 5000x5000 pixel images, 4 bands and 24 observations in time. The observations are stored as 8byte values in R. Just to keep this data in memory you would need 17GB. If you take into account data handling and classification, you can multiply this by at least a factor of 3. It is not feasible to classify a cube all at once on most computers, and in this case the block processing is the most feasible strategy.
@ils15, It would make it much easier for us if you could provide your sessionInfo()
output. Thus, I can reproduce the warning message and see if it is something related to the version of the terra package.
@rolfsimoes
great if no resample, ok
i make a script with my girlfriend and we have some difficults with handle sits data with another package or export results as well for example i want do make a map with all tiles i need to merge or mosaic it for example
about processing
we used a i7 8gen and a ryzen 7 5800H with ddr4 and 32GB Ram. Processing is hard to do and your strategy is fantastic. I'm only little afraid about sits products outputs and compatibilization with another softwares.
Some things We got looking at the source code here and implementing it by ourserlves (such as handling the "sits" objects and using it with other R libraries), but I confess that the issue of block joins made me a little worried, and it caused some problems when i used functions such as extract and others directly from the generated products. However, since you assured that the "terra" warning is the result of a problem in the approach, no problem anymore. Otherwise, I think it's good to take a look at terra warnings (for example, raster didnt mosaic your outputs) and even recommending the use of gdal in the tile merge process. About AEA projection, the main problem is about compatibilize all input data with it. Cause we need to extract a crs in a sits object and reproject all input data to force coincidence. I know it's not your focus at the moment, but it's a real feedback from a user who has some familiarity with these raster data extraction packages and the difficulties in being able to insert metodology and scripts in sits reality.