Daniel Neumann
2014-07-23 22:45:46 UTC
Hello all-together,
I have got a horizontal grid with about 24x24 km^2 grid resolution over the North Sea region (Lambert conformal conic projection). I would like to calculate per grid cell: (a) the area of ocean, (b) the area of land and, (c) the area of surfzone. As surfzone I denote a stripe along the coastline of 50m width. Calculating (c) is my main problem for the moment. (a) and (b) follow if (c) is solved.
In the end, the ration of grid cell surface which is covered by these three types is needed. Therefore, the unit of the calculated area is not important.
This data is required to create a so called OCEANfile for simulations with the Community Multiscale Air Quality (CMAQ).
My approach to calculate the surfzone area is:
# set all grid points within the surfzone to 1 and all other to 0
# (can be on the water OR on the land side)
grdmath -R-5/15/40/65 -I1 mycoast LDIST 0.05 LE = close_to_coast.nc
# set all ocean points to 1 and all other to zero
grdlandmask -R-5/15/40/65 -I1.0 -N1/0/0/0/0 -Df -Gocean_only.nc
# combine the two files from above in order to set all grid points
# to 1 which are in the ocean and close to the coastline
grdmath close_to_coast ocean_only.nc MUL = surfzone.nc
# Finally we
# - get intersection between each gridcell and the surfzone and
# - calculate area of surfzone by grdvolume.
However, the 1 degree resolution (-I1) used above is to coarse because the width of the surfzone is in the order of 50 m. 0.0001 degree seems to be more appropriate. A grid with this resolution that includes the North Sea consists of in the order of 10^10 grid points. Calling the first 3 expressions for 0.1 degree resolution takes already quite a long time ... .
To get the surfzone area per 24x24 km^2 grid cell I would iterate each grid cell (by a shell script), get the above calculated grid points which intersect with the current grid cell and, calculate the area by grdvolume. We have at least in the order of 1000 grid cells (rather 10,000) which takes quite a while.
This approach seems to be too time-consuming and complicated.
I though about adding a 50 m margin to the coastline and calculating and intersecting area between margin and ocean. Is this possible with GMT?
I would be grateful for hints and alternative solutions.
Regards,
Daniel
Helmholtz-Zentrum Geesthacht
Zentrum für Material- und Küstenforschung GmbH
Max-Planck-Straße 1 I 21502 Geesthacht I Deutschland/Germany
Geschäftsführer/Board of Management: Prof. Dr. Wolfgang Kaysser, Dipl.-Ing. Michael Ganß
Vorsitzender des Aufsichtsrates/Chairman of the Supervisory Board: MinDirig Wilfried Kraus
Amtsgericht Lübeck HRB 285 GE (Register Court)
Internet: http://www.hzg.de
Mailing list for GMT discussions of all kinds. If you are not sure you have found a bug, discuss it here first.
To formally report bugs or request features, please register and add New Issue on gmt.soest.hawaii.edu
To unsubscribe, send the message "signoff gmt-help" to ***@lists.hawaii.edu
Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum on gmt.soest.hawaii.edu instead.
I have got a horizontal grid with about 24x24 km^2 grid resolution over the North Sea region (Lambert conformal conic projection). I would like to calculate per grid cell: (a) the area of ocean, (b) the area of land and, (c) the area of surfzone. As surfzone I denote a stripe along the coastline of 50m width. Calculating (c) is my main problem for the moment. (a) and (b) follow if (c) is solved.
In the end, the ration of grid cell surface which is covered by these three types is needed. Therefore, the unit of the calculated area is not important.
This data is required to create a so called OCEANfile for simulations with the Community Multiscale Air Quality (CMAQ).
My approach to calculate the surfzone area is:
# set all grid points within the surfzone to 1 and all other to 0
# (can be on the water OR on the land side)
grdmath -R-5/15/40/65 -I1 mycoast LDIST 0.05 LE = close_to_coast.nc
# set all ocean points to 1 and all other to zero
grdlandmask -R-5/15/40/65 -I1.0 -N1/0/0/0/0 -Df -Gocean_only.nc
# combine the two files from above in order to set all grid points
# to 1 which are in the ocean and close to the coastline
grdmath close_to_coast ocean_only.nc MUL = surfzone.nc
# Finally we
# - get intersection between each gridcell and the surfzone and
# - calculate area of surfzone by grdvolume.
However, the 1 degree resolution (-I1) used above is to coarse because the width of the surfzone is in the order of 50 m. 0.0001 degree seems to be more appropriate. A grid with this resolution that includes the North Sea consists of in the order of 10^10 grid points. Calling the first 3 expressions for 0.1 degree resolution takes already quite a long time ... .
To get the surfzone area per 24x24 km^2 grid cell I would iterate each grid cell (by a shell script), get the above calculated grid points which intersect with the current grid cell and, calculate the area by grdvolume. We have at least in the order of 1000 grid cells (rather 10,000) which takes quite a while.
This approach seems to be too time-consuming and complicated.
I though about adding a 50 m margin to the coastline and calculating and intersecting area between margin and ocean. Is this possible with GMT?
I would be grateful for hints and alternative solutions.
Regards,
Daniel
Helmholtz-Zentrum Geesthacht
Zentrum für Material- und Küstenforschung GmbH
Max-Planck-Straße 1 I 21502 Geesthacht I Deutschland/Germany
Geschäftsführer/Board of Management: Prof. Dr. Wolfgang Kaysser, Dipl.-Ing. Michael Ganß
Vorsitzender des Aufsichtsrates/Chairman of the Supervisory Board: MinDirig Wilfried Kraus
Amtsgericht Lübeck HRB 285 GE (Register Court)
Internet: http://www.hzg.de
Mailing list for GMT discussions of all kinds. If you are not sure you have found a bug, discuss it here first.
To formally report bugs or request features, please register and add New Issue on gmt.soest.hawaii.edu
To unsubscribe, send the message "signoff gmt-help" to ***@lists.hawaii.edu
Note: gmt-help will become obsolete on Sept 1, 2014 - please use forum on gmt.soest.hawaii.edu instead.