Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Closed
Rapsodia86 opened this issue Nov 2, 2022 · 6 comments
Labels

Comments

@Rapsodia86
Copy link

Rapsodia86 commented Nov 2, 2022

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 yeesian changed the title Convert GDAL.OGREnvelope object to array or dataframe Should ArchGDAL.envelope(geom) return Vector{Float64} instead of GDAL.OGREnvelope? Nov 3, 2022
@yeesian
Copy link
Owner

yeesian commented Nov 3, 2022

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

@evetion
Copy link
Collaborator

evetion commented Nov 3, 2022

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
Copy link
Author

Rapsodia86 commented Nov 3, 2022

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
Copy link
Author

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
Copy link
Contributor

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
Copy link
Author

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!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

4 participants