Add projection information to DataArrays read from GRD files

[santisoler]

Description of the desired feature:

Right before merging #348 @mdtanker noticed that the load_oasis_grd() function doesn't preserve information about the grid projection. Even though the file header has a few bytes reserved for storing reference to map projection, this bytes seem to be always equal to zero. Nevertheless, the projection information is available in the .xml files that accompany the .grd files. @mdtanker wrote a snippet that can read those files and extract the information about the grid projection.

Would be nice to add the load_oasis_grd() function the capability to read the xml files and preserve the information about the grid projection. The two things that should be preserved are:

• they projection code/type (EPSG code, UTM zone, etc)
• the datum (e.g. WGS84)

Since the xml should have the same filename as the .grd file we can assume their name. But the function should be able to read the .grd file even if the xml file doesn't exist. Maybe raising a warning would be nice.

Are you willing to help implement and maintain this feature? Yes, but would love to see contributors to tackle this one, specially if they have access to Oasis Montaj license and can generate sample files to test this feature.

[markjessell]

It may be worth pointing out that the *.xml is an optional part of the geosoft system, montaj reads the projection info from the *.gi file, not the *.xml file. This doesn't matter for import, but if anyone wants to make an exporter to *.grd format then this will be an issue. In any case, for import, the xml reader would be very handy when that file is available, but we shouldn't assume it will always be available. Evren Pakyuz Charrier of Oslandia is still talking to Seequent about getting the *.gi format.

[santisoler]

Thanks @markjessell for your input. As you mentioned, the .gi format is not available and the .xml file is the easiest way to get information about the projection. I totally agree that the .xml file should be optional: the function should work as it currently does even if the .xml file doesn't exist.

Since we are not planning on making a writer to GRD files, not having the .gi format is not a deal breaker, but would be nice to have it so we could extract more information while reading the files. We don't plan on making a writer because we think there are way better file formats for multidimensional arrays out there: like netCDF which is an open format and it's being used by many other pieces of software (GMT, Xarray, rasterio, gdal, etc), allowing a smoother interaction between our software and the rest of the ecosystem. It was a major obstacle having GRD files and not being able to read them, but once we have them in Python we already have amazing tools to store them in files.

Evren Pakyuz Charrier of Oslandia is still talking to Seequent about getting the *.gi format.

I hope these conversations to be fruitful!

[mtb-za]

Briefly spoke to @santisoler about getting involved again, and agreed that this issue was a sensible one to get out of the way.

Did Evren get somewhere with Seequent regarding reading .gi files as well, or are we sticking with just the xml (if it exists) for now? Do we have some sample data to use for testing?

[markjessell]

I can ask Evren again, but no he didn’t give me any info. Funny, I was just going to look at the binary of the gi file to see if the epsg info is in a fixed location or easy to identify. I’ll let you know how I go. I have lots of sample data that is public domain what are you looking for specifically. I also have a Geosoft licence so can generate files as needed. Mark Jessell Professor/ Western Australian Fellow Centre for Exploration Targeting/ SES • M006, Perth WA 6009 Australia T +61 8 6488 5803 • E ***@***.***

…

________________________________ From: Martin Bentley ***@***.***> Sent: Saturday, March 9, 2024 8:53:24 AM To: fatiando/harmonica ***@***.***> Cc: Mark Jessell ***@***.***>; Mention ***@***.***> Subject: Re: [fatiando/harmonica] Add projection information to DataArrays read from GRD files (Issue #365) Briefly spoke to @santisoler<https://github.com/santisoler> about getting involved again, and agreed that this issue was a sensible one to get out of the way. Did Evren get somewhere with Seequent regarding reading .gi files as well, or are we sticking with just the xml (if it exists) for now? Do we have some sample data to use for testing? — Reply to this email directly, view it on GitHub<#365 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABUBSYGS4V3Q4NUGTPCTQE3YXJMQJAVCNFSM6AAAAAASLSPEU2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSOBWGY2DGOJZGM>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[mtb-za]

Do not need anything dramatic, just a representative sample of the xml files to make sure that we are reading in the right bits of it. If the projection information is consistent then I only need one or two to make sure that I am reading the correct fields, but if there are possible variations in the fields, then examples of those (or what the variation might be).

[santisoler]

I remember that @LL-Geo included some .xml files along with the .grd files he added for testing purposes. You can find them here: https://github.com/fatiando/harmonica/tree/main/harmonica/tests/data

[mtb-za]

Took a look at those, and none of them have projection data. They all have

<georeference>
  <dataextents>
	  <extent2d minx="1" miny="-24" maxx="50" maxy="24"/>
  </dataextents>
  <projection type="UNKNOWN">
	  <units name="m" unit_scale="1"/>
  </projection>
</georeference>

or similar. So we will need to get some example xml files that have this information, or re-export these to have it.

It looks like we will be checking (at least?) three things in order:

1. Whether an xml file exists at all
2. Whether the xml file contains a georeference field
3. Whether the georeference field has anything useful in it. We might only care about projection specifically?

Does that seem like a sensible approach?

[santisoler]

Oh, right. I think those files didn't have a proper projection, although having "UNKNOWN" projection is something. I think we should still use them for testing this new feature.

In order to find an xml file that contains actual projections, I think NRCAN is a good resource: https://geophysical-data.canada.ca/Portal/Results

All datasets in there are open access. And actually, the load_oasis_montaj_grid function was born as a need to read a few of those files.

[mtb-za]

I grabbed a couple of GRD files from there and found the fields that we (probably) care about.

<geo:georeference>
	<geo:dataextents>
		<geo:extent2d minx="499100" miny="5402250" maxx="519625" maxy="5413700"/>
	</geo:dataextents>
	<geo:projection type="TRANSVERSE_MERCATOR" name="NAD83 / UTM zone 21N" wellknown_epsg="26921" projection="UTM zone 21N" ellipsoid="GRS 1980" datum="NAD83" datumtrf="NAD83 to WGS 84 (1)" datumtrf_epsg="1188" datumtrf_description="[NAD83] (4m) North America - Canada and USA (CONUS, Alaska mainland)" radius="6378137" eccentricity="0.08181919104">
		<geo:shift x="0" y="0" z="0"/>
		<geo:parameters scale="0.9996" false_east="500000" false_north="0" lat_origin="0" lon_origin="-57"/>
		<geo:units name="m" unit_scale="1"/>
	</geo:projection>
</geo:georeference>

It also looks like we should be able to save the projection as a WKT when saving a netCDF from xarray. This will need some conversion from the string(s) that look to be in the xml file, but this should be doable. Just need to get my head around the cfconventions, but at least proj should be able to generate a proper WKT. Not sure if we only need the epsg code for that, or if we will need more of this information.

Failing that, I can save the above as a sensible dictionary within the xarray and let it handle things as it wishes.

Do we want to do anything with the data extents, assuming that they exist?

[LL-Geo]

Well done! 👍

We can start with the extract of the "wellknown_epsg=26921", similar like the def extract_proj_str(fname) https://github.com/markjessell/grd_loader/blob/main/geosoft_grid_parser.py.

If that exists, we can read that. If not, we can ask an input for the function. load_oasis_grd(, epsj='')

For the data extents, these can be read from the grd file. So we can ignore these from the XML file.

[markjessell]

Here are some real world XMLs https://www.dropbox.com/scl/fi/nz56vhq5i7xrnso9b61ds/grd.xmls.zip?rlkey=kobjxpb7zf1d2a2wcyvmtx3o5&dl=0 [https://www.dropbox.com/static/metaserver/static/images/opengraph/opengraph-content-icon-file-zip-landscape.png]<https://www.dropbox.com/scl/fi/nz56vhq5i7xrnso9b61ds/grd.xmls.zip?rlkey=kobjxpb7zf1d2a2wcyvmtx3o5&dl=0> grd.xmls.zip<https://www.dropbox.com/scl/fi/nz56vhq5i7xrnso9b61ds/grd.xmls.zip?rlkey=kobjxpb7zf1d2a2wcyvmtx3o5&dl=0> Shared with Dropbox www.dropbox.com I think they all follow the same logic, but in the version I made for the qgis plugin (https://github.com/Loop3D/grd_loader) that uses extract_proj_str() I did come across one XML that was a different format (but can't find it now of course, I'll keep looking).

…

________________________________ From: Lu Li ***@***.***> Sent: Wednesday, 13 March 2024 9:21 AM To: fatiando/harmonica ***@***.***> Cc: Mark Jessell ***@***.***>; Mention ***@***.***> Subject: Re: [fatiando/harmonica] Add projection information to DataArrays read from GRD files (Issue #365) Well done! 👍 We can start with the extract of the "wellknown_epsg=26921", similar like the def extract_proj_str(fname) https://github.com/markjessell/grd_loader/blob/main/geosoft_grid_parser.py. If that exists, we can read that. If not, we can ask an input for the function. load_oasis_grd(, epsj='') For the data extents, these can be read from the grd file. So we can ignore these from the XML file. — Reply to this email directly, view it on GitHub<#365 (comment)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABUBSYAIAEN2LSETCK7ITPTYX6SYBAVCNFSM6AAAAAASLSPEU2VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSOJSHE2TKNZYG4>. You are receiving this because you were mentioned.Message ID: ***@***.***>