1

I downloaded a raster image from the following address, which corresponds to rivers and streams in the Barcelona area.

https://geoserveis.icgc.cat/icgc_geotecnicbcn25m/wms/service?REQUEST=GetMap&SERVICE=WMS&VERSION=1.3.0&LAYERS=RIER_LN&STYLES=&FORMAT=image/tiff&BGCOLOR=0xFFFFFF&TRANSPARENT=TRUE&CRS=EPSG:25831&BBOX=415134.5624,4571467.0516,444192.4298,4589092.5720&WIDTH=1634&HEIGHT=990

However, the coordinates are absolutely wrong, as the bbox has turned into (0,0,17,10.3125) instead of the (415134.5624,4571467.0516,444192.4298,4589092.5720) I was looking for. When making a quick plot with tmap, it appears correct (not deformed), exactly the same map than here. target image

Hence, is there any way I can change the coordinates, even if it is by simple direct proportionality?

1 Answer 1

2

I have no idea why the extent and CRS are both horribly wrong in the TIFF.

You can set the extent and coordinate system manually

> r = rast("download.tiff")
> crs(r) = "epsg:25831"
> ext(r) = c(415134.5624,4571467.0516,444192.4298,4589092.5720)[c(1,3,2,4)]

Note how I have pasted the extent from the request URL but have rearranged it using a subset vector because its not in the right order for the ext function.

Now, in theory you should be able to run tmap on that and get it. But:

> tm_shape(r) + tm_rgb()
stars object downsampled to 1284 by 778 cells. See tm_shape manual (argument raster.downsample)
Error in CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  : 
  OGRCreateCoordinateTransformation(): transformation not available
In addition: Warning message:
In CPL_transform(x, crs, aoi, pipeline, reverse, desired_accuracy,  :
  GDAL Error 6: Cannot find coordinate operations from `ENGCRS["unnamed",EDATUM[""],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]' to `EPSG:4326'
> 

Something in tmap is not happy with trying to transform to lat-long. Possible stars? If we try using rast's reprojection to lat-long so that tmap doesn't have to:

> rp = project(r, "epsg:4326")                                         
> tm_shape(rp) + tm_rgb()

then we get:

enter image description here

which is definitely over Barcelona... The lines do seem to line up with some features on the ground - usually roads though... maybe this is drainage of some sort... Anyway... If I load into QGIS and use "multiply" blend to see through the white...

enter image description here

1
  • Thanks! This worked great. You're not supposed to see them because they run through the sewers, but here's where they are supposed to be.
    – marcbosch
    Commented Oct 29, 2022 at 15:19

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.