2

I have data in the form of Section, Township, Range and meridian. I have 300,000 data points from all over the USA. Is there any way i can write R or python scripts to batch convert all the points to latitude and longitude? What kind of equations they use to convert it into latitude and longitude?

I found this online tool which gives latitude and longitude of centroid and four corners of the section, but I want to do the same thing using R or Python.

http://www.earthpoint.us/TownshipsSearchByDescription.aspx

I tried to search online but couldn't find any information regarding this except online tools where we put information for each point manually.

first few lines of my data:

Township    Range1  Section  QQ            Footages
154 N       96 W      3     SWNE       1980 FNL  1980 FEL
154 N       96 W      3     NWNE       660 FNL  1980 FEL
136 N       83 W     22     NENW       330 FNL  2310 FWL
141 N       81 W     18     SESE       660 FSL  660 FEL
133 N       75 W     35     NWSW       1980 FSL  660 FWL
132 N       48 W     11     NWSW       2310 FSL  330 FWL
157 N       85 W     16     SESE       60 FSL  300 FEL
140 N       77 W     6      SWSW       330 FSL  330 FWL
  • Is your data in Shapefile format? text file? – klewis Mar 1 '16 at 20:12
  • @klewis my data is in text file – liferacer Mar 1 '16 at 22:49
  • could you add the first few lines of your data to the question so we can see the format? – Spacedman Mar 3 '16 at 7:39
  • @Spacedman just added it – liferacer Mar 4 '16 at 2:18
  • I don't really have an answer, so I didn't put it in that section. I think you will have a tough time figuring this out with the data you have. If you look at PLSS outlines you'll find that they are not regular. This makes it hard to use a single equation to convert them to Lat/Long, especially if they are not in the same region. It maybe best to find shapefiles from which you can locate the centroids. You could also try some government agencies (e.g. BLM, state GIS clearinghouses) as they may have some of the data you'll need. – KevinB Mar 31 '16 at 16:12
1

You can send a text file to BLM for batch processing.
http://www.geocommunicator.gov/Geocomm/lsis_home/townshipdecoder/index.htm

We're looking for long answers that provide some explanation and context. Don't just give a one-line answer; explain why your answer is right, ideally with citations. Answers that don't include explanations may be removed.

  • Welcome to GIS SE! As a new user please take the tour to learn about our focused Question and Answer format. Please edit your answer with some more information about what actually needs to be done and how this answers the question. – Midavalo Oct 27 '16 at 17:43
  • Geocommunicator ceased functioning at the beginning of the year. Resources were moved, and your link currently redirects to the home page. It doesn't appear they take text files anymore, though maybe they do somewhere. All I can find is navigator.blm.gov/services which has a link to what is now a web service - you can still paste in, but there's also an API for submitting requests to the service. – Chris W Feb 8 '17 at 22:45
1

There is no formula/equation to convert from the PLSS to lat/lon. The PLSS is a reference grid, not a coordinate system (and note that it doesn't cover the entire US). To automate this, you would need to generate (or find - I don't know if one is out there) a coordinate lookup table for a script to reference or to use in a join.

Your lookup table would need to have a common field with your data (which may require modifying said data). At a minimum this would be the township-range-section-qq, but there's an added complication. Township and range are specified from reference baselines and meridians. In other words, there is more than one Township 1 North, Range 1 West. The example you posted doesn't give any method of identifying which baseline/meridian to use, and if you're covering such a large area you will most certainly run into more than one of them. One you had this table with every qq's centroid by name/reference in one column and the lat/lon for it in another (or others), you can join or add those coordinates to your data and then plot them.

As you've seen, it is possible to find the lat/lon of the corners/centroid for a given grid reference. From looking at your data, you would need to do this at the QQ (quarter-quarter) level. For many, but not all, areas covered by the PLSS, you can download a grid dataset to generate centroids and then coordinates from. But there's no single source for the entire system, and not all sources go past the section level. You'll also end up with multiple points in the same location, since anything in the same qq will have the same coordinate.

Alternatively you would have to find a service out there that can do this lookup for you. You've already found some tools that will do it one point at a time, but I believe some of them are also services (such as the BLM link in my comment on Dana's post). You have to learn the API and how to format your data for input, but it can be done in batch and not one at a time. Be aware most services like this have quantity limits that are either time or cost based.

1

I've had to do something similar and opted to use the BLM GIS Rest Service. You can create a URL and the service will respond back with the coordinates. Unfortunately you can only send about 25 PLSS descriptions at a time due to the constraints of how long a URL can be. So, I created a list of URLs and used wget to batch process them. Here is an example of a URL I created (WA330200N0280E0SN340 is the first PLSS description):

https://gis.blm.gov/arcgis/rest/services/Cadastral/BLM_Natl_PLSS_CadNSDI/MapServer/exts/CadastralSpecialServices/GetLatLon?trs=WA330200N0280E0SN340+%7C+WA330160N0260E0SN120+%7C+WA330170N0300E0SN170+%7C+WA330170N0280E0SN300+%7C+WA330160N0240E0SN130+%7C+WA330160N0250E0SN080+%7C+WA330160N0240E0SN130+%7C+WA330200N0280E0SN020+%7C+WA330210N0280E0SN340+%7C+WA330190N0250E0SN040+%7C+WA330200N0230E0SN040+%7C+WA330200N0280E0SN020+%7C+WA330210N0280E0SN340+%7C+WA330170N0300E0SN080+%7C+WA330180N0300E0SN240+%7C+WA330180N0300E0SN260+%7C+WA330200N0280E0SN020+%7C+WA330200N0280E0SN020+%7C+WA330190N0290E0SN280+%7C+WA330190N0290E0SN280+%7C+WA330190N0290E0SN330+%7C+WA330170N0240E0SN100+%7C+WA330170N0240E0SN110+%7C+WA330190N0270E0SN080+%7C+WA330200N0250E0SN110+%7C+&returnalllevels=&f=pjson&__ncforminfo=7x7RLQJqpDQ92T4PQI_-wSk34fhSp0krCuCkP5nOMe9rVpPf2stN7FAL0CQsbN_c62SdT15Xs7Rs5U4Xfe3epxz0dE12hq0HyWAsrAz430mfnjihcVff8g%3D%3D

After creating a list of these types of URLs and putting them into a text file, I ran this code to get my resultant list:

wget -i plss_working.txt -O out.txt

grep -A 2 "plssid" out.txt > new.txt

This is a bit cumbersome, but it worked!

  • this looks like the right answer. what's this PLSS format you are using? What'd this convert to: Township Range1 Section QQ Footages 154 N 96 W 3 SWNE 1980 FNL 1980 FEL? that's the first row in his example. – Dnaiel Sep 18 '18 at 15:31
0

I am also encountering the problem with converting the legal addresses (Section, Township, Range, Meridian) to latitude and longitude. I tried to implement the script to do this task based on some articles. Here is my Github link, here has the code for that. However, the code is the prototype version. It suffers from the accuracy of locating the position and can only query the location based on the first meridian. If you are interested in it, we can discuss more and figure out the solution. Thank you.

The method performs a calculation by using the offset of the township line and the range line. As mentioned in the here, the township line and the range line is 6 miles apart. I firstly extract the township number and the range number and then multiplying both of them by 6.0. It will yield a rough offset that how far a township line and a range line away from the start point (prime meridian & boundary of Canada and United States).

After that, we can get an approximated position of the township (at the right bottom corner of the township). And then, I extracted the Section number. In a township, there are 36 sections, and each section is 1-mile * 1-mile block. So, I created a numpy array which indicated how sections are arranged in a township. I can calculate an offset to the section from the start point of the township based on adding the indices of that array

I can then calculate an approximated bearing with python math.atan2(-range_offset, township_offset) if the location is at the west of the prime meridian. And convert the distance from miles to kilometers. Finally, I used the method in this repo to calculate the new latitude and longitude

  • for batch conversion. You can use df['legal_address'].map(function), to batch convert the address based on the address column in your dataframe. – Shijin Yang Oct 11 '18 at 17:44
  • It's better to download Cadnsdi or 24k grid and query it instead. The issues with the calculation listed above are,there are many Townships that are not 6x6 miles and many Sections that are not a sq. mile. Also, there are half Townships such as PM06, 92.5 West. Also, there are Indian Reservations inside a Meridian, such as Wind River in WY. – klewis Oct 12 '18 at 17:13
  • Hi klewis, thank you for your suggestions. Currently, I am also trying do a batch lookup. I have another dataset that contains the geospatial information, areas and the legal address information. I am exploring whether it will be a better method to match them with legal address information. – Shijin Yang Oct 12 '18 at 18:57

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.