0

I have a bunch of points in a text file that are represented using lat long (116.3158516 39.9748016). I import them into my GIS enabled database and create trajectories using st_makeline() as geometries projected in EPSG:4326. I'm using geometry here because all the GPS data falls into an area that is 50km^2: correct me if this is my first mistake.

Now I'm doing things like calculating trajectory length (st_length) and simplifying the geometry using st_simplifypreservetopology() but according to the documentation the results for the length and the tolerance for the simplification are in the coordinate system of the data; which in my case is degrees. This obviously doesn't make sense. I need it to be in meters.

Do I need to change the way I'm storing my data or do I need to convert my trajectories and points into something else before feeding them into the st_length() and st_simplifypreservetopology() functions?

Screenshot for Mike T.

enter image description here

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Przemek Lach
  • 229
  • 3
  • 8

2 Answers2

1

You are correct that if your data are in a small region that you should use a geometry type. However, you would also normally transform the Lat/Long coordinates to a spatial reference system for the region. Normally most folks choose a UTM Zone, which describes coordinates for Eastings and Northing with length units in metres.

Check out ST_Transform to convert your Lat/Long coordinates to Easting/Northings. You can either modify the geometry column in-place, or create another table or geometry column.

Mike T
  • 42,095
  • 10
  • 126
  • 187
  • Thanks for the suggestion. I'm also having projection issues which may have something to do with the way I'm storing these values. When I stitch these points together to form trajectories (st_makeline) that I then add as a layer over Google Maps the trajectories are a bit off relative to where they should be on the map. My other concern is that that data source for these lat/long describes the data to be in WGS84 which if I understand correctly could mean anything. I assumed it to be 4326: the data was collected using smartphones. – Przemek Lach Oct 30 '14 at 04:15
  • Yes, WGS84 is SRID=4326. Location estimates from smartphones are not great, so it is expected that your lines would look like zig-zags jittered around the expected locations. – Mike T Oct 30 '14 at 04:53
  • Ya I expect them to bounce around but they seem to be off almost too perfectly. I've attached a screenshot to the original post. Have a look and let me know if it looks like a projection issue to you. The two trajectories are green and red. Thanks again. – Przemek Lach Oct 30 '14 at 05:18
  • Looks like there is some bias, but I'm not convinced it is a projection issue since some parts match up. Also, I'm no expert with how smartphones obtain their position, but check out this thread. – Mike T Oct 30 '14 at 05:41
  • when I created my points using ST_MakePoint() I left the elevation in it's original units: feet. I'm not certain how the elevation plays into drawing these lines. Maybe I should be storing elevation in meters? – Przemek Lach Oct 30 '14 at 19:15
  • If you transform your coordinates into easting and northing, it would also make sense to also convert elevation into metres, so that functions like ST_3DLength work. However, it doesn't influence projection issues on a flat map. – Mike T Oct 30 '14 at 21:05
  • Ok. When I run the query "SELECT st_astext(st_transform(ST_SetSRID(ST_MakePoint('116.3196983', '39.9848666', '100'), 4326), 4326));" I get the printout in degrees. What am I doing wrong with my transform to get easting/northing? – Przemek Lach Oct 30 '14 at 21:39
  • You need to transform to a different SRID. You didn't transform anything since 4326=4326. – Mike T Oct 30 '14 at 21:41
  • Right ok, I think I'm finally getting this. So since my data is in Beijing I'm going to transform it using WGS84/UTM Zone 50s which is EPSG:32750. Do I have that right? SELECT st_astext(st_transform(ST_SetSRID(ST_MakePoint('116.3196983', '39.9848666', '100'), 4326), 32750)); gives me POINT Z (441916.593167786 14426299.1478804 100) – Przemek Lach Oct 30 '14 at 22:11
  • Since Beijing is in the Northern hemisphere, you want zone 50N, or EPSG:32650, although 50S will sort-of work. – Mike T Oct 30 '14 at 22:28
  • I used this chart to determine the correct zone: http://upload.wikimedia.org/wikipedia/commons/e/ed/Utm-zones.jpg. Am I looking at the wrong chart? 50N would put me way south. – Przemek Lach Oct 30 '14 at 22:30
  • That chart is misleading. There is only S = south and N = north. China is in the North. – Mike T Oct 30 '14 at 22:44
  • Let us continue this discussion in chat. – Przemek Lach Oct 30 '14 at 23:08
0

When using st_length() I can convert to geography using st_length(line.line_geometry :: geography) as was suggested or using ST_Length2D_Spheroid(line.geometry, 'SPHEROID["GRS 1980",6378137,298.257222101]'). Both answers seem to be the same as far as precision is concerned: 580139.3016159851 vs 580139.3016215176. I'm guessing the slight difference is due to a rounding error. Using either of these two methods produces my results in meters which is what I was after.

I verified that the length was correct by using the Measure Line tool in QGIS and measuring a given line.

Przemek Lach
  • 229
  • 3
  • 8