Search code examples
pythonpandasmatplotlibcartopy

Scatter Plot on Projected Map is not in Correct Location

I have the trajectory of a hurricaneHurrican Florence correct location. In a "conventional" map CRS, the trajectory of the hurricane is correct

However, when I use a unconventional projection, the location of the hurricane is off. It starts East of North America and then goes West over the land. Why is that?

Hurrican Florence Wrong Location

Here is how I try it

data = {'Lat':{0:12.9,2:12.9,4:13.2,6:13.6,8:13.7,10:13.8,12:14.2,14:14.5,16:14.8,18:16.0,20:17.0,22:17.9,24:18.3,26:18.9,28:19.7,30:20.7,32:22.0,34:23.4,36:24.6,38:25.1,40:25.0,42:24.8,44:24.5,46:24.6,48:24.4,50:24.6,52:25.0,54:25.9,56:26.4,58:27.1,60:28.0,62:28.5,64:29.4,66:30.4,68:31.5,70:32.5,72:33.1,74:33.6,76:33.9,78:34.1,80:34.1,82:34.0,84:33.9,86:33.7,88:33.6,90:33.6,92:33.6,94:33.7,96:34.0,98:35.5,100:38.5,102:39.8,104:42.6},'Long':{0:-18.4,2:-19.4,4:-20.9,6:-21.4,8:-22.7,10:-24.7,12:-25.5,14:-26.7,16:-27.8,18:-30.2,20:-33.2,22:-35.9,24:-38.7,26:-41.0,28:-42.5,30:-43.9,32:-45.7,34:-47.2,36:-48.6,38:-49.8,40:-51.8,42:-53.2,44:-54.3,46:-55.2,48:-56.3,50:-57.7,52:-60.0,54:-62.4,56:-64.6,58:-66.2,60:-67.9,62:-69.5,64:-70.7,66:-71.8,68:-73.2,70:-74.3,72:-75.1,74:-76.0,76:-76.4,78:-77.2,80:-77.9,82:-78.4,84:-78.8,86:-79.3,88:-79.5,90:-79.8,92:-80.1,94:-80.8,96:-81.8,98:-82.1,100:-82.9,102:-80.2,104:-71.9}}

df = pd.DataFrame.from_dict(data)

# Invert longitude for plotting in western Hemisphere
df['Long'] = df['Long'] * -1

newCrs = ccrs.RotatedPole(pole_latitude=37.5, pole_longitude=177.5)
fig = plt.figure(figsize=(20, 10))
ax = plt.subplot(projection = newCrs)
ax.coastlines(resolution = '110m')
ax.scatter(df['Long'], df['Lat'], color = 'red', transform = newCrs)
Solution

  • The issue is that the transform keyword in ax.scatter should be the DATA CRS. So, assuming that your indata is fine with PlateCarree crs:

    df = pd.DataFrame(dict(Lat=[i for i in data['Lat'].values()], Long=[i for i in data['Long'].values()]))
 
newCrs = ccrs.RotatedPole(pole_latitude=37.5, pole_longitude=177.5)
dataCrs = ccrs.PlateCarree()

fig = plt.figure(figsize=(20, 10))
ax = plt.axes(projection = newCrs)
ax.set_global()
ax.coastlines()
ax.scatter(df['Long'], df['Lat'], color = 'red', transform = dataCrs)

    gives: plot showing alternative projection with data points correctly projected

    Please read this page explaining transforms for details and a tutorial.

    Alternative solution (recommending GeoPandas)

    Perhaps not a direct answer to your question, but if I do it through GeoPandas (highly recommended), I don't seem to get any problems:

    import pandas as pd
import geopandas as gpd
import contextily as ctx
import cartopy.crs as ccrs

data = {'Lat':{0:12.9,2:12.9,4:13.2,6:13.6,8:13.7,10:13.8,12:14.2,14:14.5,16:14.8,18:16.0,20:17.0,22:17.9,24:18.3,26:18.9,28:19.7,30:20.7,32:22.0,34:23.4,36:24.6,38:25.1,40:25.0,42:24.8,44:24.5,46:24.6,48:24.4,50:24.6,52:25.0,54:25.9,56:26.4,58:27.1,60:28.0,62:28.5,64:29.4,66:30.4,68:31.5,70:32.5,72:33.1,74:33.6,76:33.9,78:34.1,80:34.1,82:34.0,84:33.9,86:33.7,88:33.6,90:33.6,92:33.6,94:33.7,96:34.0,98:35.5,100:38.5,102:39.8,104:42.6},'Long':{0:-18.4,2:-19.4,4:-20.9,6:-21.4,8:-22.7,10:-24.7,12:-25.5,14:-26.7,16:-27.8,18:-30.2,20:-33.2,22:-35.9,24:-38.7,26:-41.0,28:-42.5,30:-43.9,32:-45.7,34:-47.2,36:-48.6,38:-49.8,40:-51.8,42:-53.2,44:-54.3,46:-55.2,48:-56.3,50:-57.7,52:-60.0,54:-62.4,56:-64.6,58:-66.2,60:-67.9,62:-69.5,64:-70.7,66:-71.8,68:-73.2,70:-74.3,72:-75.1,74:-76.0,76:-76.4,78:-77.2,80:-77.9,82:-78.4,84:-78.8,86:-79.3,88:-79.5,90:-79.8,92:-80.1,94:-80.8,96:-81.8,98:-82.1,100:-82.9,102:-80.2,104:-71.9}}
df = pd.DataFrame.from_dict(data)
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(x=df.Long,y=df.Lat), crs='wgs84')

newCrs = ccrs.RotatedPole(pole_latitude=37.5, pole_longitude=177.5)
gdf_proj = gdf.to_crs(newCrs)

ax = gdf_proj.plot(figsize=(20,8))
minx, miny, maxx, maxy = gdf_proj.total_bounds
ax.set_xlim(minx-50, maxx+75)
ax.set_ylim(miny-50, maxy+75)
ctx.add_basemap(ax, crs=gdf_proj.crs.to_string(), source=ctx.providers.Stamen.TonerLite)

    plot showing the correct points on the map using the alternative solution

    General tips/comments

    • You should not need to flip the Longitude coordinates at all.
    • Did you check if there are any issues with Matplotlib/Cartopy?
    • If you have the Proj4 definition string you do not even need cartopy if using GeoPandas. In this case the newCrs variable would become newCrs = '+ellps=WGS84 +a=6378137.0 +proj=ob_tran +o_proj=latlon +o_lon_p=0.0 +o_lat_p=37.5 +lon_0=357.5 +to_meter=111319.4907932736 +no_defs'.
    • You should probably clean up the use of Pandas DataFrame, storing a dictionary in a column of a dataframe is not really recommended, and probably makes your life more difficult.