2

I need to count the number of points (LAT, LON) from a cvs file that exist inside an arc defined by 2 azimuths. Ultimately it will need to be iterative where the arc defined by the 2 azimuths rotates a full 360 degrees counting the # of points in each arc. I have installed geopandas as a tool box.

3
  • The problem is that GeoPandas use shapely which uses euclidean distances. – gene Mar 9 '16 at 17:44
  • Can you recommend a better library? – Shonn McNeill Mar 9 '16 at 21:59
  • Search for GeoPy with pandas – gene Mar 10 '16 at 18:09
1
import matplotlib.pyplot as plt
import mpl_toolkits.basemap.pyproj as pyproj
import pandas as pd
import geopandas as gp
from shapely import geometry
import requests
import numpy as np
%matplotlib inline

# Data file path
csv_geo_25 = 
csv_geo_50 = 

# Read in data into Pandas dataframe remove NaN values where no weights were collected
# then renindex the dataframe where were NaN's were removed.
df_geo_25 = pd.read_csv(csv_geo_25)
df_geo_25 = df_geo_25.dropna(how='any').reset_index(drop=True)
df_geo_50 = pd.read_csv(csv_geo_50)
df_geo_50 = df_geo_50.dropna(how='any').reset_index(drop=True)

# Setup Proj class initialized for Alabama West, Zone 0102
myFrag = pyproj.Proj(init='EPSG:2760') #      http://spatialreference.org/ref/epsg/2760/, "Alabama West", Zone 0102

# Convert lat lon to easting northings for 2.5 lb TNT
easting_25, northing_25 = myFrag(df_geo_25['LON'].values,     df_geo_25['LAT'].values)
df_25 = pd.DataFrame.from_items([('Easting', easting_25), ('Northing',     northing_25)]).join(df_geo_25)

# Convert lat lon to easting northings for 5.0 lb TNT
easting_50, northing_50 = myFrag(df_geo_50['LON'].values, df_geo_50['LAT'].values)
df_50 = pd.DataFrame.from_items([('Easting', easting_50), ('Northing', northing_50)]).join(df_geo_50)


# Get ground zero to blast lat and lon
gz_lat =  34.59765667
gz_lon = -86.65476833

# Convert ground zero lat and lon to easting and northing
gz_easting, gz_northing = myFrag(gz_lon, gz_lat)

# Subtract ground zero easting and northing from fragment location data
df_25['X'] = df_25['Easting'] - gz_easting
df_25['Y'] = df_25['Northing'] - gz_northing
df_50['X'] = df_50['Easting'] - gz_easting
df_50['Y'] = df_50['Northing'] - gm_northing

# R
df_25['R'] = (df_25['X']**2 + df_25['Y']**2)**0.5
df_50['R'] = (df_50['X']**2 + df_50['Y']**2)**0.5

# Theta
df_25['Angle_Radians'] = np.arctan2(df_25['Y'], df_25['X'])
df_25['Angle_Degrees'] = df_25['Angle_Radians']*(180.0/np.pi)
df_25['Azmuth'] = np.where(df_25['Angle_Degrees']>=0.0,     df_25['Angle_Degrees'], 360.0+df_25['Angle_Degrees'])
df_50['Angle_Radians'] = np.arctan2(df_50['Y'], df_50['X'])
df_50['Angle_Degrees'] = df_50['Angle_Radians']*(180.0/np.pi)
df_50['Azmuth'] = np.where(df_50['Angle_Degrees']>=0.0,  df_50['Angle_Degrees'], 360.0+df_50['Angle_Degrees'])

# Plot fragmen locations and scale point by mass size
plt.figure(figsize=(8, 8))
ax = plt.axes(polar=True)
ax.set_theta_zero_location('E')
ax.scatter(df_25['R'], df_25['Azmuth'], s=df_25['WT_g'])
ax.scatter(df_50['R'], df_50['Azmuth'], s=df_50['WT_g'], color='r')

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.