0

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')

Read in data

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

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?

1
  • 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.
    – Vince
    Commented May 6, 2022 at 0:36

2 Answers 2

3

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

import geopandas
import pandas as pd

df = pd.read_csv("data_for_question.csv").set_index("Sample_ID")
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.

1
  • 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
1

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!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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