Closed evanthomas closed 1 year ago
How big is nationalz55_ag.tif ? If’s 641646 * 137974 pixels and each pixel is 4 bytes (i.e. your elevation is encoded as floating point data) it would be 354 Gb uncompressed, and you’re running out of memory when that geotiff is read in.
On Mar 3, 2023, at 02:25, evanthomas @.**@.>> wrote:
I'm trying to run TT in the docker notebook against an Aussie DEM. The DEM was sourced from Geosciences Australia and seems to load into the map without issues.
Here is a code that reproduces the problem:
from geojson import Polygon poly = Polygon([[[145.147961, -37.874657], [145.146262, -37.874848], [145.144699, -37.875406], [145.143395, -37.876287], [145.142454, -37.87742], [145.141952, -37.878715], [145.141928, -37.880069], [145.142384, -37.881375], [145.143285, -37.882528], [145.144557, -37.883437], [145.146101, -37.884029], [145.147792, -37.884257], [145.149497, -37.884103], [145.15108, -37.88358], [145.152415, -37.882728], [145.153394, -37.881616], [145.153942, -37.880332], [145.154013, -37.878979], [145.153602, -37.877664], [145.152743, -37.876491], [145.151503, -37.875555], [145.149981, -37.87493], [145.148298, -37.874665], [145.147961, -37.874657]]]) args = {'importedDEM': '/data/nationalz55_ag.tif', 'DEM_name': 'JAXA/ALOS/AW3D30/V2_2', 'bllat': -37.88327658196921, 'bllon': 145.13953117216852, 'trlat': -37.87247077823912, 'trlon': 145.15559581602838, 'tilewidth': 200, 'printres': 0.4, 'ntilesx': 1, 'ntilesy': 1, 'basethick': 1, 'zscale': 3, 'fileformat': 'STLb', 'zip_file_name': 'myterrain', 'tile_centered': False, 'CPU_cores_to_use': 0, 'max_cells_for_memory_only': 5002, 'no_bottom': False, 'no_normal': True, 'bottom_image': None, 'ignore_leq': None, 'lower_leq': None, 'unprojected': False, 'projection': None, 'only': None, 'importedGPX': [], 'polygon': poly} from touchterrain.common import TouchTerrainEarthEngine as TouchTerrain
totalsize, full_zip_file_name = TouchTerrain.get_zipped_tiles(**args) # args are in a dict print("\nDONE!\n\nCreated zip file", full_zip_file_name, "%.2f" % totalsize, "Mb")
This produces:
INFO:root:EE init() worked with .config/earthengine/credentials INFO:root:Using GeoJSON polygon for masking with 24 points INFO:root:Log for creating 1 x 1 3D model tile(s) from nationalz55_ag.tif
INFO:root:started: 08:18:47.116651 Log for creating 1 x 1 3D model tile(s) from nationalz55_ag.tif
MemoryError Traceback (most recent call last) Cell In[1], line 6 3 args = {'importedDEM': '/data/nationalz55_ag.tif', 'DEM_name': 'JAXA/ALOS/AW3D30/V2_2', 'bllat': -37.88327658196921, 'bllon': 145.13953117216852, 'trlat': -37.87247077823912, 'trlon': 145.15559581602838, 'tilewidth': 200, 'printres': 0.4, 'ntilesx': 1, 'ntilesy': 1, 'basethick': 1, 'zscale': 3, 'fileformat': 'STLb', 'zip_file_name': 'myterrain', 'tile_centered': False, 'CPU_cores_to_use': 0, 'max_cells_for_memory_only': 5002, 'no_bottom': False, 'no_normal': True, 'bottom_image': None, 'ignore_leq': None, 'lower_leq': None, 'unprojected': False, 'projection': None, 'only': None, 'importedGPX': [], 'polygon': poly} 4 from touchterrain.common import TouchTerrainEarthEngine as TouchTerrain ----> 6 totalsize, full_zip_file_name = TouchTerrain.get_zipped_tiles(**args) # args are in a dict 7 print("\nDONE!\n\nCreated zip file", full_zip_file_name, "%.2f" % totalsize, "Mb")
File /TouchTerrain/standalone/touchterrain/common/TouchTerrainEarthEngine.py:973, in get_zipped_tiles(DEM_name, trlat, trlon, bllat, bllon, polygon, polyURL, poly_file, importedDEM, printres, ntilesx, ntilesy, tilewidth, basethick, zscale, fileformat, tile_centered, CPU_cores_to_use, max_cells_for_memory_only, temp_folder, zip_file_name, no_bottom, bottom_image, ignore_leq, lower_leq, unprojected, only, original_query_string, no_normals, projection, use_geo_coords, importedGPX, gpxPathHeight, gpxPixelsBetweenPoints, gpxPathThickness, map_img_filename, smooth_borders, offset_masks_lower, fill_holes, min_elev, **otherargs) 971 dem = gdal.Open(importedDEM) 972 band = dem.GetRasterBand(1) --> 973 npim = band.ReadAsArray().astype(numpy.float64) 975 # Read in offset mask file 976 # offset_npim = [] # Define offset_npim in parent scope before importedDEM "if statement" to avoid python unbound error 977 if offset_masks_lower is not None:
File /usr/local/lib/python3.10/site-packages/osgeo/gdal.py:4142, in Band.ReadAsArray(self, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_type, buf_obj, resample_alg, callback, callback_data) 4137 """ Reading a chunk of a GDAL band into a numpy array. The optional (buf_xsize,buf_ysize,buf_type) 4138 parameters should generally not be specified if buf_obj is specified. The array is returned""" 4140 from osgeo import gdal_array -> 4142 return gdal_array.BandReadAsArray(self, xoff, yoff, 4143 win_xsize, win_ysize, 4144 buf_xsize, buf_ysize, buf_type, buf_obj, 4145 resample_alg=resample_alg, 4146 callback=callback, 4147 callback_data=callback_data)
File /usr/local/lib/python3.10/site-packages/osgeo/gdal_array.py:428, in BandReadAsArray(band, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_type, buf_obj, resample_alg, callback, callback_data) 426 if buf_type == gdalconst.GDT_Byte and band.GetMetadataItem('PIXELTYPE', 'IMAGE_STRUCTURE') == 'SIGNEDBYTE': 427 typecode = numpy.int8 --> 428 buf_obj = numpy.empty([buf_ysize, buf_xsize], dtype=typecode) 430 else: 431 if len(buf_obj.shape) not in (2, 3):
MemoryError: Unable to allocate 330. GiB for an array with shape (641646, 137974) and data type float32
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5LPQSTPZFHUZZF3SGDW2GTF3ANCNFSM6AAAAAAVOJ55I4. 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
I should also say that my work around for this would be to externally resample this DEM using a GIS (QGIS is free) to larger cellsize (i.e. lower number of cells) that is closer to your nozzle size. Ex.: if your 3D print is 400 mm wide and you nozzle is 0.4 mm, I would resample it to about 4000 pixels.
On Mar 3, 2023, at 11:41, Harding, Chris [GE AT] @.**@.>> wrote:
How big is nationalz55_ag.tif ? If’s 641646 * 137974 pixels and each pixel is 4 bytes (i.e. your elevation is encoded as floating point data) it would be 354 Gb uncompressed, and you’re running out of memory when that geotiff is read in.
On Mar 3, 2023, at 02:25, evanthomas @.**@.>> wrote:
I'm trying to run TT in the docker notebook against an Aussie DEM. The DEM was sourced from Geosciences Australia and seems to load into the map without issues.
Here is a code that reproduces the problem:
from geojson import Polygon poly = Polygon([[[145.147961, -37.874657], [145.146262, -37.874848], [145.144699, -37.875406], [145.143395, -37.876287], [145.142454, -37.87742], [145.141952, -37.878715], [145.141928, -37.880069], [145.142384, -37.881375], [145.143285, -37.882528], [145.144557, -37.883437], [145.146101, -37.884029], [145.147792, -37.884257], [145.149497, -37.884103], [145.15108, -37.88358], [145.152415, -37.882728], [145.153394, -37.881616], [145.153942, -37.880332], [145.154013, -37.878979], [145.153602, -37.877664], [145.152743, -37.876491], [145.151503, -37.875555], [145.149981, -37.87493], [145.148298, -37.874665], [145.147961, -37.874657]]]) args = {'importedDEM': '/data/nationalz55_ag.tif', 'DEM_name': 'JAXA/ALOS/AW3D30/V2_2', 'bllat': -37.88327658196921, 'bllon': 145.13953117216852, 'trlat': -37.87247077823912, 'trlon': 145.15559581602838, 'tilewidth': 200, 'printres': 0.4, 'ntilesx': 1, 'ntilesy': 1, 'basethick': 1, 'zscale': 3, 'fileformat': 'STLb', 'zip_file_name': 'myterrain', 'tile_centered': False, 'CPU_cores_to_use': 0, 'max_cells_for_memory_only': 5002, 'no_bottom': False, 'no_normal': True, 'bottom_image': None, 'ignore_leq': None, 'lower_leq': None, 'unprojected': False, 'projection': None, 'only': None, 'importedGPX': [], 'polygon': poly} from touchterrain.common import TouchTerrainEarthEngine as TouchTerrain
totalsize, full_zip_file_name = TouchTerrain.get_zipped_tiles(**args) # args are in a dict print("\nDONE!\n\nCreated zip file", full_zip_file_name, "%.2f" % totalsize, "Mb")
This produces:
INFO:root:EE init() worked with .config/earthengine/credentials INFO:root:Using GeoJSON polygon for masking with 24 points INFO:root:Log for creating 1 x 1 3D model tile(s) from nationalz55_ag.tif
INFO:root:started: 08:18:47.116651 Log for creating 1 x 1 3D model tile(s) from nationalz55_ag.tif
MemoryError Traceback (most recent call last) Cell In[1], line 6 3 args = {'importedDEM': '/data/nationalz55_ag.tif', 'DEM_name': 'JAXA/ALOS/AW3D30/V2_2', 'bllat': -37.88327658196921, 'bllon': 145.13953117216852, 'trlat': -37.87247077823912, 'trlon': 145.15559581602838, 'tilewidth': 200, 'printres': 0.4, 'ntilesx': 1, 'ntilesy': 1, 'basethick': 1, 'zscale': 3, 'fileformat': 'STLb', 'zip_file_name': 'myterrain', 'tile_centered': False, 'CPU_cores_to_use': 0, 'max_cells_for_memory_only': 5002, 'no_bottom': False, 'no_normal': True, 'bottom_image': None, 'ignore_leq': None, 'lower_leq': None, 'unprojected': False, 'projection': None, 'only': None, 'importedGPX': [], 'polygon': poly} 4 from touchterrain.common import TouchTerrainEarthEngine as TouchTerrain ----> 6 totalsize, full_zip_file_name = TouchTerrain.get_zipped_tiles(**args) # args are in a dict 7 print("\nDONE!\n\nCreated zip file", full_zip_file_name, "%.2f" % totalsize, "Mb")
File /TouchTerrain/standalone/touchterrain/common/TouchTerrainEarthEngine.py:973, in get_zipped_tiles(DEM_name, trlat, trlon, bllat, bllon, polygon, polyURL, poly_file, importedDEM, printres, ntilesx, ntilesy, tilewidth, basethick, zscale, fileformat, tile_centered, CPU_cores_to_use, max_cells_for_memory_only, temp_folder, zip_file_name, no_bottom, bottom_image, ignore_leq, lower_leq, unprojected, only, original_query_string, no_normals, projection, use_geo_coords, importedGPX, gpxPathHeight, gpxPixelsBetweenPoints, gpxPathThickness, map_img_filename, smooth_borders, offset_masks_lower, fill_holes, min_elev, **otherargs) 971 dem = gdal.Open(importedDEM) 972 band = dem.GetRasterBand(1) --> 973 npim = band.ReadAsArray().astype(numpy.float64) 975 # Read in offset mask file 976 # offset_npim = [] # Define offset_npim in parent scope before importedDEM "if statement" to avoid python unbound error 977 if offset_masks_lower is not None:
File /usr/local/lib/python3.10/site-packages/osgeo/gdal.py:4142, in Band.ReadAsArray(self, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_type, buf_obj, resample_alg, callback, callback_data) 4137 """ Reading a chunk of a GDAL band into a numpy array. The optional (buf_xsize,buf_ysize,buf_type) 4138 parameters should generally not be specified if buf_obj is specified. The array is returned""" 4140 from osgeo import gdal_array -> 4142 return gdal_array.BandReadAsArray(self, xoff, yoff, 4143 win_xsize, win_ysize, 4144 buf_xsize, buf_ysize, buf_type, buf_obj, 4145 resample_alg=resample_alg, 4146 callback=callback, 4147 callback_data=callback_data)
File /usr/local/lib/python3.10/site-packages/osgeo/gdal_array.py:428, in BandReadAsArray(band, xoff, yoff, win_xsize, win_ysize, buf_xsize, buf_ysize, buf_type, buf_obj, resample_alg, callback, callback_data) 426 if buf_type == gdalconst.GDT_Byte and band.GetMetadataItem('PIXELTYPE', 'IMAGE_STRUCTURE') == 'SIGNEDBYTE': 427 typecode = numpy.int8 --> 428 buf_obj = numpy.empty([buf_ysize, buf_xsize], dtype=typecode) 430 else: 431 if len(buf_obj.shape) not in (2, 3):
MemoryError: Unable to allocate 330. GiB for an array with shape (641646, 137974) and data type float32
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5LPQSTPZFHUZZF3SGDW2GTF3ANCNFSM6AAAAAAVOJ55I4. 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/
Chris Harding Associate Professor Department of Geological & Atmospheric Sciences Touchterrain.geol.iastate.eduhttp://Touchterrain.geol.iastate.edu
Yes, the file is around 350GB. There are not many machines that have that much RAM so I assumed that the software would just read the data for the polygon from the geotif. I'm not a geospatial expert but presumbaly I could use QGIS to extract the bit I'm interested and save that in a smaller file with losing resolution?
Yes, ideally I would not read in anything outside of the polygon but that’s not implemented. IN fact I’m pretty sure that using a polygon with a geotiff won’t actually clip it b/c I’d have to do this via GDAL (a common raster library) myself. It currently only works for online data b/c I can tell Google Earth Engine to clip it for me. I could look at implementing this but even then your raster would be far too large. But, I can show you how to clip your raster with a kms file (and even downsample it) in QGIS. Why don’t you send my the kms file and I’ll work up an example workflow?
On Mar 3, 2023, at 17:32, evanthomas @.**@.>> wrote:
Yes, the file is around 350GB. There are not many machines that have that much RAM so I assumed that the software would just read the data for the polygon from the geotif. I'm not a geospatial expert but presumbaly I could use QGIS to extract the bit I'm interested and save that in a smaller file with losing resolution?
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1454263546, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5NFWGOB3O4IK3E4TRTW2J5QPANCNFSM6AAAAAAVOJ55I4. 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
That's very generous of you but I'm only this for a bit of fun, I don't want to waste your time unless it was going to be more generally useful. I'm happy to bumble around in qgis (also I probably have access to some large memory nodes at work ;).
No worries! Here’s a folder with example data and instructions:
On Mar 3, 2023, at 21:02, evanthomas @.**@.>> wrote:
That's very generous of you but I'm only this for a bit of fun, I don't want to waste your time unless it was going to be more generally useful. I'm happy to bumble around in qgis (also I probably have access to some large memory nodes at work ;).
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1454351265, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5NKVOPHCW2KSRWTQHLW2KWDZANCNFSM6AAAAAAVOJ55I4. 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
I think the folder is missing.
Maybe I can’t attach a zip folder with GitHub. Here are the single files:
On Mar 5, 2023, at 14:43, evanthomas @.**@.>> wrote:
I think the folder is missing.
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1455200154, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5KEQ4KMK55EMQ7S753W2T3HTANCNFSM6AAAAAAVOJ55I4. 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
Sorry for being obtuse - I can't any files attached or associated with this issue. It looks like your updating this issue via email. Maybe attachments are being stripped?
No need to apologize! Yes, I’m using email and I guess I never realized this! I’ve put the zipped folder on my google drive here:
https://drive.google.com/file/d/1BfMz4Z8lpKvpn9aorZIbJ-4PlAzqolOW/view?usp=share_link.
Cheers
Chris
On Mar 6, 2023, at 17:47, evanthomas @.**@.>> wrote:
Sorry for being obtuse - I can't any files attached or associated with this issue. It looks like your updating this issue via email. Maybe attachments are being stripped?
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1457226768, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5IXKSXA5XEA6PQTNJTW2ZZPJANCNFSM6AAAAAAVOJ55I4. 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
Ok, that worked perfectly and save me hours of bumbling. I actually used geojson instead of kml since that was easy to export from the notebook and qgis ingested it without issue.
Thanks!
Cool! Did you actually load your huge file into QGIS or did you just run the tool w/o first viewing it?
On Mar 8, 2023, at 02:25, evanthomas @.**@.>> wrote:
Ok, that worked perfectly and save me hours of bumbling. I actually used geojson instead of kml since that was easy to export from the notebook and qgis ingested it without issue.
Thanks!
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1459716107, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5OTG5L44WKYU3JT3S3W3A7BHANCNFSM6AAAAAAVOJ55I4. 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
QGIS was able to open the file without problems. I thinking it was just loading the data for the view port. Scrolling and zooming was slow, even on a reasonably spec’ed computer but good enough for this purpose. The actual extraction was fast. Not sure if the gdal executable is available, or could be included, in the docker but you might be able to use that in the notebook.
I just dug a bit into the Python bindings for GDAL warp and was able to re-create what QGIS produces. So I will put that in to the standalone code base eventually to clip the geotiff with a kml and also Geojson. So thanks for making me investigate it :) Is what is the terrain you’ve clipped out? And what is it for?
On Mar 8, 2023, at 13:58, evanthomas @.**@.>> wrote:
QGIS was able to open the file without problems. I thinking it was just loading the data for the view port. Scrolling and zooming was slow, even on a reasonably spec’ed computer but good enough for this purpose. The actual extraction was fast. Not sure if the gdal executable is available, or could be included, in the docker but you might be able to use that in the notebook.
— Reply to this email directly, view it on GitHubhttps://github.com/ChHarding/TouchTerrain_for_CAGEO/issues/70#issuecomment-1460771294, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEYDF5NVHD6J3S5OBUMKBILW3DQG3ANCNFSM6AAAAAAVOJ55I4. 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
I just wanted a model of my neighbourhood that I could print. I thought sourcing data from Geosciences Australia would give me better resolution but it didn't. The standalone version did let create a create a circular model though!
closed
I'm trying to run TT in the docker notebook against an Aussie DEM. The DEM was sourced from Geosciences Australia and seems to load into the map without issues.
Here is a code that reproduces the problem:
This produces: