Draw circle with longitude, latitude and radius (km) in cartopy of python

user_6781624 0
Viewed 746
#!/usr/bin/env python

import os, sys
import pandas as pd
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
import numpy as np
from cartopy.geodesic import Geodesic


if __name__ == '__main__':

    stn = pd.read_csv('obs_station.csv')
    gd = Geodesic()

    lcc = ccrs.LambertConformal(central_longitude=126., central_latitude=38.)
    fig = plt.figure(figsize=(7,7))
    ax = fig.add_subplot(111, projection=lcc)
    ax.coastlines(resolution='50m')
    geoms = []
    for lon, lat in zip(stn['longitude'], stn['latitude']):
        cp = gd.circle(lon=lon, lat=lat, radius=250000.)
        geoms.append(sgeom.Polygon(cp))
    ax.add_geometries(geoms, crs=lcc, edgecolor='r')
    ax.set_extent([120., 133., 30., 43.])
    plt.show()

The file 'obs_station.csv' contain several coordinates of longitudes and latitudes.

Using code above, I try to draw circles with specific radius (250 km). But, nothing is on the map as below. Only show the map with coastlines.

I don't know what is the problem. Help please.

Result

cartopy geometry python
1 Answers

You did not get the plots of the circles because of wrong coordinate transformation you specifies in .add_geometries() statement.

To get it right, suppose I use this data file: 'obs_station.csv':

longitude,latitude
127.603897,36.932988
126.505337,38.555939

And the modified code:-

#import os, sys
import pandas as pd
import cartopy
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import shapely.geometry as sgeom
import numpy as np
from cartopy.geodesic import Geodesic

if __name__ == '__main__':

    stn = pd.read_csv('obs_station.csv')
    gd = Geodesic()

    # This is long-lat coordinate system for use in ..
    # .. coordinate transformation options
    src_crs = ccrs.PlateCarree()

    lcc = ccrs.LambertConformal(central_longitude=126., central_latitude=38.)
    fig = plt.figure(figsize=(7,7))
    ax = fig.add_subplot(111, projection=lcc)
    ax.coastlines(resolution='50m')
    geoms = []
    for lon, lat in zip(stn['longitude'], stn['latitude']):
        cp = gd.circle(lon=lon, lat=lat, radius=250000.)

        #x,y = lcc.transform_point(lon, lat, src_crs)
        #cp = gd.circle(lon=x, lat=y, radius=250000.)
        geoms.append(sgeom.Polygon(cp))

    # Note the specification of coordinate transformation, using the
    # .. correct parameter:  crs=src_crs
    ax.add_geometries(geoms, crs=src_crs, edgecolor='r', alpha=0.5)
    ax.set_extent([120., 133., 30., 43.])

    plt.show()

The output:

korea01

In conclusion, your line of code:

ax.add_geometries(geoms, crs=lcc, edgecolor='r')

needs correct CRS. And the correct CRS is

ccrs.PlateCarree()
user_21774130
Related