Hello,
I am trying to create a mask using a custom projection based on data from a shapefile. The script has the form like this:
scale=1:100000000
center_longitude=-94
center_latitude=90
J_options="-Js${center_longitude}/${center_latitude}/90/60/${scale}"
west_latitude=25.5090414309402647
west_longitude=-129.8926656877455912
east_latitude=52.1656584759616919
east_longitude=-0.6335393365702657
R_options_ll="-R${west_longitude}/${west_latitude}/${east_longitude}/${east_latitude}r"
resolution=80
fileout=basal_mask_${resolution}.nc
ogr2ogr -f "GMT" sed_mask.gmt sediment_properties_edit.shp
gmt mapproject sed_mask.gmt ${J_options} ${R_options_ll} -C -Fe > sed_mask_proj.gmt
gmt grdmask sed_mask_proj.gmt ${R_options} -G${fileout}=ni -Nz -aZ=dis_mas
The issue I am having is that the mapproject call completely strips out the aspatial data from the input gmt formatted file. This means that the grdmask call does not work, because the information from the “dis_mas” column in the original shapefile is not preserved in the mapproject call. Perhaps I am missing something, but is there a flag that ensures the aspatial data is preserved?
Many thanks for any advice.
Hi,
I don’t think mapproject “sees” aspatial fields the way these are presented in .gmt files produced by ogr2ogr. Basically, mapproject is unable to preserve full .gmt file structure as far as I can understand.
One way would be to reproject directly in ogr2ogr call, using the -t_srs option as ogr2ogr preserves full source shapefile structure. A problem here is that org2ogr does not recognize target projection specified as -Js${center_longitude}/${center_latitude}/90/60/${scale}, it rather needs an EPSG code or a proj4 string or a WKT (well known text) specification. I am not aware of a method of automagically converting gmt projection to a format that gdal tools can recognize and not experienced enough to write a proj4 projection manually from your -Js... specification manually, but this is possible.
another possibility could be to populate z values directly in the ogr2ogr call if this is all you need here:
ogr2ogr -f "GMT" sed_mask.gmt sediment_properties_edit.shp -zfield dis_mas
this way mapproject should be able to preserve z values as the third coordinate on the data input lines rather than in the aspatial fields.
Hi, thank you for the response.
In the end I reprojected the shapefile to what I am using within QGIS, so that I could use the GMT formatted file directly in grdmask.