1

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

Subset of DEM plotted in Julia

This is the result of the gdal_warp function

Result of the cropping

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
1
  • Can you edit your question and include the output of command-line gdalinfo INPUT.tif.
    – user2856
    Apr 4 at 22:38

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.