Create polygons for sites table in BETYdb

[dlebauer]

Description

The spring planting plots are described in this spreadsheet on google drive

• the plot layout is described in the 'experimental plot plan' table, the important bits are:

dimension	value
# rows	32
# rows / plot	2
# plots (2 rows ea)	864
# ranges	54
# columns	16
row width (m)	0.762
plot length (m)	4
row length (m)	3.5
alley length (m)	0.5

I began working with these data in R. The following code uses the above dimensions to generate the points and then convert them to polygons in WKT that would be suitable for entering into the database
(20160421_plots.csv.txt)
except that we need to

• put them in gantry coordinate system and
• convert them to a geographic coordinate system.
• insert them into the database
  • ultimately we need something like the follw
    The integer values for 'plot' are 1-864 are associated with plots in the terraref instance of BETYdb named "MAC Field Scanner Field Plot i"

Here are the sites in TERRA-REF BETYdb that we are trying to update: https://terraref.ncsa.illinois.edu/bety/sites.json
Here are the sql statements in a system that puts (0,0) at the bottom right of the field:
insert_site_geometries.sql.txt

# code in https://github.com/terraref/reference-data/blob/master/scripts/maricopa_field_metadata.Rmd
require(udunits2)
library(dplyr)
library(tidyr)
library(ggplot2)

x_origin = 0 
y_origin = 0
length = 3.5
alley = 0.5
row_spacing = 30
ranges = 54
columns = 16

row_spacing_cm <- ud.convert(row_spacing, 'in', 'm')

dx <- 2 * row_spacing_cm
dy <- length

x <- x_origin + (1:columns - 1) * dx
y <- y_origin + (1:ranges - 1) * (length + alley) 

grid <- data.frame(expand.grid(x,y)) %>% 
  transmute(plot = 1:(ranges*columns), x1 = Var1, y1 = Var2, x2 = x1 + dx, y2 = y1 + dy) #%>% 
  #tidyr::gather(key = 'lat', )

theme_set(theme_bw())
ggplot(data = grid) + 
  geom_text(aes(x1+dx/2, y1+dy/2, label = plot), size = 2) +
  geom_point(aes(x1,y1), size = 1, alpha = 0.4) + 
  geom_point(aes(x2,y2), col = 'red', size = 1, alpha = 0.4) +
  ggtitle('gantry plot layout') + xlab('x') + ylab('y')


grid2 <- grid %>% gather(key, value = x, starts_with("x")) %>% gather(key2, value = y, starts_with("y"))   %>% select(-starts_with("key"))  %>% arrange(plot) %>% setDT

library(spatstat)
library(wellknown)
write.csv(merge(grid, plots, by = 'plot'), file = '20160421_plots.csv.txt')

insert <- plots[,list(sql = paste0('update sites set geometry = ', geometry, 'updated_at = now(), where name = ', paste0("'MAC Field Scanner Field Plot ", plot, "';")))]

writeLines(insert$sql, con = file('insert_site_geometries.sql.txt'))

Related issues:

• Converting lemnatec scanner metadata to absolute geographic positions reference-data#32
• Convert gantry coordinate system from right handed (x positive to left) to conventional (x positive to right) #109
• Create custom spatial reference systems for PostGIS reference-data#51

[dlebauer]

@JeffWhiteAZ or @smarshall-bmr

1. could you please provide information for mapping these plots to the gantry (0,0) coordinates?
2. do you have similar dimensions for the second planting?

[smarshall-bmr]

https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=883764630

https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=1433617856

These are the plot locations as given in the gantry coordinate system for the first and second planting respectively. The y coordinates change a negligible amount as to the point that the y coordinates are given as the same between plantings. The x axis is offset to the north approximately 50cm for the second planting.

[tingli3]

@dlebauer I have a few questions related to the R code:

1. Q1: the plot number in the picture given by ggplot2 in the R code does not match with the excel below (is it important?)
   https://docs.google.com/spreadsheets/d/1QQaWc0UaQQKfEtnSO1G2za8tKU2huC0_VYMBqm5CKAo/edit
   For example, the location of plot 17 is not consistent with the excel.

A1 yes, this is important. Please use the spreadsheet; there is probably a bug in my code

1. Q2: what is the correct lat lon for the four corners that we should build the projection on? I have designed a projection system using the lat lon of the corners provided before. That is

S-E 33.07459197 -111.97480975
S-W 33.07459172 -111.97504692
C-W 33.07547558 -111.97504675
C-E 33.07547550 -111.97481017
N-W 33.07642747 -111.97504728
N-E 33.07642819 -111.97481064

However, I found that the corners points I used are actually the six points SW SE CW CE NW NE in the sheet below. I am not sure what is the relationship between the three lat&lons in the sheet and which one is the truth.
https://docs.google.com/spreadsheets/d/1pTgmpf9kltYW1SHo_K51cynVTYs6-JdS6rt4usfwV10/edit

A2 We are assuming USDA values are the truth. They are within ~2cm of reality.

1. Q3: why the bounding box is much larger. The bounding box I got from the grid calculated in the R code is much larger than the bounding box of the corner points in the google sheet mentioned in Q2.

boundingbox=c(min(grid$x1),min(grid$y1),max(grid$x2),max(grid$y2))
boundingbox
[1] 0.000 0.000 24.384 215.500

A3 By 'much larger' do you mean 24 x 215m in the R code instead of 22 x 204 m in the reference points? This difference is because the planted area is bigger than the area that the sensors travel (there are border rows that the sensors don't see, but that the tractor and the UAV will see. So the purpose of those six points is to map the gantry and MAC coordinate reference systems.

1. Q4: I saw many people talking about gantry coordinates and lemnatec coordinates. I wonder what's the relationship between these different coordinates. Also, how they are related to coordinates with SW corner as (0,0) .

A4: The terms 'Lemnatec', 'Gantry', 'Field Scanner' and sometimes 'Wagner Mueller' are used to refer to the same (x, y) coordinate system with x decreasing going west.

1. I failed to run the R code from grid2 <- grid %>% gather(key, value = x, starts_with("x")) %>% gather(key2, value = y, starts_with("y")) %>% select(-starts_with("key")) %>% arrange(plot) %>% setDT to the end. I am not sure what that part means.

A5 don't worry about this. all you need is in the file attached 20160421_plots.csv.txt. The code is just for reference.

Here is an example of how to add lat1 and lon1 for X1 and Y1, add lat2 and lon2 for X2 and Y2 to grid data frame using the projection I provided before.

require(proj4)
options(digits = 12)
x1y1=cbind(grid$x1,grid$y1)
x2y2=cbind(grid$x2,grid$y2)
proj= '+proj=tmerc +lat_0=0 +lon_0=-111.843479 +k=1 +x_0=12285.092664 +y_0=-3661028.344732 +datum=WGS84 +units=m +no_defs '
latlon1=project(x1y1, proj, inverse = T, degrees = TRUE, silent = FALSE, ellps.default="sphere")
latlon2=project(x2y2, proj, inverse = T, degrees = TRUE, silent = FALSE, ellps.default="sphere")
grid$lon1=latlon1[,1]
grid$lat1=latlon1[,2]
grid$lon2=latlon2[,1]
grid$lat2=latlon2[,2]

A6 that is very helpful! I've amended it to the documentation

[dlebauer]

@tingli3 I answered your questions inline above

[tingli3]

@dlebauer Thank you.

The true SW and SE lat & lon will be very important in the projection definition. If USDA is the truth lat & lon for the SW and SE, I need to update the projection based on it and make sure after projection the SW will be (0,0) and SE will be something like (27.***, 0). Please tell me whether my understanding is correct.

However, here the coordinate system we generated using R code start at (0, 0) and ends at (24.384, 215.500). The four points given by USDA is approximately 27.5m * 110m. If that is the case, I wonder whether the four corners of the gantry plot layout should match with the 4 points of USDA. Specifically, whether the SW point of the plot1 is the SW point of USDA points and SE point of the plot16 is the SE point of USDA points, and how to handle the difference.

[dlebauer]

@tingli3 there are two different issues.

First, is to determine the coordinate reference systems. That should be covered in issue terraref/reference-data#32. There are four (or more) points used to align the two (MAC and USDA) GPS systems. These are in concrete posts around the edge of the field. There are also points used to align the gantry coordinate system with the MAC RTK GPS.

Second is this issue, to determine the bounding boxes of the individual plots in these different coordinate reference systems.

I will send you an email with a complete discussion of how these points were collected.

[tingli3]

@dlebauer So, what is the X and Y in the R code corresponds to ? Where is the true (0,0) ?

[dlebauer]

@tingli3 First, I apologize on behalf of all involved that this is such a mess.

But to try and clarify, there are a three (X,Y) coordinate systems

1. the gantry / lemnatec system that has reversed X. These are
2. the 'corrected' X,Y with (0,0) on the SW reference point. This does not yet exist, but is covered in (Convert gantry coordinate system from right handed (x positive to left) to conventional (x positive to right) #109) that this would be helpful for plotting and because humans are used to X increasing to the right / East
3. the values used in the R code and experimental design. These are in relation to the corner of the planted area, which will change each season.

To map the plots (3) to the gantry (1), it is necessary to align the X,Y pairs in the R code with the X,Y measured coordinates provided by @smarshall-bmr in this spreadsheet: https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=883764630

[tingli3]

@dlebauer
So, the first step is to convert the gantry system based on its bounding box to make sure that X starts at 0 and increases to the East using something like x'=-y+22.135 and y'=x-3.8. But this assumes that the six point given is the bounding box (22.135m*203.5m) of the gantry field. It is important to know the exact boundary of the gantry system and how the plots are placed in it.

If this is not the case, I think the first step is to find out the bounding box of gantry field , and how the gantry should be converted (set 0,0 X-axis Y-axis for it), and how it is relate the to the XY of the plots. Is it possible to use the min and max in the https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=883764630 to get bounding box of the gantry coordinates? If this is the case, which sheet (RangeCentersSorghum1 or RangeCentersSorghum2) should be used? And it seems that range 1 and range 54 does not have value in RangeCentersSorghum1 and the value in this two sheet are different.

If it is hard to know the exact bounding box of gantry, will it possible for us to define a coordinate for it? That is, the (0,0) of the plots may related to a (small value, small value) of the defined coordinate with X increasing to the East.

[dlebauer]

@tingli3 I think it would help if you could create a map with the different coordinates reference systems and available reference points. I am confused that x'=f(y) instead of f(x), and it would really help to clarify. Otherwise, I am afraid that the discussion we are having here will add to the confusion.

... and as you sort this out please document the issues and solutions terraref/documentation#7

[dlebauer]

@tingli3 RangeCentersSorghum1 and RangeCentersSorghum2 are for the first and second planting seasons, respectively.

I believe the GPS reference points that compare UA MAC and USDA GPS's correspond to concrete pylons that are in the ground a few meters outside of the planted area.

[tingli3]

@dlebauer
Using x'=f(y) instead of f(x) is to transpose X and Y as the X goes S to N in the gantry coordinates (the X and Y of the gantry should be inversed).

For the two planting seasons, are the coordinates recorded in the same gantry coordinate system? Which one should be used? Finding the exact bounding box of the gantry data is very important.

[dlebauer]

@tingli3

For the two planting seasons, the coordinates are recorded in the same system. We will need to recreate the plot boundaries each time a new crop is planted, and the bounding box of the planted area as well as the plots themselves will change.

[tingli3]

@dlebauer
From the information, I assume we are going to use one system and put all the data with different bounding box into it. Also, the two planting seasons at least must have to have the same (0,0) point (though we do not know where it is ) and XY axes (I think this should be true). That means the same gantry coordinate in both two seasons should exactly point to the same point. Is this correct?

And, as the bounding box is changing and (0,0) point is unknown, the new gantry coordinates with X increasing to the east will also not start exactly at (0,0). The (0,0) point of the plots may corresponding to a gantry coordinate not exactly (0,0). Is that OK?

[dlebauer]

Yes that is okay. Only the plots change. I think we should discuss in
person, but please map what you have so it is easier to understand.
On Fri, Sep 30, 2016 at 5:51 PM Ting Li notifications@github.com wrote:

@dlebauer https://github.com/dlebauer

From the information, I assume we are going to use one system and put all
the data with different bounding box into it. Also, the two planting
seasons at least must have to have the same (0,0) point (though we do not
know where it is ) and XY axes (I think this should be true). That means
the same gantry coordinate in both two seasons should exactly point to the
same point. Is this correct?
And, as the bounding box is changing and (0,0) point is unknown, the new
gantry coordinates with X increasing to the east will also not start
exactly at (0,0). The (0,0) point of the plots may corresponding to a
gantry coordinate not exactly (0,0). Is that OK?

—
You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub
#174 (comment),
or mute the thread
https://github.com/notifications/unsubscribe-auth/AAcX5zUnkX-4CWN05Er4XA2Z_G9y9ozLks5qvZJpgaJpZM4KFRZU
.

[dlebauer]

I was very confused and when you (@tingli3) pointed out x'=f(y) I realized that I misunderstood the gantry coordinate system.

But drawing a map helped me. This may be incorrect, but it is my current working hypothesis of how the different coordinate systems overlay.

If x goes positive to the north and y goes positive to the east, we don't need a coordinate system with x going east. That was only because if x increases to the west / to the left, images, plots, and perspective are inverted.

So, we should review these maps and be able to encode this in a database and use it to subset data and visualize the experimental design.

[p.s. Did anyone else ( @pless @LTBen @rjstrand @smarshall-bmr @doctorroboto et al ) catch on that x is increasing to the north?]

[tingli3]

This R code corrects the order of plotID and add column and range ID to each plot. Grid[i,] corresponds to plot i. The X1, X2, Y1 and Y2 need to be updated using the true gantry coordinates.

require(ggplot2)
require(udunits2)
x_origin = 0 
y_origin = 0
length = 3.5
alley = 0.5
row_spacing = 30
ranges = 54
columns = 16
row_spacing_cm <- ud.convert(row_spacing, 'in', 'm')
dx <- 2 * row_spacing_cm
dy <- length
x <- x_origin + (1:columns - 1) * dx
y <- y_origin + (1:ranges - 1) * (length + alley) 

grid <- data.frame(expand.grid(x,y))
names(grid)=c('x1','y1')
plot = 1:(ranges*columns)
grid=cbind(plot,grid)
grid$x2=grid$x1+dx
grid$y2=grid$y1+dy

grid$column=(grid$plot-1)%%16+1
grid$range=floor((grid$plot-1)/16)+1

for ( i in (1:27)*2 )
{
  grid[((i-1)*16+1):(i*16), 1]=c((i*16):((i-1)*16+1))
}

grid=grid[order(grid$plot),]

theme_set(theme_bw())
ggplot(data = grid) + 
  geom_text(aes(x1+dx/2, y1+dy/2, label = plot), size = 3.5) +
  geom_point(aes(x1,y1), size = 1, alpha = 0.4) + 
  geom_point(aes(x2,y2), col = 'red', size = 1, alpha = 0.4) +
  ggtitle('gantry plot layout') + xlab('x') + ylab('y')

[tingli3]

This is the R code for projection. Lat&lon can be returned through gantry2latlon(gantry_x, gantry_y) e.g. gantry2latlon(3.8,0). How the choose of the point (currently CE) and projection will effect the error needs to be experimented later.

require('proj4')
options(digits=15)
gantry2latlon <- function(gx,gy)
{
  #you can replace CE with other points such as SW SE ...
  CE=c(-111.974810166667,33.07547550)
  CE_gantry=c(101.750,0)
  dx_g=gx-CE_gantry[1]
  dy_g=gy-CE_gantry[2]
  proj= '+proj=tmerc +lat_0=0 +lon_0=-111.843479 +k=1 +x_0=12285.092664 +y_0=-3661028.344732 +datum=WGS84 +units=m +no_defs '
  #or you can use
  #proj= '+proj=tmerc +lat_0=0 +lon_0=-111.974928458333 +k=1 +x_0=50000 +y_0=0 +datum=WGS84 +units=m +no_defs '
  #proj='+proj=utm +zone=12 +ellps=GRS80 +datum=NAD83 +units=m +no_defs'

  CE_UTM=project(CE,proj)
  x_UTM=CE_UTM[1]-dy_g
  y_UTM=CE_UTM[2]+dx_g
  #return (round(project(c(x_UTM,y_UTM),proj,inverse=T),15))
  return (project(c(x_UTM,y_UTM),proj,inverse=T))
}

[yanliu-chn]

Next:

• generate polygons to be inserted to betydb in postgis sql statements. send to Scott.
• two site tables for two seasons. tables use USDA GPS reference points.
  https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=883764630

https://docs.google.com/spreadsheets/d/1eQSeVMPfrWS9Li4XlJf3qs2F8txmddbwZhjOfMGAvt8/edit#gid=1433617856

• work out sql query statements on querying the tables to output MAC and gantry coordinates.

Due Oct 14th.

[yanliu-chn]

From Andy French on Sept 2, 2016:

1. Use the USDA RTK system as the base reference, it is likely more accurate than the M&S (this is not CASE IH as understand). The USDA RTK system was checked a while ago (2007 so not so recently) three AZ base stations  (designations AZMO, AZSC, AZCLL) and differed from the NAD83 reference frame by 2, 1.2, 1.9 cm for the X,Y, Z coordinates respectively. We routinely check our locations against the USGS National Map and locations are typically within ~2cm
   

2.   We (Bill Luckett and I) checked the location of the M&S base station with the USDA rover (with some error because the M&S antenna is behind a locked fence to the west of the scanalyzer field). Thanks to Pedro for additional information we’ve confirmed the 2.5 m discrepancy is due to an inaccurate position used in commissioning the M&S base. This discrepancy appears to be constant across MAC fields. I suspect the M&S system precision is excellent even though its accuracy is not. Because of these factors, moving between the ‘best’ coordinates and the M&S coordinates is a matter of book-keeping to make sure one records your GPS source GPS.
   

3. Bill and Bob Strand have collected locations of the concrete piers using the USDA RTK system. Thus a linear transformation matrix can be easily created to move between scanalyzer coordinates and the currently best possible UTM12 coordinates. This matrix and the offset matrix from Existing software that can be leveraged or coupled to the pipeline #2 can be stored in your PostGIS system.

4.   At some unknown time in the future the inaccuracy of the M&S system will have to be addressed, but as this is not under UA control, cannot be forced. To me this represents a risk, maybe a small one, that one day we will be surprised to find the bias go away. My opinion is that it would be better to maintain more control over your geolocation positioning accuracy.
   

[yanliu-chn]

@tingli3 Upon Andy comment 3: need to confirm with him the EPSG number:
http://spatialreference.org/ref/epsg/26948/
http://spatialreference.org/ref/epsg/26949/
http://spatialreference.org/ref/epsg/26950/

[yanliu-chn]

@tingli3 Upon Andy comment 2: Need to clarify if the 2.5m is on x or y or both directions by going through the email. Ask if needed.

[yanliu-chn]

@tingli3 Upon Andy comment 2: need to ask Andy or Bob whether the Lematec coord derived using the linear regression results for MAC coord system has been corrected in the definition of b0, b1 and b2.

[yanliu-chn]

@tingli3 Need to confirm with Andy and Bob or @dlebauer whether the lat lon info associated with Lemnatec gantry raw data is the corrected values based on Andy's comment 1. This is related to the USDA point coord in column F-G row 4-7 in Coordinate Reference System Summary

[dlebauer]

I tweaked the code to insert polygons for season 2 by adding

  betydb_insert[[i]] <- paste("INSERT INTO sites (sitename, geometry, updated_at, created_at) VALUES (", 
                         paste0("'MAC Field Scanner Field Plot ", grid$plot[i], " Season 2', "),
                         geom_latlon, ", now(), now());")
}

to the final for loop.

and then filled out the information by running

update sites set city = 'Maricopa', state = 'Arizona', country = 'USA', time_zone = 'America/Phoenix', greenhouse = 'f', user_id = 6000000004, updated_at = now() where sitename like 'MAC Field%';

the geom_gantry information will be useful for extractors that use the gantry coordinate system, but it does not fit into BETYdb. Where can it be stored? Or can we modify the geometry column so that it can optionally return these points?

[gsrohde]

@dlebauer (and @tingli3) Please don't include the created_at and/or updated_at columns in the insert and update statements! A trigger function handles these automatically so that created_at and updated_at are set to the current time when a row is inserted, and updated_at is reset to the current time when a row is updated. Setting them explicitly overrides this behavior, and using now() sets the column to machine time, not UTC time. (There is a function I defined, utc_now(), which will return UTC time, but you don't need to use it in this case.)

[gsrohde]

@yanliu-chn I think the lines in the MAC --> USDA (M2U) section of https://github.com/terraref/documentation/blob/version1/spatial-reference-system-data.md are mis-labelled: Mx/Ux is the longitude and My/Uy is the latitude, correct? So is 0.000020308287 the longitude offset or the y-offset (which is really latitude)? Also, are you still assuming Mx/Ux are in UTM units here, or are these in degrees? @tingli3's script seems to assume they are degrees; he applies the offsets after doing the proj4 projection. If they really are degrees, could you clarify this in the documentation?

[yanliu-chn]

@gsrohde Thanks for pointing this out. I just fixed it. It was a typo. You are right. Mx/Ux is longitude in degrees; My/Uy is latitude in degrees. I have specified in the original doc that it was to convert from MAC lat,lon to USDA lat,lon.

[gsrohde]

@yanliu-chn Thanks for fixing. I guess using the terms "lat" and "lon" makes it clear that in the MAC --> USDA (M2U) section we are no longer working in UTM coordinates. What was confusing to me is that in the Gantry --> MAC section, the variable Mx is used to represent the easting coordinate in the UTM system whereas in the MAC --> USDA section, Mx is used to represent the longitude. (Similarly for My.)

[gsrohde]

@dlebauer, @yanliu-chn I'm not sure which SQL statements I should be looking at. The statements generated by the code in SQL.zip mentioned in comment #174 (comment) doesn't seem to directly relate to BETYdb, but @dlebauer references SQL statements designed for BETYdb in his comment #174 (comment).

[dlebauer]

@gsrohde I unassigned you when I realized more work needed to be done. But I went ahead and inserted these ....
insert_season2_sites.sql.txt

for this case the incorrect use of now() is trivial, but now that I know I won't specify this.

[tingli3]

@dlebauer Please give me an example of the sql statement I should generate.

[dlebauer]

@tingli3 please see #174 (comment)

[gsrohde]

@tingli3 ... but leave out the created_at and updated_at columns.

[tingli3]

@gsrohde Can you give me an example sql?

[gsrohde]

@tingli3 The file David attached to the comment #174 (comment) has several example insertion statements. Here's the first (but with the timestamp columns omitted, as they should be):

INSERT INTO sites (sitename, geometry) VALUES ( 'MAC Field Scanner Field Plot 1 Season 2',  ST_GeomFromText('POLYGON((-111.975049874375 33.0745312921391 353, -111.975033517034 33.0745313124814 353, -111.975033529346 33.0745670737771 353, -111.975049886694 33.0745670534346 353, -111.975049874375 33.0745312921391 353))', 4326));

[tingli3]

season1_sql.txt
season2_sql.txt
Could someone please help me to verify these results?

[dlebauer]

@tingli3

• I've already generated and inserted the season2 plots.
• The season 1 plots should have update rather than insert statements; the sitenames do not include Season 1` for the first season
• please add your code to the github repository here: https://github.com/terraref/computing-pipeline/blob/master/scripts/geospatial/field_scanner_plots.R

[dlebauer]

@tingli3 issue #187 will check for consistency with other field measurements.

[tingli3]

@dlebauer The season 1 and season 2 sql statements are based on the example given. If season 1 is different from season 2, please give me another example. Also, I have added the code.

[dlebauer]

@tingli3 I already inserted season 2

For season 1, the first line would be

UPDATE sites 
  SET geometry = 
    ST_GeomFromText('POLYGON((
      -111.975049873054 33.0745274579941 353, 
      -111.975033515714 33.0745274783365 353, 
      -111.975033525032 33.0745545428892 353, 
      -111.975049882377 33.0745545225468 353, 
      -111.975049873054 33.0745274579941 353))', 4326));
  WHERE sitename = 
       'MAC Field Scanner Field Plot 1';

P.S. I reformatted your R code following http://adv-r.had.co.nz/Style.html so ... please pull changes before updating.

[tingli3]

season1_sql.txt

[dlebauer]

@tingli3 all set - thanks! Closing - can you please update documentation (terraref/documentation#7)?