Setting color mapping/transparency when writing a kmz file

[stefanomoro]

I have a DEM in tif format. I can read it and create a Raster object. I can write it as a kmz file.

dem = Raster("dem_file.tif")
write("output.kmz",dem)

The result is a black-and-white image. I want to set the color mapping, as I can do with a plot function.
Is it possible?

[rafaqz]

How does GDAL do it? You can pass an options dict as a keyword to writeif that helps. It will be sent to GDAL.

And I guess you need three bands for color, when DEM is usually one?

[stefanomoro]

Thank you for the quick reply!
My idea was to produce a kmz file with the same palette of :terrain, as done when plotting. I would also like to represent the NaN values as transparent.
I will look into GDAL to see what it's available.

[rafaqz]

Maybe ArchGDAL already handles this? It already has imread. We could add writing of a Color matrix to GDAL bands.

Did you try just converting the raster to a colorscheme? (first thing you need to do)

Like colorraster = map(x -> ColorSchemes.terrain[x], raster) should do it, once you have normalised. But yeah, write probably wont work. You may need to split it into three bands, which you can also just do with similar broadcasts and cat(r, g, b; dims=Band).

[stefanomoro]

Thank you again for the tips! I'll try it and report what I can get, hoping it can help others.

[stefanomoro]

I'm having some issues with the output kmz transparency.

I have scaled the single Band to (0,255) first. Then, I replicated the same layer on 3 Bands to get a grey-scale output (to simplify the procedure).
Then, I created a 4 Bands Raster based on a guess that the 4th should be mapped to the Alpha channel.

alpha = map(x -> isnan(x) ? 0 : 255 , A)
B = cat(A, A, A, alpha, dims=Band)

but the result doesn't have any transparency. Any idea?

[rafaqz]

I have never done anything like this in my life ;)

[stefanomoro]

I came up with a possible (dumb) solution.
If I know the bounds of the interesting area, I could cut the Raster object and save it to kmz only that area.
I'm trying to use warp to resample the area I wish, but it doesn't for this case. It is assuming the bounds are aligned with the X,Y axis.
Like in this case, the red cross are the vertex of the bounds:

Any idea on how to crop using a generic polygon?

[rafaqz]

crop ?

You cant crop a raster on anything but the x and y axis... if you want to rotate it, you will need to rotate with warp as you will be resampling.

[stefanomoro]

Searching online, I've found a possible solution to the problem.
The GDAL command for cropping a raster to a specific polygon is like this

gdalwarp -srcnodata <in> -dstnodata <out> -crop_to_cutline -cutline INPUT.shp INPUT.tif OUTPUT.tif

where <in> is the value in the dataset associated to NoData and the <out> is the one for the output.
-crop_cut_line specify to cut the dataset on the shapefile, and INPUT.shp is the cutline to be used.

I tried to replicate this with Rasters.warp(), having a dataset A with the NaN for the missing values:

flags = Dict(
        "-srcnodata"=> NaN, "-dstnodata" =>0, 
        "-crop_to_cutline" => "", "-cutline" =>"./shape.json",
        )
B = warp(A, flags)   
write("out.kmz",B)

I used a GEOJSON file for the shape, defined like this

{
 "coordinates": [
  [
   [-118.08863908451265,    34.265122559375975   ],
   [    -118.089138690936,    34.26432188457904   ],
   [    -118.08489975017962,    34.26249852189752   ],
   [    -118.08440011340001,    34.26329918021013   ],
   [    -118.08863908451265,    34.265122559375975   ]
  ]
 ],
 "type": "Polygon"
}

The results is not what I expected. The cut is the same I had when defining a bounding box with only 2 corners.
Moreover, it seems the nodata fields are not working correctly.

Any idea?

[rafaqz]

Crop usually means to a square? I think GDAL will does exactly what crop does here. So maybe what you want is crop and mask? Check the docs for syntax, but can just pass you geojson to them and it should work.

Rasters doesnt use gdal for anything but loading files and warp, which are both hard to implemeny natively.

And warp is probived as a power user tool, its not something I can provide help with or ever use myself.

To be clear: these questions are too detailed and case spevific to ask a single package auther on github: if the functions her dont work for you, a gdal forum is the place for you to get answers

[stefanomoro]

Thank you for the answer. Appreciate it!
I'll try it, but I'll probably have to roll back to ArchGDAL for the kmz manipulation.

[rafaqz]

I had an idea about the missing values... if you can post a file with download in the MWE I could actually look at your problem ;)

(You will always get 1000x better help if there is actually a file and an MWE that downloads it so I can just run it and look myself, otherwise it's wild speculation)

[stefanomoro]

Sure, thank you for the help. I can share the DEM Raster object with the missing values. I saved it as GeoTiff.
I have mapped the values to 0-255 and set the missing values as NaN

Here is a small example, just trying to save the DEM to kmz, with the missing values being mapped not correctly.

A = Raster("https://drive.google.com/uc?export=download&id=1xwbaGLT-mY8v-PC0PkPUKrtGMOEFiSBH")
plot(A)

write("out.kmz",A)