I'm trying to save a
.kmz file from a dataset representing a DEM of an area stored in a GeoTiff file.
My code is reading the DEM, transforming the CRS to UTM, and defining a subset based on coordinates. On this subset I did some computation based on the DEM altitude.
I'd like to save the result in a
kmz file to be opened in Google Earth using Julia.
What I've tried up to now is to use ArchGDAL together with Rasters.jl. However the results is not what I wanted, I'm still missing the transparency for the missing values. I've found a similar question, but it didn't explain how he fixed.
The gdal command that should do the job is
gdalwarp -srcnodata <in> -dstnodata <out> -crop_to_cutline -cutline INPUT.shp INPUT.tif OUTPUT.tif
I've tried to replicate it with Rasters.jl, having defined the the cropline as a GeoJSON
A = Raster(dem_file) flags = Dict( "-srcnodata"=> NaN, "-dstnodata" =>0, "-crop_to_cutline" => "", "-cutline" =>"./shape.json", ) B = warp(A, flags) write("out.kmz",B)
The subset I'm working on in Julia is a structure with the same dimensions of the original DEM, but to define the subset I have set all the values out of the bounds to
NaN. This is the subset plotted in Julia with the red crosses represetning the 4 points in Polygon of the cutline
This is the result of the
What I can see is that
gdal_warp is cropping correctly the original DEM file but it is not considering the
nodata field, and it is mapping it to a the black color. I've tried changing the values but it seems to not affect the output.
I include the output of
gdalinfo INPUT.tif as requested by @user2856
Driver: GTiff/GeoTIFF Files: copernicus_california_wgs84_small.tif Size is 3601, 3601 Coordinate System is: GEOGCRS["WGS 84", DATUM["World Geodetic System 1984", ELLIPSOID["WGS 84",6378137,298.257223563, LENGTHUNIT["metre",1]]], PRIMEM["Greenwich",0, ANGLEUNIT["degree",0.0174532925199433]], CS[ellipsoidal,2], AXIS["geodetic latitude (Lat)",north, ORDER, ANGLEUNIT["degree",0.0174532925199433]], AXIS["geodetic longitude (Lon)",east, ORDER, ANGLEUNIT["degree",0.0174532925199433]], ID["EPSG",4326]] Data axis to CRS axis mapping: 2,1 Origin = (-118.743838888888888,34.552338888888890) Pixel Size = (0.000277777777778,-0.000277777777778) Metadata: AREA_OR_POINT=Area Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left (-118.7438389, 34.5523389) (118d44'37.82"W, 34d33' 8.42"N) Lower Left (-118.7438389, 33.5520611) (118d44'37.82"W, 33d33' 7.42"N) Upper Right (-117.7435611, 34.5523389) (117d44'36.82"W, 34d33' 8.42"N) Lower Right (-117.7435611, 33.5520611) (117d44'36.82"W, 33d33' 7.42"N) Center (-118.2437000, 34.0522000) (118d14'37.32"W, 34d 3' 7.92"N) Band 1 Block=3601x1 Type=Float32, ColorInterp=Gray