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 gdal_warp function
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.
EDIT:
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[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
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
gdalinfo INPUT.tif
.