Forum Replies Created
-
AuthorPosts
-
May 25, 2018 at 7:39 pm in reply to: Obtain array of Station/Elevation for a cross section in RAS Controller #13590TobyParticipant
Hi,
You can accomplish this by using the .g01.hdf or .p01.hdf file.
Here is an R script to get you started. This works in 5.0.3. You would need to update a line or two to get it to work in 5.0.4.
cheers,
################################################
# To access HDF5 files in R, we will use the rhdf5 library which is part of the
#Bioconductor suite of R libraries.#install rhdf5 package
#source(“http://bioconductor.org/biocLite.R”)
#biocLite(“rhdf5”)# load libraries
library(raster)
library(rhdf5)
library(rgdal)
library(sp)
library(maptools)
library(rgeos)# #####################################################################################################
# # Cross-section compairson
#
# # set working directory
# setwd(“C:\\Users\\yourComp\\Documents\\myRAS”)
#
# # define the file name as an object
# plan<-"model.g01.hdf" # # proj 4 string - look to http://www.spatialreference.org
# sp.Ref<- "+proj=lcc +lat_1=38.33333333333334 +lat_2=39.83333333333334 +lat_0=37.66666666666666 +lon_0=-122 # +x_0=2000000 +y_0=500000.0000000001 +datum=NAD83 +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0" # # ##################################################################################################### # # Use of function... # # ras.data<-RAS_XS_to_R(plan, SP.Ref=NULL) # # #Unpack xs data # xs_lines<-ras.data[[1]] # # #Unpack profile data # profile<-ras.data[[2]] # ##################################################################################################### #function RAS_XS_to_R<-function(plan, sp.Ref=NULL){ ##################################################################################################### #Extract xs metadata path<-"Geometry/Cross Sections/River Names" riverNames <- t(h5read(plan, path)) #trim trailing whitespace riverNames<-sub("\\s+$", "", riverNames) path<-"Geometry/Cross Sections/Reach Names" reachNames <- t(h5read(plan, path)) #trim trailing whitespace reachNames<-sub("\\s+$", "", reachNames) path<-"Geometry/Cross Sections/River Stations" riverStations <- t(h5read(plan, path)) #trim trailing whitespace riverStations<-sub("\\s+$", "", riverStations) path<-"Geometry/Cross Sections/Node Descriptions" nodeDesc <- t(h5read(plan, path)) #remove \r\n nodeDesc<-as.character(gsub("\r\n"," ",nodeDesc)) # create reach names xsNames<-paste0(riverNames," ",reachNames," ",riverStations) #find the River channel reach lengths path<-"Geometry/Cross Sections/Lengths" reachLengths <- t(h5read(plan, path)) ##################################################################################################### # Get geometry # get xs point count path<-"Geometry/Cross Sections/Polyline Parts" xs.parts <- h5read(plan, path) xs.parts.sum<-cumsum(xs.parts[2,]) # get xs coords path<-"Geometry/Cross Sections/Polyline Points" xs.points <- t(h5read(plan, path)) nrow(xs.points) ##################################################################################################### #Create XS lines tmp<-xs.points[1:xs.parts.sum[1],] xs_lines<-SpatialLines(list(Lines(list(Line(tmp)), ID = "1"))) for(i in 2:length(xs.parts.sum)){ tmp<-xs.points[(xs.parts.sum[i-1]+1):xs.parts.sum[i],] L.tmp <- SpatialLines(list(Lines(list(Line(tmp)), ID = as.character(i)))) xs_lines <- spRbind(xs_lines, L.tmp) }# close i loop #create the spatial lines xs_lines<-SpatialLinesDataFrame(xs_lines, data.frame(xsNames=xsNames,rivNames=t(riverNames),reaNames=t(reachNames), rivSta=t(riverStations),Desc=nodeDesc,ReachLen=reachLengths)) # Assign the sp ref if provided if(!is.null(sp.Ref)){ crs(xs_lines)<-sp.Ref } ##################################################################################################### # get cross section station and elevation data from RAS path<-"Geometry/Cross Sections/Station Elevation Info" eleInfo <- t(h5read(plan, path)) path<-"Geometry/Cross Sections/Station Elevation Values" eleValue <- t(h5read(plan, path)) # convert profile values to dataframe profile<-data.frame(station=eleValue[,1], elev=eleValue[,2]) # assign each profile an id profile$id<-NA for(i in 1:nrow(eleInfo)){ startID<-eleInfo[i,1] endID<-startID + eleInfo[i,2] profile$id[(startID+1):endID]<-i }#close i loop # give profile names ids xsNames<-data.frame(id=1:length(xsNames), names=xsNames) profile$xsName<-xsNames[match(profile$id,xsNames$id),"names"] # xs stationing may not start at zero!!! # correct for xs which do not start at 1 for(i in 1:length(unique(profile$id))){ index<-which(profile$id==i) if(profile[index,1][1]!=0){profile[index,1]<-profile[index,1]-profile[index,1][1]} }# close i loop return(list(xs_lines, profile)) }# close functionTobyParticipantHi Tong Zhai,
Thanks for sharing this info! Have you been modifying the *.rasmap file to update your terrain and land class data?
thanks,
TobyParticipantIF you link the cell face points with their indexes you will get the right count found in
Faces FacePoint Indexes. Found in Geom/2D flow/mesh name
TobyParticipantYou can get facepoint coords in the HDF file at “Geometry/2D Flow Areas/”,your mesh name,”/FacePoints Coordinate”
TobyParticipantHowdy,
Has anyone in the last few years attempted this?
It would be great to use RAS mapper’s utilities in another program, such as program R.
thanks,
TobyParticipantI have tricked RAS before with the issue, by deleting the spatial reference in RAS mapper, and then creating the terrain, then assigning the sp. ref. again. I have had your error, when I know the terrain data and XS data should match, and they end up getting moved elsewhere. Just don’t let RAS mapper associate a projection with the document as you add your raster files.
Hope this helps.
TobyParticipantYou could try and draw your cross sections so they snap together at the ends and put a lateral structure connecting the two reaches.
TobyParticipantThe WSE values match the cell centers data.
I use program R to make shapefiles from the HDF files, but that gets a bit tricky if you have never used R.
In the HDF viewer you can export a text file. Play around with it. Then merge the xy data and the velocity or WSE data in excel. This is a straight forward process you should be able to figure out. You will need to transpose the data to get the xy data to line up with the velocity or WSE data. Off the top of my head I can’t remember, but your runtime length should match your columns or rows and the data should match the other column or row count. Explore the HDF files….lots of good stuff in there!
runtime = “Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series/Time Date Stamp”
cellXY = “Geometry/2D Flow Areas/”,mesh,”/Cells Center Coordinate”
TobyParticipantHi Tom,
To get the velocity data make sure you set it in the plan output options. Under the plan window -> Options -> Output Options… -> HDF5 Write Parameters -> check the Write Velocity data box.
In the HDF viewer you can Export the data to a text file. You can then use excel to get your data formatted into the right setup for an XY data import to Arc Map.
http://support.esri.com/en/technical-article/000012745
Good Luck,
TobyParticipantYou might consider using the HDF viewer and copying csv tables from there to create a shapefile in Arc. You will have to set the HDF parameters in your plan options to write the velocity data.
The data you want is under the following drop downs in the HDF file:
Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series/2D Flow Areas/”your mesh”/Node X Vel
Results/Unsteady/Output/Output Blocks/Base Output/Unsteady Time Series/2D Flow Areas/your mesh/Node Y Vel
Geometry/2D Flow Areas/”your mesh”/FacePoints Coordinate
https://support.hdfgroup.org/downloads/index.html
you want the viewer – HDFView 2.14
TobyParticipantI bet you are using Crayfish?….
Make your max depth map in RAS Mapper.
Select depth in RAS mapper right click, edit map parameters -> select maximum, then save to disk. It will save in a new folder in your model directory.
cheers,
TobyParticipantGlad the fix worked for you too. I have not heard anything back from HEC about the bug, but hopefully they are fixing it or, set the pre-processor to trigger an error so you know if you have a split cell in your mesh/1d boundary.
Cheers!
TobyParticipantI have had lots of issues with the 1d / 2d connections.
Check out this post I did and check for these errors first. Your LS’s look like they are snapped to the XS as mine where.
I have also found areas with high velocities at the flow boundaries, can cause mass balance issues.
Good luck. It would be interesting to know if you find the same issues I had. I have been in contact with HEC and they think there might be a bug, so I sent them my model files.
TobyParticipantIf the direction if flow is correct (eg water is flowing uphill), what is your inflow?
You may want to use the unsteady engine and run your inflow hydrograph as a “fake” steady state in the unsteady flow editor. Then run your model over a time window and see if water level changes. Right now your WSE is sloped at the downstream end, but you expect that to translate upstream further is what I am guess you are thinking.
TobyParticipantYour river direction is backwards. In RAS station 0 is the downstream station. You can see what looks like a normal depth BC at station 0. So water has to fill the cross sections, like it would in a lake, before spilling out. Hope this fixes it for you.
-
AuthorPosts