Should ArchGDAL.envelope(geom) return Vector{Float64} instead of GDAL.OGREnvelope?

[Rapsodia86]

UPDATE! Ok, I figured it out. I used the original GDAL function as I was working with raster, and created an array as a pointer:

shp_ds = GDAL.gdalopenex(shp_dir, GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)
layer1 = GDAL.gdaldatasetgetlayer(shp_ds , 0)
transformE = Array{Cdouble}(undef, 4)
GDAL.ogr_l_getextent(layer1,pointer(transformE),false)

Hi! I am a bit struggling to get shapefile extent coordinates and projection. If had a raster, I would do:

dt= GDAL.gdalopen(raster_in, GDAL.GA_ReadOnly)
transform = Array{Cdouble}(undef, 6)
result = GDAL.gdalgetgeotransform(dt, pointer(transform))
xmin = transform[1]
ymax = transform[4]
xmax = xmin + transform[2] * GDAL.gdalgetrasterxsize(dt)
ymin = ymax + transform[6] * GDAL.gdalgetrasterysize(dt)
pr_m = GDAL.gdalgetprojectionref(dt)

For shapefile, I do:

shp_in=ArchGDAL.read(shp_dir)
layer = ArchGDAL.getlayer(shp_in, 0)
#to get projection in WKT format 
spatial_ref = ArchGDAL.getfeature(layer, 0) do feature
    ArchGDAL.getgeom(feature, 0) do geom
        source = ArchGDAL.getspatialref(geom) do source
        spref = ArchGDAL.toWKT(source)
        end
    end
end

Without do-blocks approach, I am getting error:

geom = ArchGDAL.getfeature(layer, 0)
ERROR: MethodError: no method matching getfeature(::ArchGDAL.IFeatureLayer, ::Int64)

To get the extent of the shapefile (it is not rectangle but administration unit, but I need extent)

ext = ArchGDAL.getfeature(layer, 0) do feature
    ArchGDAL.getgeom(feature, 0) do geom
        envelope_ext = ArchGDAL.envelope(geom)
    end
end

And here is the main problem: I cannot convert GDAL.OGREnvelope object to array nor dataframe. What am I missing? Sorry if that is basic!

[yeesian]

Thank you for filing this issue, and proposing an approach (in https://github.com/yeesian/ArchGDAL.jl/issues/343#issue-1433468095) that I haven't thought of:

shp_ds = GDAL.gdalopenex(shp_dir, GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)
layer1 = GDAL.gdaldatasetgetlayer(shp_ds , 0)
transformE = Array{Cdouble}(undef, 4)
GDAL.ogr_l_getextent(layer1,pointer(transformE),false)

My main concern with returning a vector rather than GDAL.OGREnvelope is that users will need help interpreting it as [MaxX, MaxY, MinX, MinY], rather than being able to inspect it in the terminal with tab-completion like so:

julia> envelope = GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)
GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)

julia> envelope.M
MaxX  MaxY  MinX  MinY

julia> [envelope.MaxX, envelope.MaxY, envelope.MinX, envelope.MinY]
4-element Vector{Float64}:
 100.0
  70.0
 100.0
  70.0

[evetion]

May I kindly ask why you are working directly with GDAL?

I ask because we've abstracted away many of these C like calls in ArchGDAL, and further abstractions exist in GeoInterface, Rasters, GeoArrays and GeoDataFrames. (disclosure, I'm the author of the latter two). Are we missing functionality in these higher level packages?

For example, there's GeoInterface.extent, which gives you a more generic extent (namedtuple like) from the Extents package, that should work on Archgdal geometries.

[Rapsodia86]

Thank you for filing this issue, and proposing an approach (in #343 (comment)) that I haven't thought of:

shp_ds = GDAL.gdalopenex(shp_dir, GDAL.GDAL_OF_VECTOR, C_NULL, C_NULL, C_NULL)
layer1 = GDAL.gdaldatasetgetlayer(shp_ds , 0)
transformE = Array{Cdouble}(undef, 4)
GDAL.ogr_l_getextent(layer1,pointer(transformE),false)

My main concern with returning a vector rather than GDAL.OGREnvelope is that users will need help interpreting it as [MaxX, MaxY, MinX, MinY], rather than being able to inspect it in the terminal with tab-completion like so:

julia> envelope = GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)
GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)

julia> envelope.M
MaxX  MaxY  MinX  MinY

julia> [envelope.MaxX, envelope.MaxY, envelope.MinX, envelope.MinY]
4-element Vector{Float64}:
 100.0
  70.0
 100.0
  70.0

Oh, I wish I knew about the envelope.MaxX, then I would not ask and bother you. I should have checked the attributes on GDAL OGREnvelope class. I am sorry for that! I never before used that function. That really solves my issue without making any changes and makes totally sense.

However, I am getting an error when try envelope.M

julia> envelope = GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)
GDAL.OGREnvelope(100.0, 100.0, 70.0, 70.0)

julia> envelope.M
ERROR: type OGREnvelope has no field M
Stacktrace:
 [1] getproperty(x::GDAL.OGREnvelope, f::Symbol)
   @ Base .\Base.jl:42
 [2] top-level scope
   @ REPL[45]:1

julia> [envelope.MaxX, envelope.MaxY, envelope.MinX, envelope.MinY]
4-element Vector{Float64}:
 100.0
  70.0
 100.0
  70.0

Forgot to mention, my ArchGDAL version is v0.8.4 and GDAL v1.3.0. Looks like it is time to make an update! May I close the issue?

[Rapsodia86]

May I kindly ask why you are working directly with GDAL?

I ask because we've abstracted away many of these C like calls in ArchGDAL, and further abstractions exist in GeoInterface, Rasters, GeoArrays and GeoDataFrames. (disclosure, I'm the author of the latter two). Are we missing functionality in these higher level packages?

For example, there's GeoInterface.extent, which gives you a more generic extent (namedtuple like) from the Extents package, that should work on Archgdal geometries.

I started learning Julia and working with, by rewriting my python scripts where GDAL was the library I was using. Therefore, I have been using ArchGDAL.jl and GDAL.jl since that time. I never had a chance or need to explore other libraries. Most of my work is in R or python, unless I need some heave duty jobs to be done. Then Julia steps in. So, I am not familiar with those other packages you mentioned. When have an opportunity, I will surely take a look! Thank you for pointing them out!

[felixcremer]

However, I am getting an error when try envelope.M You need to press tab after you wrote envelope.M to see all fields of the envelope type, which start with M

[Rapsodia86]

However, I am getting an error when try envelope.M You need to press tab after you wrote envelope.M to see all fields of the envelope type, which start with M

Yes, right! My mistake! Thanks!