import rhinoscriptsyntax as rs import scriptcontext as sc from time import time import Rhino from Rhino.ApplicationSettings.AppearanceSettings import DefaultLayerColor """Imports a standard ASCII Grid file and creates a point cloud, mesh and/or surface Script by Mitch Heynick 19 November 2011 Updated 03.12.14 - added check for Mac, if yes turn off Status Bar progress Updated 25.02.16 - non-meter unit warning, default=mesh, save output preferences, zoom object at end,set Perspective to Rendered if output is surface or mesh Updated 08-11.17, added try/except for Z values in case of read failure, plus whitespace strip for lines that might have accidental spaces at beginning/end.""" def PC(): return Rhino.Runtime.HostUtils.RunningOnWindows def SurfaceType(rows,cols): #asks user for input; returns degree plus True for Interpolate, False for Control str1 = "InterpolatePoints" ; str2 = "ControlPoints" srfType=rs.GetBoolean("Surface type",[("Type",str2,str1)],(True)) if srfType==None: return srfDeg = rs.GetInteger("Surface Degree?", 3, 1, 5) if srfDeg==None: return if srfDeg > min(rows,cols)-1 : srfDeg = min(rows,cols)-1 return (srfType[0],srfDeg) def CreateFaceList(rows,cols): #returns list of mesh faces for a rectangular mesh grid faceList=[] for j in range(rows-1): for i in range(cols-1): face=Rhino.Geometry.MeshFace(cols*j+i,cols*(j+1)+i,cols*(j+1)+i+1,cols*j+i+1) faceList.append(face) return faceList def ObjectsToNewLayer(objs,name,color=DefaultLayerColor): if not rs.IsLayer(name): rs.AddLayer(name,color) rs.ObjectLayer(objs,name) def HeaderInfo(data): xExt=data[4]*(data[0]-1) yExt=data[4]*(data[1]-1) origin=str(data[2])+" , "+str(data[3]) keys=["Columns (X):","Rows (Y):","Cell size:","Origin:","Extents in X:","Extents in Y:"] vals=[str(data[0]),str(data[1]),str(data[4]),origin,str(xExt),str(yExt)] HeaderMsg="" for n in range(6): HeaderMsg+="{0} {1}\n".format(keys[n], vals[n].rjust(30)) return HeaderMsg def ImportASCIIGridFile(): #unit system warning if rs.UnitSystem() != 4: us=rs.UnitSystemName(singular=False,abbreviate=False).upper() msg="Document units are {}, not meters - continue?".format(us) resp=rs.MessageBox(msg,33) if resp !=1: return #File open strPath=rs.OpenFileName("ASC file to import", "ASCII Grid Files (*.asc)|*.asc") if not strPath: return objFile=open(strPath) if not objFile: return #Get first 6 lines in file (header) gv=[] for n in range(5): line=objFile.readline() if line == '\n': print "Unable to read header correctly" ; return if n < 2: gv.append(int(line.split()[-1])) else: gv.append(float(line.split()[-1])) """The following is necessary because the next line is optional and may not be present in file""" pos=objFile.tell() #mark current file position line=objFile.readline() if line == '\n': print "Unable to read header correctly" ; return elif str.lower(line.split()[0])=="nodata_value": gv.append(float(line.split()[-1])) else: gv.append(-9999.0) objFile.seek(pos,0) #go back to previous line #Present message box with header info cancel=rs.MessageBox(HeaderInfo(gv),1,"ASCII Grid file data") if cancel !=1: return #User input section msg="Point grid density? (1 = every point, 2 = every other point, etc.)" gridStep=rs.GetInteger(msg,1,1) if gridStep==None: return ndValue=rs.GetReal("Z height to use if data is missing",0) if ndValue==None: return #calculate numerical values for file header items cols=((gv[0]-1)//gridStep)+1 #integer division... rows=((gv[1]-1)//gridStep)+1 cellSize=gv[4] ulCorner=[gv[2],((gv[1]-1)*cellSize)+gv[3]] #Note: ascii Grid file DATA starts in the UPPER left corner !! #ask user to choose what type of output if "ASCII_Output_Type" in sc.sticky: user_output = sc.sticky["ASCII_Output_Type"] else: user_output = (False,True,False) items=[("PointCloud","No","Yes"),("Mesh","No","Yes"),("Surface","No","Yes")] geoType=rs.GetBoolean("Objects to output:",items,user_output) if geoType == None or sum(geoType)==0: return if geoType[0] and not geoType[1] and not geoType[2]: resp=rs.GetBoolean("Separate NoData values?",[("Separate","No","Yes")],[(False)]) if resp==None: return sepND=resp[0] else: sepND=False if geoType[2]: srfSpecs=SurfaceType(cols,rows) if srfSpecs==None: return #srfSpecs is a tuple:(type,degree) type: True=Interpolated; False=Control """Main file read section starts here Iterate over file line by line, iterate over each line extract Z values, add X and Y values for points, append to master list""" start=time() pts=[] ; ndPts=[] ; read_errors=0 print "Reading in file data and building point grid" rs.Prompt("Processing file...") j=0 if PC(): rs.StatusBarProgressMeterShow("File read",0,rows,False,True) for line in objFile: if j%gridStep==0: if line != '\n': #strip any whitespace at beginning/end lineDataList=line.strip().split() #create list from line data for i in range(0,len(lineDataList),gridStep): x=ulCorner[0]+i*cellSize y=ulCorner[1]-j*cellSize #trying to get a float can cause errors good_value=True try: z=float(lineDataList[i]) except: good_value=False if good_value: if z==gv[5]: if sepND: ndPts.append(Rhino.Geometry.Point3d(x,y,ndValue)) else: pts.append(Rhino.Geometry.Point3d(x,y,ndValue)) else: pts.append(Rhino.Geometry.Point3d(x,y,z)) else: read_errors+=1 #print lineDataList[i] #debug pts.append(Rhino.Geometry.Point3d(x,y,ndValue)) j+=1 #print j if PC(): rs.StatusBarProgressMeterUpdate(j,True) objFile.close() if PC(): rs.StatusBarProgressMeterHide() readTime=str(round(time()-start,2)) print "Done reading in file and building point grid after "+readTime+" seconds" #Create and add geometry objects to file and finish start=time() if geoType[0]: #add a point cloud rs.Prompt("Creating point cloud") pc = Rhino.Geometry.PointCloud(pts) obj=sc.doc.Objects.AddPointCloud(pc) ObjectsToNewLayer(obj,"ASCII_Point_Cloud") if sepND and ndPts != []: pc = Rhino.Geometry.PointCloud(ndPts) obj=sc.doc.Objects.AddPointCloud(pc) ObjectsToNewLayer(obj,"NoData Points",(255,0,0)) if geoType[1]: #add a mesh rs.Prompt("Creating mesh") faces=CreateFaceList(rows,cols) mesh=Rhino.Geometry.Mesh() for pt in pts: mesh.Vertices.Add(pt) for face in faces: mesh.Faces.AddFace(face) mesh.Normals.ComputeNormals() mesh.Compact() obj=sc.doc.Objects.AddMesh(mesh) ObjectsToNewLayer(obj,"ASCII_Quad_Mesh") if geoType[2]: #add a surface rs.Prompt("Creating surface") deg=srfSpecs[1] if srfSpecs[0]: srf=Rhino.Geometry.NurbsSurface.CreateThroughPoints(pts,rows,cols,deg,deg,False,False) obj=sc.doc.Objects.AddSurface(srf) ObjectsToNewLayer(obj,"ASCII_IP_Grid_Srf") else: srf=Rhino.Geometry.NurbsSurface.CreateFromPoints(pts,rows,cols,deg,deg) obj=sc.doc.Objects.AddSurface(srf) ObjectsToNewLayer(obj,"ASCII_CP_Grid_Srf") rs.ZoomBoundingBox(rs.BoundingBox(obj),all=True) if geoType[1] or geoType[2]: try: rs.ViewDisplayMode("Perspective","Rendered") except: pass sc.doc.Views.Redraw() sc.sticky["ASCII_Output_Type"] = (geoType[0],geoType[1],geoType[2]) msg="Done creating geometry after "+str(round(time()-start,2))+" seconds" if read_errors: msg+=" | Found {} read errors!".format(read_errors) print msg ImportASCIIGridFile()