2

Can I read a raster map which has been constructed in mapserver as WMS with rasterio?

rasterio.open

('http://...&service=WMS&request=GetMap&...FORMAT=image/tiff')

I tried :

from rasterio.io import MemoryFile

def print_metadata(dataset):
     print(dataset.profile)


url ='http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

tif_bytes = open(url,'rb').read()

with MemoryFile(tif_bytes) as memfile:
     with memfile.open() as dataset:
         print_metadata(dataset)

and it gives :

FileNotFoundError: [Errno 2] No such file or directory: 'http://localhost:8080/?map=/maps/ivm3.map&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=tr_ivm_gray&STYLES=&SRS=epsg:3857&BBOX=2857613.741389,4274927.875099,4989229.633477,5176940.449967&WIDTH=1160&HEIGHT=540&FORMAT=image/tiff'

another code :

raster=rasterio.open(url)
print (type(raster))

print ( raster.meta )
print ( raster.count )

array=raster.read()

this gives:

<class 'rasterio.io.DatasetReader'>
{'driver': 'WMS', 'dtype': 'uint8', 'nodata': None, 'width': 1073741824, 'height': 454363579, 'count': 3, 'crs': CRS({'init': 'epsg:3857'}), 'transform': Affine(0.001985222000710666, 0.0, 2857613.741389,
       0.0, -0.0019852220040462356, 5176940.449967)}
3
Traceback (most recent call last):
  File "rio1.py", line 20, in <module>
    array=raster.read()
  File "rasterio/_io.pyx", line 322, in rasterio._io.DatasetReaderBase.read
MemoryError

raster.read() tries to read too much data. What can I do here?

  • from rasterio.io import MemoryFile def print_metadata(dataset): print(dataset.profile) url ='localhost:8080/?map=/maps/…' #tif_bytes = open('deptiled.tif', 'rb').read() tif_bytes = open(url,'rb').read() with MemoryFile(tif_bytes) as memfile: with memfile.open() as dataset: print_metadata(dataset) >> gives FileNotFoundError: ..such file or directory – atemiz88 Aug 28 '18 at 6:46
  • 1
    'width': 1073741824, 'height': 454363579 you're trying to read way too much data. – user2856 Aug 28 '18 at 8:22
0

I had the same problem, and I solved it mostly thanks to this question.

My XML file looks like this:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>http://a.tiles.wmflabs.org/bw-mapnik/${z}/${x}/${y}.png</ServerUrl>
        <SRS>EPSG:3857</SRS>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>18</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <Cache>
        <Path>../.gdalwmscache</Path>
        <Extension>.png</Extension>
    </Cache>
</GDAL_WMS>

But any GDAL_WMS file should do. Note the TileLevel parameter, this sets the highest zoom level you want to go to. So even if this tileset supported zoom level 20, GDAL/rasterio will only see as far as 18.

The basic Python code I have is:

TMS_dataset = rasterio.open('templates/gdal_tms.xml')
# Calculate the pixel positions of the desired bounding box at the highest zoom level
# as specified in the XML file.
bl = TMS_dataset.index(west, south, op=math.floor)
tr = TMS_dataset.index(east, north, op=math.ceil)
# image_size is a tuple (h, w, num_bands)
output_dataset = np.empty(shape=image_size, dtype=TMS_dataset.profile['dtype'])

# Read each band
TMS_dataset.read(1, out=output_dataset[:, :, 0], window=((tr[0], bl[0]), (bl[1], tr[1])))
TMS_dataset.read(2, out=output_dataset[:, :, 1], window=((tr[0], bl[0]), (bl[1], tr[1])))
TMS_dataset.read(3, out=output_dataset[:, :, 2], window=((tr[0], bl[0]), (bl[1], tr[1])))

# Create an output image dataset
output_image = rasterio.open('flange3.png', 'w', driver='png', width=image_size[1], height=image_size[0], count=3,
                             dtype=output_dataset.dtype)
# Write each band
output_image.write(output_dataset[:, :, 0], 1)
output_image.write(output_dataset[:, :, 1], 2)
output_image.write(output_dataset[:, :, 2], 3)

output_image.close()

This particular code doesn't bother with georeferencing the output because I didn't need it for this application, but the rasterio docs tell you how to achieve that if you need it. It will emit a NotGeoReferencedwarning and will complain that the output file doesn't exist, but they are harmless as far as I can tell.

Note that it doesn't read in every tile at every level, GDAL/rasterio is clever enough to read just the tiles that meet the desired resolution. In addition, adding the Cache option to the XML caches tiles locally. This is only valid for TMS-like datasets, but is absolutely perfect for my use-case.

In addition, despite the docs saying that you can't resample with "decimated reads", you can add resampling=Resampling.cubic to the dataset read() function and it does the right thing. Although this might only apply to this use-case.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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