How to drop the Z-dim of a shape file

Specific Domains Geo
1

I would like to rasterise and draw border of a polygon from a Shapefile as in this example from Rasters.jl.

I load my Shapefile as following

julia> vd_bord = Shapefile.Handle("data/Shapefiles/Vaud.shp").shapes[1]
Shapefile.PolygonZ(Shapefile.Rect(494319.65625, 115136.65125000104, 585425.3537500016, 198589.62999999896), Int32[0, 1615, 1633, 1663], Shapefile.Point[Shapefile.Point(519503.5712499991, 176197.8575000018), Shapefile.Point(519799.44249999896, 176509.0162499994), Shapefile.Point(519919.7650000006, 176744.0450000018), Shapefile.Point(519980.78500000015, 176863.1887499988), Shapefile.Point(519971.28125, 177153.04374999925), Shapefile.Point(519748.66750000045, 177428.7412500009), Shapefile.Point(520216.8062499985, 177940.73999999836), Shapefile.Point(520693.0062500015, 178248.2674999982), Shapefile.Point(520965.76000000164, 178169.1712500006), Shapefile.Point(521062.1087500006, 178194.1950000003)  …  Shapefile.Point(551995.2474999987, 176404.64375000075), Shapefile.Point(551935.2474999987, 176374.6424999982), Shapefile.Point(551815.245000001, 176284.6387500018), Shapefile.Point(551755.245000001, 176264.63749999925), Shapefile.Point(551656.2412500009, 176217.63500000164), Shapefile.Point(551495.2462499999, 176400.63500000164), Shapefile.Point(551038.2274999991, 176067.625), Shapefile.Point(550923.2375000007, 176332.625), Shapefile.Point(550798.2287500016, 176179.62124999985), Shapefile.Point(550702.2312499993, 176260.62124999985)], Shapefile.Interval(-1.4551915228366852e-11, -1.4551915228366852e-11), [-1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11  …  -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11, -1.4551915228366852e-11], Shapefile.Interval(-1.7976931348623157e308, -1.7976931348623157e308), [-1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308  …  -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308, -1.7976931348623157e308])

Then rastserize it as following (dim I get from a geotif file that has correct extent)

vd = rasterize(vd_border; to=dim, missingval=-9999, fill=1, order=(X, Y), shape=:polygon, boundary=:touches)

Plotting vd looks good:

p = plot(vd)

Screenshot 2022-12-07 at 21.52.33
Screenshot 2022-12-07 at 21.52.33830×682 64.6 KB

But when I try adding the border, then it shows a buggy behaviour:

plot!(p, vd_border; fillalpha=0, linewidth=0.6)

Screenshot 2022-12-07 at 21.52.59
Screenshot 2022-12-07 at 21.52.59830×682 88.1 KB

I guess it’s because my polygon vd_border has a Z-dim

julia> GeoInterface.extent(vd_border)
Extent(X = (494319.65625, 585425.3537500016), Y = (115136.65125000104, 198589.62999999896), Z = (-1.4551915228366852e-11, -1.4551915228366852e-11))

Any hint how I could drop the Z-dim of the polygon? Thanks

2

I’m not sure there is an easy way to drop the third dimension (but there could be something like GeoInterface.dropz or similar).

If you post your shapefile here we can work out one of the more difficult ways of doing it.