1

I am trying to plot the Pumps.shp data on top of the OSMap.tif file from this website on the same figure.

I tried using rasterio.plot() and geopandas.plot() methods, with matplotlibs subplots.

The problem is that the plots don't match, the raster file gets plotted in the range(0,1000) for both axis and the shp gets plotted in the actual coordinates range(around 50000 on the x axis and around).

The crs are equal in both objects and the coordinates are in the same range. Why is this? What am I doing wrong?

Here is my code

   import rasterio as rast
   import rasterio.plot as rsplot
   import geopandas as gpd
   src=rast.open("OSMap.tif")
   data=gpd.read_file("Pumps.shp")
   fig,ax=plt.subplots()
   rsplot.show(src,ax=ax)
   data.plot(ax=ax)
   plt.show()

This is the result of calling src.bounds:

BoundingBox(left=528765.0, bottom=180466.0, right=529934.0, top=181519.0)

This is the result of data.bounds

(528765.0, 180466.0, 529934.0, 181519.0)

This is crs of both:

CRS({'lon_0': -2, 'y_0': -100000, 'k': 0.9996012717, 'lat_0': 49, 'proj': 'tmerc', 'wktext': True, 'datum': 'OSGB36', 'no_defs': True, 'x_0': 400000, 'units': 'm'})

Absalon Castañon
  • 359
  • 5
  • 18
  • I can't find the two input files on the page you link to. – ImportanceOfBeingErnest Dec 06 '17 at 18:13
  • It's on the first zip, where it says : A zip file with the Vector data as Shapefiles... @ImportanceOfBeingErnest – Absalon Castañon Dec 06 '17 at 19:31
  • 1
    https://stackoverflow.com/questions/34458251/plot-over-an-image-background-in-python/34459284#34459284 – Paul H Dec 06 '17 at 20:51
  • What version of `rasterio` are you using? Not that this seems particularly likely to be the problem, but `rasterio.plot.show` only began using the extent of the raster in 0.32. – jdmcbr Dec 07 '17 at 03:42
  • @jdmcbr I am using version 0.36 – Absalon Castañon Dec 07 '17 at 04:00
  • @AbsalonCastañon I downloaded the example data and was able to plot it all correctly. Can you confirm that `with_bounds=True` (which changes the plotting extent to match the raster extent) is the default in `rasterio.plot.show` for you? – jdmcbr Dec 07 '17 at 06:05

2 Answers2

1

I had the same problem with rasterio 0.36.0. I first tried to translate and scale the raster but than prefered to translate the shapefile.

My code looks like:

import geopandas as gpd
import matplotlib.pyplot as plt
import rasterio

image = rasterio.open('input.tif') # with tgw world file
shapefile = gpd.read_file('input.shp')

# coordinates and scaling factors
scale_x = image.transform[1]
scale_y = image.transform[5]
x0 = image.transform[0]
y0 = image.transform[3]

# translates back shapefile
shapefile.geometry = shapefile.translate(-x0, -y0)
shapefile.geometry = shapefile.scale(-1.0/scale_x, -1.0/scale_y, origin=(0, 0, 0))

# plots both elements
fig, ax = plt.subplots()
ax = rasterio.plot.show(image.read(), with_bounds=True, ax=ax)
shapefile.plot(ax=ax)
kogexo
  • 53
  • 6
0

Use matplotlib imshow instead of rasterio show. Pass the bounds of raster as "extent" parameter of imshow.

cass
  • 1