0

I am trying to draw some "station" coordinates on the NYC map. I use Geopandas(gpd) to create a GeoDataFrame and create a NYC map from gpd.datasets with EPSG of 4326.

In addition, I wanted to add a background image to my map following the Geopandas tutorial (https://geopandas.org/gallery/plotting_basemap_background.html). However, this works on EPSG of 3857. Is there any way to maintain the information and change the EPSG?

Here is my codebase

# Creating GeoDataFrame
gdf = gpd.GeoDataFrame(stations, geometry=gpd.points_from_xy(stations.Longitude, stations.Latitude))
# Find the NYC Borough map from contextily pkg
nyc = gpd.read_file(gpd.datasets.get_path('nybb')).to_crs(epsg=4326)
ax = nyc.plot(figsize=(10,10), alpha=0.3, edgecolor="k")

# Contextily adding background image
def add_basemap(ax, zoom, url='http://tile.stamen.com/terrain/tileZ/tileX/tileY.png'):
    xmin, xmax, ymin, ymax = ax.axis()
    basemap, extent = ctx.bounds2img(xmin, ymin, xmax, ymax, zoom=zoom, url=url)
    ax.imshow(basemap, extent=extent, interpolation='bilinear')
    # restore original x/y limits
    ax.axis((xmin, xmax, ymin, ymax))

add_basemap(ax, zoom=10, url='http://tile.stamen.com/toner-lite/tileZ/tileX/tileY.png')
# nyc.to_crs(epsg=4326)
gdf.plot(ax=ax, color="red")
ax.set_axis_off()
Maeror
  • 93
  • 2
  • 11
  • You may refer to the following answer. https://gis.stackexchange.com/questions/48949/epsg-3857-or-4326-for-googlemaps-openstreetmap-and-leaflet To summarise: "So if you are making a web map, which uses the tiles from Google Maps or tiles from the Open Street Map web service, they will be in Spherical Mercator (EPSG 3857 or srid: 900913) and hence your map has to have the same projection." – Debjit Bhowmick Mar 24 '20 at 05:45

1 Answers1

0

Yes there is. You can usually pass the EPSG information of your data to the add_basemap() function of contextily, by

add_basemap(ax, crs=nyc.crs.to_string(), zoom=10, url='...')

The only thing required is, that you have set a proper CRS on your data.

In your case you probably should add the contextily.wrap_tiles() function.

tr_basemap, tr_extent = ctx.wrap_tiles(basemap, extent, t_crs='EPSG:4326')

This should give you, what you want.

For further reading you can look here:

https://contextily.readthedocs.io/en/latest/warping_guide.html#Convert-the-tiles-to-your-data%E2%80%99s-CRS https://contextily.readthedocs.io/en/latest/reference.html

J. P.
  • 306
  • 1
  • 8