# Convert from latitude/longitude into feet and acres using Python

I have a field that I am trying to convert into acreage using Python. My understanding is that the process depends on your location in the world, as this process can be different based on location. The coordinates I have are for data that were collected in the field. I have interpolated them, and am trying to determine the field acreage size.

What I have done involves obtaining the lon/lat max and min, subtracting them, and multiplying each by a factor of 110,000. Then, I divided this answer by the number of feet in a meter (0.3048).

Is this the correct process?

Note: For reproducibility, the input file can be downloaded from here (file size: 1 KB). The code is set up to run in a Jupyter Lab/Notebook file.

## Here is the code:

Import modules

``````import numpy as np
import pandas as pd
from scipy.interpolate import griddata, interp2d
import holoviews as hv
hv.extension('bokeh')
``````

``````df = pd.read_csv('data_for_question.csv')
``````

Generate grid

``````#Minimum and maximum longtude values
lon_min = df['longitude'].min()
lon_max = df['longitude'].max()

#Create range of nitrogen values
lon_vec = np.linspace(lon_min, lon_max, 100) #Set grid size

#Find min and max latitude values
lat_min = df['latitude'].min()
lat_max = df['latitude'].max()

# Create a range of N values spanning the range from min to max latitude
# Inverse the order of lat_min and lat_max to get the grid correctly
lat_vec = np.linspace(lat_max,lat_min,100,)

# Generate the grid
lon_grid, lat_grid = np.meshgrid(lon_vec,lat_vec)
``````

Specify points and targets to interpolate over

``````points = (df['longitude'], df['latitude'])
targets = (lon_grid, lat_grid)
map_bounds=(lon_min,lat_min,lon_max,lat_max)
``````

Generate map

``````grid = griddata(points, df['n_ppm'], targets, method='cubic', fill_value=np.nan)

interp_map = hv.Image(grid, bounds=map_bounds).opts(aspect='equal',
colorbar=True,
title='Graph Title',
xlabel='Longitude',
ylabel='Latitude',
cmap='Reds')

interp_map
``````

Compute pixel area

``````lat_span = (lat_max - lat_min) * 111_000 / 0.3048
lon_span = (lon_max - lon_min) * 111_000 / 0.3048

# Convert to acres
field_area = (lat_span * lon_span) / 43560
cell_area = field_area / lon_grid.size

field_pix = (~np.isnan(grid)).sum()

print("Total grid size is:", (round(field_area, 2)), "acres")
print("Total field size is:", (round(cell_area * field_pix, 2)), "acres")
``````

I happen to know that these values are close to the actual field values, but is there a better way to do this?

• Hard-coding the conversion of degrees to distance is pretty much always going to be wrong. Doing it right involves spheroidal trigonometry and a partial differential equation only solvable by iterative means. I suggest using a code library. Commented May 6, 2022 at 0:36

Use GeoPandas (and pyproj) to transform these to a better CRS:

``````import geopandas
import pandas as pd

lat_lng_gdf = geopandas.GeoDataFrame(
df,
geometry=geopandas.points_from_xy(
df.longitude, df.latitude, crs="EPSG:4326"))
xy_gdf = lat_lng_gdf.to_crs(lat_lng_gdf.estimate_utm_crs())
print(repr(xy_gdf.crs))
# <Derived Projected CRS: EPSG:32614>
# Name: WGS 84 / UTM zone 14N
# Axis Info [cartesian]:
# - E[east]: Easting (metre)
# - N[north]: Northing (metre)
# Area of Use:
# - name: Between 102°W and 96°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Manitoba; Nunavut; Saskatchewan. Mexico. United States (USA).
# - bounds: (-102.0, 0.0, -96.0, 84.0)
# Coordinate Operation:
# - name: UTM zone 14N
# - method: Transverse Mercator
# Datum: World Geodetic System 1984 ensemble
# - Ellipsoid: WGS 84
# - Prime Meridian: Greenwich

xmin, ymin, xmax, ymax = xy_gdf.total_bounds
width = xmax - xmin
height = ymax - ymin
area = width * height
print(f"Width: {width:.2f} m, {width / 0.3048:,.0f} ft")
# Width: 693.47 m, 2,275 ft
print(f"Height: {height:.2f} m, {height / 0.3048:,.0f} ft")
# Height: 696.42 m, 2,285 ft
print(f"Area: {area:.1f} m^2, {area / 4046.8564:.2f} acres")
# Area: 482,943.3 m^2, 119.34 acres
``````

As you can see, the points are all in a small region, and could easily fit within a UTM zone. Having length units in meters (or feet) should solve your questions.

• Thanks! That is very helpful. As a beginner in geopandas, I am increasingly amazed with what you can do with it. Commented May 9, 2022 at 14:27

You will face a challenge in that the distance spanned (in ground coordinates, e.g. metres) by one degree of longitude is not constant. One degree longitude near the poles is vanishingly small. One degree longitude spans the same distance as one degree latitude only at the equator.

As an approximation, if you assume the earth is spherical, then the distance spanned by one degree is `cos(latitude) * 111,000` (in metres), so at least you need to adjust your area calculation by this amount. If your latitude span is small, i.e. you are measuring distances in the metres or hundreds of meters, then `cos` of any latitude in the range `lat_min` to `lat_max` will be about the same and it's a simple scaling. If you're spanning several degrees, you'll have to reflect that a "rectangle" in lat-long will have different north and south side lengths.

If you want to be more precise, I would not reinvent the wheel but use the `pyproj` library to convert to, or work in, projected coordinates in metres. That library will take care of it (if you use an appropriate projection), as well as not make the simplification that the earth is fully spherical.

As to what is the best projection, if your objects are pretty small, your local UTM zone is probably a sensible choice. If your objects are more than a couple of degrees in lat/long, you'll have to find an appropriate equal-area projection. Do not use "Google Mercator" EPSG:3857, which distorts scales in a way inimical to area calculations!