2

I am able to export data from a WMS service using QGIS. I would like to be able to do the same using GDAL command line tools.

Below is a reproducible setup of my current method, which I hope to replace.


QGIS setup

The Swedish Geological Survey provides this WMS: https://resource.sgu.se/service/wms/130/jorddjupsmodell

To visualize it in QGIS, I do: Layer -> Add Layer -> Add WMS/WMTS Layer... -> New, name layer as soildepth, provide above URL, Connect, then select Jorddjup 10x10m raster, skattat jorddjup till berg (m) from list, Add.

To have an example polygon, I create testArea.geojson in Python:

import geopandas as gpd
from shapely.geometry import Polygon

gpd.GeoSeries(Polygon([[15, 60],[15, 60.01], [15.01, 60.01], [15.01, 60]])).to_file("testArea.geojson")

Adding this file to QGIS canvas too: Layer -> Add Layer -> Add Vector Layer..., select testArea.geojson, Add.

Result:

enter image description here


Manual data export

Right-click on the Jorddjup 10x10m raster, skattat jorddjup till berg (m) layer -> Export -> Save As....

Then:

  1. Unclick "Create VRT" checkbox
  2. Give a File name, such as /path/to/file/soilTest.tif
  3. Select CRS: Project CRS: EPSG:3006 - SWEREF99 TM
  4. Extent: click on Calculate from Layer, choose testArea
  5. Resolution: write 10 for both Horizontal and Vertical fields

enter image description here

Then click OK. Result is as expected:

enter image description here


How to export without QGIS?

I would like to perform the above data export using command line tools only. It is important that the final result should have projection EPSG:3006, it's extent should be the extent of testArea.geojson, and it's resolution should be 10 (m) in both directions. (Ie QGIS steps 3, 4 & 5 above.)

How can I perform the above data export using command line tools?


Attempts using gdal_translate

I.

The suggestion of this answer is answer:

gdal_translate -of GTiff "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262" -tr 10 10 /tmp/test.tif

Checking the result with Python:

import rasterio
with rasterio.open("/tmp/test.tif") as src:
    src_read = src.read()

src_read is:

array([[[255, 255]],
       [[252, 255]],
       [[255, 255]]], dtype=uint8)

so the command needs some improvement.

II.

A somewhat naive attempt:

gdal_translate -of GTiff -tr 10 10 -projwin 15.0 60.01 15.01 60.0 -projwin_srs EPSG:4326 "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262" /tmp/test.tif

Here I use -projwin to specify the bounding box of the area I am interested in, and -projwin_srs to specify the projection of this bounding box. The [-tr] is there to specify that I want 10 by 10 m pixels. This is what the documentation says about tr:

set target resolution. The values must be expressed in georeferenced units. Both must be positive values. This is mutually exclusive with -outsize and -a_ullr.

I don't understand how to specify the units for this.

The above command returns:

Input file size is 1073741824, 826238570
ERROR 5: Invalid output size: 0.501 x 0.501

III.

Instead of using geographical coordinates, can try using pixel coordinates using the srcwin option.

gdalinfo "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262"

returns:

Driver: WMS/OGC Web Map Service
Files: none associated
Size is 1073741824, 826238570
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ID["EPSG",6326]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]
Data axis to CRS axis mapping: 1,2
Origin = (7.395557961066897,69.413672039902622)
Pixel Size = (0.000000017633713,-0.000000017633713)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   7.3955580,  69.4136720) (  7d23'44.01"E, 69d24'49.22"N)
Lower Left  (   7.3955580,  54.8440179) (  7d23'44.01"E, 54d50'38.46"N)
Upper Right (  26.3296135,  69.4136720) ( 26d19'46.61"E, 69d24'49.22"N)
Lower Right (  26.3296135,  54.8440179) ( 26d19'46.61"E, 54d50'38.46"N)
Center      (  16.8625857,  62.1288450) ( 16d51'45.31"E, 62d 7'43.84"N)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  Overviews: 536870912x413119285, 268435456x206559643, 134217728x103279821, 67108864x51639911, 33554432x25819955, 16777216x12909978, 8388608x6454989, 4194304x3227494, 2097152x1613747, 1048576x806874, 524288x403437, 262144x201718, 131072x100859, 65536x50430, 32768x25215, 16384x12607, 8192x6304, 4096x3152, 2048x1576, 1024x788
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  Overviews: 536870912x413119285, 268435456x206559643, 134217728x103279821, 67108864x51639911, 33554432x25819955, 16777216x12909978, 8388608x6454989, 4194304x3227494, 2097152x1613747, 1048576x806874, 524288x403437, 262144x201718, 131072x100859, 65536x50430, 32768x25215, 16384x12607, 8192x6304, 4096x3152, 2048x1576, 1024x788
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Blue
  Overviews: 536870912x413119285, 268435456x206559643, 134217728x103279821, 67108864x51639911, 33554432x25819955, 16777216x12909978, 8388608x6454989, 4194304x3227494, 2097152x1613747, 1048576x806874, 524288x403437, 262144x201718, 131072x100859, 65536x50430, 32768x25215, 16384x12607, 8192x6304, 4096x3152, 2048x1576, 1024x788

3rd line says: Size is 1073741824, 826238570. Let's extract a 500 by 500 raster from the middle (1073741824/2=536870912):

gdal_translate -of GTiff -srcwin 536870912 536870912 500 500 "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262" /tmp/test.tif

Check results:

import rasterio
import numpy as np
import matplotlib.pyplot as plt

with rasterio.open("/tmp/test.tif") as src:
    src_read = src.read()

plt.imshow(np.swapaxes(src_read,0,2))

enter image description here

The output is actually a raster now, not an error, which is good.

First problem with this, is that if I change the command to:

gdal_translate -of GTiff -srcwin 100000000 100000000 500 500 "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262" /tmp/test.tif

I get the same result. Not good. Second, I am not able to translate the geographical coordinates of the bounding box to pixel coordinates of the WMS (I know how to do it on a raster image, see this, but I haven't yet done pixel coordinates on WMS). I think I need to understand better the output of gdalinfo in order to do that.

3
  • I think the problem is the server seems to lie about it's coordinates and/or SRS. QGis is unable to import the map either.
    – Ian Turton
    Oct 20, 2022 at 10:59
  • For me QGIS does import it I think. Export also works within QGIS (first and second screenshot).
    – zabop
    Oct 20, 2022 at 11:05
  • Maybe I misunderstand what you mean by "import".
    – zabop
    Oct 20, 2022 at 11:05

1 Answer 1

2

Your data is a WMS so you need to use gdal rather than ogr but it will work, subject to the usual caveats about WMS producing a picture of the data not the actual data.

First I did a gdalinfo on your capabilities link:

gdalinfo "WMS:https://resource.sgu.se/service/wms/130/jorddjupsmodell"

which gave me:

Driver: WMS/OGC Web Map Service
Files: none associated
Size is 512, 512
Subdatasets:
  SUBDATASET_1_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=JORDDJUPSMODELL&CRS=CRS:84&BBOX=7.695846187558585,54.789792944659375,25.49767199163619,69.14432572671839
  SUBDATASET_1_DESC=Jorddjupsmodell
  SUBDATASET_2_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262
  SUBDATASET_2_DESC=Jorddjup 10x10m raster, skattat jorddjup till berg (m)
  SUBDATASET_3_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=JORDDJUPSMODELL_UNDERLAG&CRS=CRS:84&BBOX=7.695846187558585,54.789792944659375,25.49767199163619,69.14432572671839
  SUBDATASET_3_DESC=Utnyttjat underlag till jorddjupsmodellen
  SUBDATASET_4_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.JORDARTSKARTOR&CRS=CRS:84&BBOX=7.695846187558585,54.789792944659375,25.49767199163619,69.14432572671839
  SUBDATASET_4_DESC=Jordartskartor
  SUBDATASET_5_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.SPRICKZONER&CRS=CRS:84&BBOX=9.573739996397562,55.15360120621276,24.754326365873002,69.08179777642593
  SUBDATASET_5_DESC=Sprickzoner
  SUBDATASET_6_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.JORDDJUPSUPPGIFT_MINVARDE&CRS=CRS:84&BBOX=9.367378360212268,55.17706082150909,25.115002345782752,68.60472942460113
  SUBDATASET_6_DESC=Jorddjupsuppgift, minimi-värde
  SUBDATASET_7_NAME=WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.JORDDJUPSUPPGIFT&CRS=CRS:84&BBOX=9.166007780988775,55.17269811450856,25.38084453449201,69.1395851193348
  SUBDATASET_7_DESC=Jorddjupsuppgift
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

then I did a gdal_translate on the link I found in the results:

gdal_translate -of GTiff "WMS:https://maps3.sgu.se/geoserver/misc/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetMap&LAYERS=SE.GOV.SGU.MISC.JORDDJUPSMODELL.RASTER_INTERVALL&CRS=CRS:84&BBOX=7.395557961066897,54.844017949679454,26.329613471079828,69.41367203990262" -tr 10 10 /tmp/test.tif

to generate a tif, but I can't get that WMS to provide a map anywhere near Sweden.

2
  • Thank you! Gonna remove ogr2ogr references from my question not to confuse future readers. Also adding attempt relying on this answer.
    – zabop
    Oct 20, 2022 at 9:03
  • Question updated :)
    – zabop
    Oct 20, 2022 at 10:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.