GMT Image to ImageView

Specific Domains Geo
,
1

Thankyou to @joa-quim I can read Geotiff files using GMT.jl
When I run image = gmtread(“world.tif”,img=true)
the image is of type GMT.GMTimage
GMT.imshow(image) produces a Postscript file

Rather stupidly I would like to use ImageView to display the image quickly.
https://github.com/JuliaImages/ImageView.jl
ImageView.imshow() does not have a method for the GMTImage type

Can we work on a conversion to the array types which Image and ImageView like to work with?
I am probably missing the fact that the GMTImage format has lots of metadata about the coordinates etc.

2

Yes, being able to display with ImageView would be very nice. I even started to look at it once on how to add to Images (condition to display with ImageView) but I think I didn’t really understand what I had to do.

If you look at this line of the source code it even anticipates the memory layout that would be needed to integrate with Images. If someone with more knowledge of Images wanted to give an help, it would be fantastic.

3

Does gmt only output ps?

4

It produces lots of things but in terms of mapping, basically yes. One of it’s modules (grdimage) is able to write raster files (many of the formats supported by GDAL) but without any vector content. However, it has a module that converts (via ghostscript) the ps into may other formats.

5

I know this example is not using GMT but ArchGDAL to read the GeoTIFF, but perhaps a similar method can be used for GMT:

# add ArchGDAL#idataset
using ArchGDAL
using ImageCore
using ImageView
using ColorTypes

# read the file into a dataset
filename = joinpath(dirname(pathof(ArchGDAL)), "../test/pyrasterio/RGB.byte.tif")
dataset = ArchGDAL.read(filename)

# read all of the 3 bands as UInt8 2D arrays and view them together as colors
# see also https://juliaimages.org/latest/conversions_views.html
img = permuteddimsview(colorview(RGB,
    normedview(ArchGDAL.read(dataset, 1)),
    normedview(ArchGDAL.read(dataset, 2)),
    normedview(ArchGDAL.read(dataset, 3))
), (2,1))

imshow(img)
6

It looks like it should do it, thanks. (have to go now for some hours).

7

I’m having this error

julia> filename = joinpath(dirname(pathof(ArchGDAL)), "../test/pyrasterio/RGB.byte.tif")
"C:\\Users\\j\\.julia\\packages\\ArchGDAL\\tUwz2\\src\\../test/pyrasterio/RGB.byte.tif"

julia> dataset = ArchGDAL.read(filename)
ERROR: MethodError: no method matching read(::String)
You may have intended to import Base.read
Closest candidates are:
  read(::ArchGDAL.Dataset) at C:\Users\j\.julia\packages\ArchGDAL\tUwz2\src\raster\rasterio.jl:308
  read(::ArchGDAL.Dataset, ::Array{Int32,1}) at C:\Users\j\.julia\packages\ArchGDAL\tUwz2\src\raster\rasterio.jl:301
  read(::ArchGDAL.RasterBand) at C:\Users\j\.julia\packages\ArchGDAL\tUwz2\src\raster\rasterio.jl:195
8

You’ll need to be on the idataset branch in Provide Interactive Datasets for REPL usage by yeesian · Pull Request #76 · yeesian/ArchGDAL.jl · GitHub

9

Ah, guessed something like that. Thanks.

10

@johnh, this is undocumented territory and things may change (for example, I don’t like the idea of making a image copy), but for now you can do this

julia> using GMT, ImageCore, ColorTypes, ImageView
julia> I = gmtread("the_file_name.tif", img=true, layout="TCP");
julia> img = colorview(RGB, normedview(I.image));
julia> ImageView.imshow(img);
11

I don’t like the idea of making a image copy

The view in both normedview and colorview is a hint that both are non-copying.

12

That’s right, but on the other hand the normed part makes think on a true normalization, which would imply the creation of a new array … but it doesn’t although it shows up at such in the REPL. Find this confusing to grok.

julia> c = zeros(UInt8,3,1,1);

julia> c[1] = 255;

julia> cv = colorview(RGB, normedview(c));

julia> cv[1]
RGB{N0f8}(1.0,0.0,0.0)

julia> c[1] = 127;

julia> cv[1]
RGB{N0f8}(0.498,0.0,0.0)

So when is the normalization taking place? Only for printing in REPL?

13

It refers to normed numbers, GitHub - JuliaMath/FixedPointNumbers.jl: fixed point types for julia. In most image-processing frameworks the value of white (for a grayscale pixel) is representation dependent:

white == 255 #= 8-bit camera =# == 65535 #= 16-bit camera =# == 4095 #= 12-bit camera =# == 1.0 #= floating-point representation =#

While many people who work with images have gotten used to thinking this way, it’s fair to say it’s rather bad mathematics: I’m not aware of any other domain of computation where the above equalities hold (255 == 1.0, seriously?). It also violates our general expectation that “meaning” is independent of “representation”: just do convert(Float64, x) and suddenly you’re very confused. Particularly for cameras with 10-, 12-, and 14-bit depths, it’s very hard to automatically detect pixel-saturation because the maximum output is less than typemax(UInt16), unless you always carry around an extra bit of metadata about the camera.

We’ve gone to the trouble to fix this by introducing a set of numbers, the normed numbers, that allow you to represent 1.0 with an 8-bit number, or a 16-bit number, or the first 10 bits of a 16-bit number, etc. These numbers are just reinterpretations of UInt8, UInt16, etc, but they provide consistent rules of arithmetic in which they behave like their floating-point cousins. Consequently, they fix image mathematics, and white always means 1.0 no matter with what precision you choose to represent it.

5 Likes
14

Thanks for the detailed explanation. I guess that the harder part for me to understand is how this happens without creating a new array.

15

ReinterpretArray and MappedArrays are magical. In Julia we can:

  • change the values or apparent element type of an array
  • choose a subset of an array
  • reshape an array
  • permute the dimensions of an array
  • change the “location” of an array by adjusting the indices to access its elements

all without taking a copy. This makes Julia quite amazing for big data, where raw arrays are too big to hold in memory yet you can do some remarkable transformations efficiently without ever creating a new disk file or running out of memory.

2 Likes
16

And the magician works for free?
Never mind this is all good for me and it means that images read with GMT can be displayed/manipulated directly without further work.

17

And the magician works for free?

Mostly, yes. For example:

julia> using MappedArrays

julia> A = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> B = mappedarray(x->cbrt(x), y->y^3, A)
2×2 mappedarray(x->(Main.cbrt)(x), getfield(Main, Symbol("##26#28"))(), ::Array{Int64,2}) with eltype Float64:
 1.0      1.25992
 1.44225  1.5874 

julia> getfirst(C) = @inbounds C[1,1]
getfirst (generic function with 1 method)

julia> @code_llvm debuginfo=:none getfirst(B)

define double @julia_getfirst_12660(%jl_value_t addrspace(10)* nonnull align 8 dereferenceable(8)) {
top:
  %1 = addrspacecast %jl_value_t addrspace(10)* %0 to %jl_value_t addrspace(11)*
  %2 = bitcast %jl_value_t addrspace(11)* %1 to %jl_value_t addrspace(10)* addrspace(11)*
  %3 = load %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)* addrspace(11)* %2, align 8
  %4 = addrspacecast %jl_value_t addrspace(10)* %3 to %jl_value_t addrspace(11)*
  %5 = bitcast %jl_value_t addrspace(11)* %4 to i64 addrspace(13)* addrspace(11)*
  %6 = load i64 addrspace(13)*, i64 addrspace(13)* addrspace(11)* %5, align 8
  %7 = load i64, i64 addrspace(13)* %6, align 8
  %8 = sitofp i64 %7 to double
  %9 = call double @julia_cbrt_12551(double %8)
  ret double %9

This is exactly the same amount of work that a fully-optimized “eager” implementation of the cube-root would require. There is literally zero overhead (other than being somewhat harder on the poor compiler). Now, there can be advantages and disadvantages to lazy computation, but Julia’s implementations are magical in the sense that they usually have no cost for the abstraction itself.

1 Like