ChHarding / TouchTerrain_for_CAGEO

Touch Terrain: A python app to create 3D printable terrain models (STL/OBJ) from only elevation data (via Google Earth Engine) or from a local geotiff. Has been used for CNC terrain models. Runs as a web app (http://touchterrain.org), as .py file (standalone.py) or as jupyter notebook. Docker image: https://github.com/ChHarding/TouchTerrain_jupyter_docker
http://touchterrain.geol.iastate.edu
191 stars 44 forks source link

Converting Lat/Lng to X/Y position on tile #58

Open lukezirngibl opened 1 year ago

lukezirngibl commented 1 year ago

I am trying to project a Lat/Lng pin onto the 3D printed tile. Given the topright & bottomleft coordinates.. How would you suggest mapping any given latlng to XY position on the print?

ChHarding commented 1 year ago

Luke,

To do this properly you’d have to project your pin location and the corner coords into UTM. The UTM zone (e.g. UTM 32N) and its EPSG number (32632) is in the logfile. After that you can just use proportions and you true mm size to get to your pin location in mm space.

This https://products.aspose.app/gis/transformation/convert-to-utm is useful for the conversion. Your input system is 4326. But, I’m not sure about the order on the output “lat/long” i.e. is it x/y or y/x. (I’m guessing y/x)

[cid:FC21A8D7-1B80-43CD-B2F8-C6F1F0FB69BB]

Is there a requirement to add your pin to the STL? Or could you simply put a physical pin on your model?

Cheers

Chris

On Aug 23, 2022, at 04:55, Luke Zirngibl @.**@.>> wrote:

I am trying to project a Lat/Lng pin onto the 3D printed tile. Given the topright & bottomleft coordinates.. How would you suggest mapping any given latlng to XY position on the print?

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5JNG2V34JLXNBVF7BLV2SN2VANCNFSM57KZUK3A. You are receiving this because you are subscribed to this thread.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

lukezirngibl commented 1 year ago

Hey Chris,

Thanks for the speedy response. The goal is to project a pin onto the tile (not integrate into the model). My current approach is as follows:

My touchterrain coordinates are:

const trlat = 47.851128701997794
const trlon = 10.846901928710908
const bllat = 45.79197038566249
const bllon = 5.514870678710881

My logfile does indeed show the EPSG:32632

I first converted both my TR and BL coordinates into UTM. Since they lie in different UTM Zones (31 & 32), I projected the BL to Zone 32. This yielded:

BL - UTM = [229133.28260624083, 5076843.084432984]
TR - UTM = [638164.1, 5301405.3]

My test target pin LatLng is:

const target = {
  lat: 47.366935,
  lng: 8.543069,
}

Converted to UTM is: [465500.1, 5246042.8]

And now i simply map that coordinate to my computed UTM coordinate space:

const convert = (latlng: LatLng) => { ... converts to correct UTM zone (32) }  

// map proportions are 16:9

const WIDTH = 1600
const HEIGHT = 900

const bottomLeftUTM = convert(bottomLeft)
const topRightUTM = convert(topRight)

const UTM_X_ZERO = bottomLeftUTM[0]
const UTM_Y_ZERO = bottomLeftUTM[1]

const UTM_WIDTH = topRightUTM[0] - bottomLeftUTM[0]
const UTM_HEIGHT = topRightUTM[1] - bottomLeftUTM[1]

const getXY = (latlng: LatLng): Coordinates => {
  const targetUTM = convert(latlng)
  return {
    x: ((targetUTM[0] - UTM_X_ZERO) / UTM_WIDTH) * WIDTH,
    y: ((targetUTM[1] - UTM_Y_ZERO) / UTM_HEIGHT) * HEIGHT,
  }
}

This gives me: {x: 924.5926999821821, y: 678.1182828367208} for the pin...

But its still about 1 cm off when projected onto the physical tile! Is simply converting UTM not enough? Is there some sort of distortion happening when you generate the STL? I'm not quite sure what you mean by "use proportions and your true mm size".

FYI - The map is 1.2m wide (6 tiles wide, 4 tiles high). 16:9 proportions.

Thank you for your help!

lukezirngibl commented 1 year ago

To give you a sense of it, here is a screenshot (right now i'm just doing the projection in pixel space).

Red = actual Black = expected

Screen Shot 2022-08-26 at 09 15 18
ChHarding commented 1 year ago

Luke,

I re-did you math using ArcGIS and I can confirm x=924.59 y=678.11 mathematically correct. ( Also, the height/width ratios don’t match, 900mm / 1600mm is 0.5625 whereas the “box” in utm is 224562m / 409030m = 0.549011 If you look in the log file I do make some minor adjustments but I doubt it’d lead to that much distortion. Can you send my that log file? Or maybe we should just use the ratio between these (0.976019) as correction factor? That would make it x=902.4 and y=661.8. How would that look on your tile?

Cheers

Chris

On Aug 26, 2022, at 02:10, Luke Zirngibl @.**@.>> wrote:

Hey Chris,

Thanks for the speedy response. The goal is to project a pin onto the tile (not integrate into the model). My current approach is as follows:

My coordinates are:

const trlat = 47.851128701997794 const trlon = 10.846901928710908 const bllat = 45.79197038566249 const bllon = 5.514870678710881

My logfile does indeed show the EPSG:32632

I first converted both my TR and BL coordinates into UTM. Since they lie in different UTM Zones (31 & 32), I projected the BL to Zone 32. This yielded:

BL - UTM = [229133.28260624083, 5076843.084432984] TR - UTM = [638164.1, 5301405.3]

My test target ping LatLng is:

const target = { lat: 47.366935, lng: 8.543069, }

Converted to UTM is: [465500.1, 5246042.8]

Any now i simply do

// map proportions are 16:9

const WIDTH = 1600 const HEIGHT = 900

const bottomLeftUTM = convert(bottomLeft) const topRightUTM = convert(topRight)

const UTM_X_ZERO = bottomLeftUTM[0] const UTM_Y_ZERO = bottomLeftUTM[1]

const UTM_WIDTH = topRightUTM[0] - bottomLeftUTM[0] const UTM_HEIGHT = topRightUTM[1] - bottomLeftUTM[1]

const getXY = (latlng: LatLng): Coordinates => { const targetUTM = convert(latlng) console.log(targetUTM) return { x: ((targetUTM[0] - UTM_X_ZERO) / UTM_WIDTH) WIDTH, y: ((targetUTM[1] - UTM_Y_ZERO) / UTM_HEIGHT) HEIGHT, } }

This gives me: {x: 924.5926999821821, y: 678.1182828367208} for the pin...

But its still about 1 cm off when projected onto the physical tile! Is simply converting UTM not enough? Is there some sort of distortion happening when you generate the STL? I'm not quite sure what you mean by "use proportions and you true mm".

FYI - The map is 1.2m wide (6 tiles wide, 4 tiles high). 16:9 proportions.

Thank you for your help!

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1228136022, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5MEL5WQGPTB7QRJK4LV3BUVZANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

lukezirngibl commented 1 year ago

Just using the "correction factor" doesn't seem to do it since the distortion doesn't appear to be linear. For the single point it's OK, but for others, its equally inaccurate.

Here is my log file:

Log for creating 6 x 4 3D model tile(s) from switzerland-16-9.tif 

started: 09:43:52.647827 
source raster upper left corner (x/y):  229110.0 5305650.0 
source raster cells size 30.0 m  (7769, 13815) 
projection: WGS 84 / UTM zone 32N 
z-scale: 1.3 
basethickness: 15 
fileformat: STLb 
tile_centered: False 
no_bottom: False 
no_normals: True 
ignore_leq: None 
lower_leq: None 
importedGPX: None 
Warning: Given outline polygon will be ignored when using local raster file! 
undefined GDAL value: None 
tile_width: 200 
tile_height: 112.4719507781397 
source raster width 414450.0 m, cell size: 30.0 m, elev. min/max is 40.0 4809.0 m 
source raster 3D print resolution would be 0.08686210640608034 mm 
re-sampling switzerland-16-9.tif :
  (13815, 7769) 0.08686210640608034 mm  30.0 m  40.0 - 4809.0 m to 
  (5999, 3374) 0.20003333888981498 mm  69.075 m  51.31341 - 4796.6377 m 
after resampling, requested print res was adjusted from 0.2 to 0.20003333888981498 to ensure correct model dimensions 
Cropping for nice fit of 6 (width) x 4 (height) tiles, removing: 5 columns, 2 rows 
 cropped (5999, 3374) to (5994, 3372) 
 cropping changed physical size from 200 mm x 168.72812135355892 mm 
 to 199.83330555092516 mm x 168.62810468411402 mm 
map scale is 1 : 345317.4375 
Cells per tile (x/y) 999 x 843 
using single-core only 

I suppose the real question is: Why does the computed UTM coordinate space not have the same proportions as the tiles?

Could it be that the cropping is changing the actual TR / BL? Therefore shifting the entire system by a few millimeters. And could subdividing the print (into 24 tiles) somehow compound the issue?

ChHarding commented 1 year ago

Well, I’ll need to run a test to get to the bottom of this.

On Aug 27, 2022, at 03:18, Luke Zirngibl @.**@.>> wrote:

Just using the "correction factor" doesn't seem to do it since the distortion doesn't appear to be linear. For the single point it's OK, but for others, it equally inaccurate.

Here is my log file:

Log for creating 6 x 4 3D model tile(s) from switzerland-16-9.tif

started: 09:43:52.647827 source raster upper left corner (x/y): 229110.0 5305650.0 source raster cells size 30.0 m (7769, 13815) projection: WGS 84 / UTM zone 32N z-scale: 1.3 basethickness: 15 fileformat: STLb tile_centered: False no_bottom: False no_normals: True ignore_leq: None lower_leq: None importedGPX: None Warning: Given outline polygon will be ignored when using local raster file! undefined GDAL value: None tile_width: 200 tile_height: 112.4719507781397 source raster width 414450.0 m, cell size: 30.0 m, elev. min/max is 40.0 4809.0 m source raster 3D print resolution would be 0.08686210640608034 mm re-sampling switzerland-16-9.tif : (13815, 7769) 0.08686210640608034 mm 30.0 m 40.0 - 4809.0 m to (5999, 3374) 0.20003333888981498 mm 69.075 m 51.31341 - 4796.6377 m after resampling, requested print res was adjusted from 0.2 to 0.20003333888981498 to ensure correct model dimensions Cropping for nice fit of 6 (width) x 4 (height) tiles, removing: 5 columns, 2 rows cropped (5999, 3374) to (5994, 3372) cropping changed physical size from 200 mm x 168.72812135355892 mm to 199.83330555092516 mm x 168.62810468411402 mm map scale is 1 : 345317.4375 Cells per tile (x/y) 999 x 843 using single-core only

I suppose the real question is: Why does the computed UTM coordinate space not have the same proportions as the tiles?

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1229148911, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5KQKZYSPDELYXEYTK3V3HFLTANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

ChHarding commented 1 year ago

Actually, one more question: How did you get this geotiff? Did you download it via the EarthEngine using the JS export script?

On Aug 27, 2022, at 03:18, Luke Zirngibl @.**@.>> wrote:

Just using the "correction factor" doesn't seem to do it since the distortion doesn't appear to be linear. For the single point it's OK, but for others, it equally inaccurate.

Here is my log file:

Log for creating 6 x 4 3D model tile(s) from switzerland-16-9.tif

started: 09:43:52.647827 source raster upper left corner (x/y): 229110.0 5305650.0 source raster cells size 30.0 m (7769, 13815) projection: WGS 84 / UTM zone 32N z-scale: 1.3 basethickness: 15 fileformat: STLb tile_centered: False no_bottom: False no_normals: True ignore_leq: None lower_leq: None importedGPX: None Warning: Given outline polygon will be ignored when using local raster file! undefined GDAL value: None tile_width: 200 tile_height: 112.4719507781397 source raster width 414450.0 m, cell size: 30.0 m, elev. min/max is 40.0 4809.0 m source raster 3D print resolution would be 0.08686210640608034 mm re-sampling switzerland-16-9.tif : (13815, 7769) 0.08686210640608034 mm 30.0 m 40.0 - 4809.0 m to (5999, 3374) 0.20003333888981498 mm 69.075 m 51.31341 - 4796.6377 m after resampling, requested print res was adjusted from 0.2 to 0.20003333888981498 to ensure correct model dimensions Cropping for nice fit of 6 (width) x 4 (height) tiles, removing: 5 columns, 2 rows cropped (5999, 3374) to (5994, 3372) cropping changed physical size from 200 mm x 168.72812135355892 mm to 199.83330555092516 mm x 168.62810468411402 mm map scale is 1 : 345317.4375 Cells per tile (x/y) 999 x 843 using single-core only

I suppose the real question is: Why does the computed UTM coordinate space not have the same proportions as the tiles?

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1229148911, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5KQKZYSPDELYXEYTK3V3HFLTANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

ChHarding commented 1 year ago

Thinking about it some more, why don’t you convert your plate locations to UTM 32N first? As you have the x/y of the upper left corner you could get the offset of your location to that corner and, knowing that each cell is 30 m, you could then figure out at which cell your location is and lower the cell (and maybe its surrounding cell) by some amount, so it would make a hole when printed. Or how else where you planning to “mark” the location in your print?

On Aug 30, 2022, at 12:37, Harding, Chris [GE AT] @.**@.>> wrote:

Actually, one more question: How did you get this geotiff? Did you download it via the EarthEngine using the JS export script?

On Aug 27, 2022, at 03:18, Luke Zirngibl @.**@.>> wrote:

Just using the "correction factor" doesn't seem to do it since the distortion doesn't appear to be linear. For the single point it's OK, but for others, it equally inaccurate.

Here is my log file:

Log for creating 6 x 4 3D model tile(s) from switzerland-16-9.tif

started: 09:43:52.647827 source raster upper left corner (x/y): 229110.0 5305650.0 source raster cells size 30.0 m (7769, 13815) projection: WGS 84 / UTM zone 32N z-scale: 1.3 basethickness: 15 fileformat: STLb tile_centered: False no_bottom: False no_normals: True ignore_leq: None lower_leq: None importedGPX: None Warning: Given outline polygon will be ignored when using local raster file! undefined GDAL value: None tile_width: 200 tile_height: 112.4719507781397 source raster width 414450.0 m, cell size: 30.0 m, elev. min/max is 40.0 4809.0 m source raster 3D print resolution would be 0.08686210640608034 mm re-sampling switzerland-16-9.tif : (13815, 7769) 0.08686210640608034 mm 30.0 m 40.0 - 4809.0 m to (5999, 3374) 0.20003333888981498 mm 69.075 m 51.31341 - 4796.6377 m after resampling, requested print res was adjusted from 0.2 to 0.20003333888981498 to ensure correct model dimensions Cropping for nice fit of 6 (width) x 4 (height) tiles, removing: 5 columns, 2 rows cropped (5999, 3374) to (5994, 3372) cropping changed physical size from 200 mm x 168.72812135355892 mm to 199.83330555092516 mm x 168.62810468411402 mm map scale is 1 : 345317.4375 Cells per tile (x/y) 999 x 843 using single-core only

I suppose the real question is: Why does the computed UTM coordinate space not have the same proportions as the tiles?

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1229148911, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5KQKZYSPDELYXEYTK3V3HFLTANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://touchterrain.geol.iastate.edu/

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

lukezirngibl commented 1 year ago

We used your approach from: https://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/34 to generate the geotiff (yeah, from Earth Engine)

OK - so i like this idea of converting to your cell coordinates. How do I convert from UTM coord to cell coord? Or how do you go from:

source raster upper left corner (x/y): 229110.0 5305650.0

to

source raster cells size 30.0 m (7769, 13815)

End goal is to project a dot onto the model. So I think LatLng -> Cell position would make most sense

ChHarding commented 1 year ago

Lets say you know that your raster’s upper left cell is as x=1000, y=1000 in UTM and you cell size is 100 m. In other words in raster-land the origin is always in the upper left and x goes to the right and y goes down.

Say, in your raster you have 100 cells in x and 100 cells in y, so your x ranges from 1000 to 2000 (west to east) and your y starts at 1000 and goes down to 0 (north to south)

Your lat/lon location converts to. x=1505, y=807. So its x offset (to the right) from 1000 is 1505 - 1000 = 505 and it y offset (down) from 1000 is 1000 - 807 = 193. Raster index numbers (cell coordinates) usually is done as I along x (here from 0 to 99) and j along y (also 0 to 99) Or sometimes people call along x “columns" and along y “rows". Converted to raster cell index wise this would be i/column = 505 / 10 = 50.5 (so actually its sitting on 51) and j = 193 / 10 - 19.3 = 20

So you location would be sitting on cell 51/20

On Aug 31, 2022, at 10:08, Luke Zirngibl @.**@.>> wrote:

We used your approach from: #34https://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/34 to generate the geotiff (yeah, from Earth Engine)

OK - so i like this idea of converting to your cell coordinates. How do I convert from UTM coord to cell coord? Or how do you go from:

source raster upper left corner (x/y): 229110.0 5305650.0

to

source raster cells size 30.0 m (7769, 13815)

End goal is to project a dot onto the model. So I think LatLng -> Cell position would make most sense

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1233065189, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5NB2F27LFPRRZHQQDLV35YP7ANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu

ChHarding commented 1 year ago

Once we have a cell coordinate, this makes perfect sense. But it is unclear to me how we translate from LatLng to Cell [x,y].

So far I have:

Step 1: LatLng -> UTM.

const trlat = 47.851128701997794 const trlon = 10.846901928710908 const bllat = 45.79197038566249 const bllon = 5.514870678710881

This gives us a top left coordinate of [229,133, 5,301,405.3] in 32N. Which is roughly equal to yours (in the logfile)

When I use the online converter I get. x=239,293 and y=5,305,636:

[cid:C523E588-137B-494B-8846-81720E922386]

How do you get your [229133, 5301405.3]? Am I overlooking something?

source raster upper left corner (x/y): 229110.0 5305650.0

Why it's slightly different, i don't know..

Step 2: generate cell coordinate system?

source raster cells size 30.0 m (7769, 13815)

What is this?

tile_width: 200 tile_height: 112.4719507781397 source raster width 414450.0 m, cell size: 30.0 m, elev. min/max is 40.0 4809.0 m

^ I think whatever is happening here is the source of our distortion and my own misunderstanding

— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/58#issuecomment-1233903553, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5P5GAZQ4PE3M6JSMNDV4BPX3ANCNFSM57KZUK3A. You are receiving this because you commented.Message ID: @.***>

Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu