WMS - Web Map Services displays

Specific Domains Geo
1

May I ask if anyone is familiar with WMS displays of data Web Map Service - MATLAB & Simulink - MathWorks United Kingdom
Web Map Service - Wikipedia

I ask for a good reason, as I would like to do some example plots using this in Julia. Yes, I should Google before I ask. However there is sually a great deal of knowledge on this forum!

Thanks in advance

2

As far as I know there are no WMS client libraries in Julia. With PyCall you could try using this:
https://geopython.github.io/OWSLib/

If the data is avaiable as web tiles (slippy map) or vector tiles you can have a look at
https://github.com/yeesian/MapTiles.jl

3

@visr Thankyou. Good information.

4

I am wanting to display wind speeds and wave heights on plots of the ocean, on web based displays. If anyone is doing work like this please drop me a message and give me a street in the correct direction. Nautical pun slightly intended.

5

Supporting also WMTS (Web Map Tile Service) and offline maps using MBTiles format will be a great feature to have in Julia.

MBTiles files can be generated using Java software MOBAC (MOBile Atlas Creator by fetching various map sources and storing maps offline.

MBTiles files can be displayed using Qt/C++ software GPXSee which can also handle WMS (Web Map Service) and WMTS (Web Map Tile Service).

6

This thread is now the top google result for “WMS julia” :grinning:
Here is some small script that might be helpful:

using Dates
using PyPlot
using Statistics
using Images

module WMS
using HTTP
using Images

function map(baseurl,service = "WMS", version = v"1.3.0", request = "GetMap"; kwargs...)
    args = Dict(kwargs)
    args[:service] = service
    args[:request] = request
    args[:version] = string(version)
    args[:time] = string(args[:time],'Z')
    args[:layers] = join(args[:layers],',')
    args[:bbox] = join(string.(args[:bbox]),',')

    r = HTTP.get(baseurl,query = args)
    return Images.load(IOBuffer(r.body))
end
end

west, east = 9, 22
south, north = 32, 43

img = WMS.map("https://gibs-c.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi",
            time = DateTime(2021,10,29),
            layers = ["MODIS_Terra_CorrectedReflectance_TrueColor"],
            format = "image/jpeg",
            bbox = [south,west,north,east],
            width = 512,
            height = 512,
            crs = "epsg:4326")

# Matrix{RGB{N0f8}} -> Array{Float32, 3} (height x width x 3) for imshow
data = float.(permutedims(Images.channelview(img),(2,3,1)))
imshow(data, extent = [west, east, south, north])
gca().set_aspect(1/cosd(mean((south, north))))

This example produces this image with PyPlot. (WMS.map returns a image which should be useful with other plotting frameworks).

Figure_1

3 Likes
7

@Alexander-Barth, this looks great. Would you know what is the resolution of the images?

8

Apparently the maximal spatial resolution is 250 m:

https://rdrr.io/cran/ursa/man/get_earthdata.html

9

Thank you for your time and response. Fyi, searching in the NASA site blog, found this post that indicates that last year they had released 30 m resolution to the public?

The Google Earth images seem to have higher resolution than that and might be accessed throught the EarthEngine.jl package (never tried it).

1 Like
10

Funny, also found that thread and implemented a (semi) WMS in GMT based on GDAL. I’m saying semi because apparently GDAL does not read all WMS services (and I also have quite limited knowledge on the matter).

using GMT
wms = wmsinfo("http://tiles.maps.eox.at/wms?");
img = wmsread(wms, layer=3, region=(-10,-5,37,44), pixelsize=1500); # one can go up to pixelsize=10 meters)
imshow(img)

GMTjl_tmp
GMTjl_tmp1766Ă—2410 380 KB

2 Likes
11

@Alexander-Barth thanks for the pointer to the time. I’ve include that option in GMT wms functions and I can now reproduce your example with

using GMT
wms = wmsinfo("https://gibs-c.earthdata.nasa.gov/wms/epsg4326/best/wms.cgi");
img = wmsread(wms, layer="MODIS_Terra_CorrectedReflectance_TrueColor", region=(9,22,32,43), time="2021-10-29T00:00:00", pixelsize=750);
imshow(img, proj=:guess, dpi=75)  # Degraded resolution to keep it *not-too-big*

And, in this particular regard, I think this (translated/adapted from a Python one) tutorial may be interesting too.
https://www.generic-mapping-tools.org/RemoteS.jl/dev/gallery/HLS/cloud-native-hls-data/

GMTjl_tmp
GMTjl_tmp1832Ă—1762 319 KB

2 Likes
12

Thank you for sharing the tutorial ! I have also been working on a STAC client here: essentially it helps with the JSON parsing and with searching. The idea is that catalog behaves like a julia “dictionary” of subcatalogs.

using STAC
catalog = STAC.Catalog("https://cmr.earthdata.nasa.gov/stac/")
# LPCLOUD is here the id of the sub-catalog
@show catalog["LPCLOUD"]
1 Like
13

Thanks, I’ll look at it in a future revision of that tutorial. I find it to cryptic but that’s how it was in Python and mostly wanted to be able to reproduce it in Julia.

1 Like