Skip to content

geopandas: Mapping int/int64 to int32 for OGR_GMT format #2592

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
seisman merged 13 commits into from
Aug 21, 2023
Merged

geopandas: Mapping int/int64 to int32 for OGR_GMT format #2592

Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
8cec0c4
geopandas: Mapping int/int64 to int32 for OGR_GMT format
seisman Jul 4, 2023
eacf45b
Merge branch 'main' into ogrgmt-int32
seisman Aug 3, 2023
f2f10ca
Update the comment
seisman Aug 3, 2023
47bc26f
Merge branch 'main' into ogrgmt-int32
seisman Aug 5, 2023
f9b7903
Merge branch 'main' into ogrgmt-int32
seisman Aug 9, 2023
84f87a5
Fix a linting issue
seisman Aug 9, 2023
c1b3b43
Add unit test for plotting GeoDataFrame with int32 and int64 columns
weiji14 Aug 11, 2023
20891c1
Isort imports
weiji14 Aug 11, 2023
60ba52a
Change cmap from lajolla to lisbon
weiji14 Aug 19, 2023
a6d4315
Disable Int32 and Int64 dtype tests until geopandas 0.13.3 or greater
weiji14 Aug 19, 2023
1e05a29
Remove unused pandas import
weiji14 Aug 19, 2023
2097add
Force an index name and reset index
seisman Aug 20, 2023
a4ddcfe
Merge branch 'main' into ogrgmt-int32
seisman Aug 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 11 additions & 2 deletions pygmt/helpers/tempfile.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,18 +126,27 @@ def tempfile_from_geojson(geojson):
E.g. '1a2b3c4d5e6.gmt'.
"""
with GMTTempFile(suffix=".gmt") as tmpfile:
# pylint: disable=import-outside-toplevel
import geopandas as gpd

os.remove(tmpfile.name) # ensure file is deleted first
ogrgmt_kwargs = {"filename": tmpfile.name, "driver": "OGR_GMT", "mode": "w"}
try:
# Map int/int64 to int32 since OGR_GMT only supports 32-bit integer
# https://github.com/geopandas/geopandas/issues/967#issuecomment-842877704
# https://github.com/GenericMappingTools/pygmt/issues/2497
schema = gpd.io.file.infer_schema(geojson)

@weiji14 weiji14 Jul 4, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, looking into this infer_schema function and geopandas/geopandas#967 (comment) a bit more, I almost think part of the fix should go upstream into modifying infer_schema to output int32 values when the column is of int32 dtype, instead of forcing things to int/int64.

Of course, we'll still need to force the OGR_GMT output to int32 in PyGMT, but the code would be much simplified with some upstream fixes.

@seisman seisman Jul 4, 2023 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, looking into this infer_schema function and geopandas/geopandas#967 (comment) a bit more, I almost think part of the fix should go upstream into modifying infer_schema to output int32 values when the column is of int32 dtype, instead of forcing things to int/int64.

There is already an open PR geopandas/geopandas#2265, but it has been inactive for one year.

@weiji14 weiji14 Jul 4, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is already an open PR geopandas/geopandas#2265, but it has been inactive for one year.

Ok, and I saw geopandas/geopandas#2464 as well, also inactive 😆

Started a small PR at geopandas/geopandas#2950 to just handle int32, let's see if their test suite passes, I took a look at https://fiona.readthedocs.io/en/latest/manual.html#keeping-schemas-simple and fiona.FIELD_TYPES_MAP seem to have the following mapping:

{'int32': int,
 'float': float,
 'str': str,
 'date': fiona.rfc3339.FionaDateType,
 'time': fiona.rfc3339.FionaTimeType,
 'datetime': fiona.rfc3339.FionaDateTimeType,
 'bytes': bytes,
 'int64': int,
 'int': int,
 'List[str]': typing.List[str]}

So not sure if the geopandas default driver (ESRI Shapefile?) supports int32 schema. It it doesn't and the test suite breaks, then we might need to handle everything fully in PyGMT.

@seisman seisman Jul 6, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's still unclear if there will be any changes in GDAL, GeoPandas, Fiona or GMT, but as we are still using old versions of these packages, I believe we still need these workarounds in the next few years, right?

@weiji14 weiji14 Aug 11, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So the patch in geopandas has been merged at geopandas/geopandas#2950, but will need to wait for geopandas v0.13.3 or newer to be released. Even so, int64 will still need to be converted to int32, so yes, we'll need to keep this workaround for a while.

for col, dtype in schema["properties"].items():
if dtype in ("int", "int64"):
schema["properties"][col] = "int32"
ogrgmt_kwargs["schema"] = schema

@weiji14 weiji14 Aug 19, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only write schema if schema["properties"] is not None. Should fix that error about ValueError: Record does not match collection schema: KeysView(<fiona.model.Properties object at 0x7f3265f03a10>) != []

Suggested change
ogrgmt_kwargs["schema"] = schema
if schema["properties"]:
ogrgmt_kwargs["schema"] = schema

@seisman seisman Aug 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it's not the root cause of the error.

Here is a modified version of the test_geopandas_info_geodataframe test to understand why it fails:

import geopandas as gpd
import shapely

linestring = shapely.geometry.LineString([(20, 15), (30, 15)])
polygon = shapely.geometry.Polygon([(20, 10), (23, 10), (23, 14), (20, 14)])
multipolygon = shapely.geometry.shape(
    {
        "type": "MultiPolygon",
        "coordinates": [
            [
                [[0, 0], [20, 0], [10, 20], [0, 0]],  # Counter-clockwise
                [[3, 2], [10, 16], [17, 2], [3, 2]],  # Clockwise
            ],
            [[[6, 4], [14, 4], [10, 12], [6, 4]]],  # Counter-clockwise
            [[[25, 5], [30, 10], [35, 5], [25, 5]]],
        ],
    }
)
# Multipolygon first so the OGR_GMT file has @GMULTIPOLYGON in the header
gdf = gpd.GeoDataFrame(
    index=["multipolygon", "polygon", "linestring"],
    geometry=[multipolygon, polygon, linestring],
)

schema = gpd.io.file.infer_schema(gdf)
print(gdf)
print(schema)

gdf.to_file("test.gmt", driver="OGR_GMT", mode="w")

The output of the two print statements is:

                                                       geometry
multipolygon  MULTIPOLYGON (((0.00000 0.00000, 20.00000 0.00...
polygon       POLYGON ((20.00000 10.00000, 23.00000 10.00000...
linestring    LINESTRING (20.00000 15.00000, 30.00000 15.00000)
{'geometry': ['MultiPolygon', 'Polygon', 'LineString'], 'properties': OrderedDict()}

To understand what geopandas/fiona does, I also added the same two print statements to source code geopandas/io/file.py at line https://github.com/geopandas/geopandas/blob/b4b10313ab57bf2c55592a28fb99687c9a538fc2/geopandas/io/file.py#L591.

The two print statements in the geopandas source code give:

          index                                           geometry
0  multipolygon  MULTIPOLYGON (((0.00000 0.00000, 20.00000 0.00...
1       polygon  POLYGON ((20.00000 10.00000, 23.00000 10.00000...
2    linestring  LINESTRING (20.00000 15.00000, 30.00000 15.00000)
{'geometry': ['MultiPolygon', 'Polygon', 'LineString'], 'properties': OrderedDict([('index', 'str')])}

As you can see, both the gdf object and the schema are different. This is because geopandas by defaults modifies the index column (see lines https://github.com/geopandas/geopandas/blob/b4b10313ab57bf2c55592a28fb99687c9a538fc2/geopandas/io/file.py#L551-L556).

In our case, our inferred schema has a empty property dict OrderedDict(), but geopandas/fiona expects a property dict like OrderedDict([('index', 'str')]).

@seisman seisman Aug 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.to_file.html#geopandas.GeoDataFrame.to_file

index: bool, default None

If True, write index into one or more columns (for MultiIndex). Default None writes the index into one or more columns only if the index is named, is a MultiIndex, or has a non-integer data type. If False, no index is written.

I think writing index into one column or not has no effects to the output of OGR_GMT driver, so maybe we should just set index=False.

@weiji14 weiji14 Aug 19, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think writing index into one column or not has no effects to the output of OGR_GMT driver, so maybe we should just set index=False.

Setting index=True (the default) will produce this in the output *.gmt file:

# @VGMT1.0
# @R0/35/0/20                                                             
# @GMULTIPOLYGON
# @Nindex
# @Tstring
# FEATURE_DATA
>
# @Dmultipolygon
# @P
0 0
20 0
10 20
0 0
>
# @H
3 2
10 16
17 2
3 2
>
# @P
6 4
14 4
10 12
6 4
>
# @P
25 5
30 10
35 5
25 5
>
# @Dpolygon
# @P
20 10
23 10
23 14
20 14
20 10
>
# @Dlinestring
20 15
30 15

Notice the lines with the comments # @Nindex and # @Dmultipolygon which are from the index. So the 'index' is output as just a regular column in the OGR_GMT file, and could potentially be used by someone doing aspatial="Z=index". We probably do need to keep it just in case.

Is there another way to handle this? E.g. cast the 'index' column explicitly to a string type?

@seisman seisman Aug 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there another way to handle this? E.g. cast the 'index' column explicitly to a string type?

This won't always work, because the index property in a GeoDataFrame object can be integers. For example, in the above example, the object can be defined like

gdf = gpd.GeoDataFrame(
    #index=["multipolygon", "polygon", "linestring"],
    index=[0, 1, 2],
    geometry=[multipolygon, polygon, linestring],
)

@weiji14 weiji14 Aug 19, 2023 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

    if index:
        df = df.reset_index(drop=False)

Probably just need these two lines, call reset_index if an index exists?

@seisman seisman Aug 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoDataFrame.to_file.html#geopandas.GeoDataFrame.to_file

index has a default value of None. So df.reset_index is only called if the GeoDataFrame object's index is named or is non-integer.

@weiji14 weiji14 Aug 19, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah sorry, I thought index was referring to the df.index attribute. We could do something like this to force an index name:

if gdf.index.name is None:
    gdf.index.name = "index"
    # print(gdf.index.names) 
    # FrozenList(['index'])

@seisman seisman Aug 19, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It looks good.

@seisman seisman Aug 20, 2023 •
edited
Loading

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah sorry, I thought index was referring to the df.index attribute. We could do something like this to force an index name:

if gdf.index.name is None:
    gdf.index.name = "index"
    # print(gdf.index.names) 
    # FrozenList(['index'])

Patch applied in 2097add.

# Using geopandas.to_file to directly export to OGR_GMT format
geojson.to_file(**ogrgmt_kwargs)
except AttributeError:
# pylint: disable=import-outside-toplevel
# Other 'geo' formats which implement __geo_interface__
import json

import fiona
import geopandas as gpd

with fiona.Env():
jsontext = json.dumps(geojson.__geo_interface__)
Expand Down
4 changes: 4 additions & 0 deletions pygmt/tests/baseline/test_geopandas_plot_int_dtypes.png.dvc
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
outs:
- md5: 6ce1b73d25ac10900a2ce4c7148d2953
size: 43434
path: test_geopandas_plot_int_dtypes.png
44 changes: 43 additions & 1 deletion pygmt/tests/test_geopandas.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,9 @@
Test integration with geopandas.
"""
import numpy.testing as npt
import pandas as pd
import pytest
from pygmt import Figure, info
from pygmt import Figure, info, makecpt, which

gpd = pytest.importorskip("geopandas")
shapely = pytest.importorskip("shapely")
Expand Down Expand Up @@ -131,3 +132,44 @@ def test_geopandas_plot3d_non_default_circle():
style="c0.2c",
)
return fig


@pytest.mark.parametrize("dtype", ["int32", "int64", pd.Int32Dtype(), pd.Int64Dtype()])
@pytest.mark.mpl_image_compare(filename="test_geopandas_plot_int_dtypes.png")
def test_geopandas_plot_int_dtypes(dtype):
"""
Check that plotting a geopandas GeoDataFrame with integer columns works,
including int32 and int64 (non-nullable), Int32 and Int64 (nullable).

This is a regression test for
https://github.com/GenericMappingTools/pygmt/issues/2497
"""
# Read shapefile in geopandas.GeoDataFrame
shapefile = which(
fname="@RidgeTest.shp @RidgeTest.shx @RidgeTest.dbf @RidgeTest.prj",
download="c",
)
gdf = gpd.read_file(shapefile[0])

# Reproject geometry and change dtype of NPOINTS column
gdf["geometry"] = (
gdf.to_crs(crs="EPSG:3857")
.buffer(distance=100000)
.to_crs(crs="OGC:CRS84") # convert to lon/lat to prevent @null in PROJ CRS
)
gdf["NPOINTS"] = gdf.NPOINTS.astype(dtype=dtype)

# Plot figure with three polygons colored based on NPOINTS value
fig = Figure()
makecpt(cmap="lajolla", series=[10, 60, 10], continuous=True)

@weiji14 weiji14 Aug 11, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We might need to check what happens when a big int64 number $2^{32} &gt; x &gt; 2^{64}$ is passed in. Don't want to accidentally convert and plot data incorrectly!

@weiji14 weiji14 Aug 20, 2023 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@seisman, could you add a test for a large int64 number too? Can use a simple info just to check the cast to int32 behaviour, want to make sure it doesn't truncate int64 numbers to int32.

@seisman seisman Aug 20, 2023

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Modified the existing test with large numbers:

import geopandas as gpd
from pygmt import which, info

shapefile = which(                                                          
    fname="@RidgeTest.shp @RidgeTest.shx @RidgeTest.dbf @RidgeTest.prj",       
    download="c",                                                           
)                                                                           
gdf = gpd.read_file(shapefile[0])

gdf["geometry"] = (                                                         
    gdf.to_crs(crs="EPSG:3857")                                             
    .buffer(distance=100000)                                                
    .to_crs(crs="OGC:CRS84")  # convert to lon/lat to prevent @null in PROJ CRS
)

gdf["NPOINTS"] *= 1.0e10
gdf["NPOINTS"] = gdf.NPOINTS.astype(dtype="int64")
info(gdf)

It raises an OverflowError error:

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
Cell In[12], line 1
----> 1 info(gdf)

File ~/OSS/gmt/pygmt/pygmt/helpers/decorators.py:598, in use_alias.<locals>.alias_decorator.<locals>.new_module(*args, **kwargs)
    591     msg = (
    592         "Parameters 'Y' and 'yshift' are deprecated since v0.8.0. "
    593         "and will be removed in v0.12.0. "
    594         "Use Figure.shift_origin(yshift=...) instead."
    595     )
    596     warnings.warn(msg, category=SyntaxWarning, stacklevel=2)
--> 598 return module_func(*args, **kwargs)

File ~/OSS/gmt/pygmt/pygmt/helpers/decorators.py:738, in kwargs_to_strings.<locals>.converter.<locals>.new_module(*args, **kwargs)
    736             kwargs[arg] = separators[fmt].join(f"{item}" for item in value)
    737 # Execute the original function and return its output
--> 738 return module_func(*args, **kwargs)

File ~/OSS/gmt/pygmt/pygmt/src/info.py:85, in info(data, **kwargs)
     83 file_context = lib.virtualfile_from_data(check_kind="vector", data=data)
     84 with GMTTempFile() as tmpfile:
---> 85     with file_context as fname:
     86         lib.call_module(
     87             module="info",
     88             args=build_arg_string(kwargs, infile=fname, outfile=tmpfile.name),
     89         )
     90     result = tmpfile.read()

File ~/opt/miniconda/lib/python3.10/contextlib.py:135, in _GeneratorContextManager.__enter__(self)
    133 del self.args, self.kwds, self.func
    134 try:
--> 135     return next(self.gen)
    136 except StopIteration:
    137     raise RuntimeError("generator didn't yield") from None

File ~/OSS/gmt/pygmt/pygmt/helpers/tempfile.py:147, in tempfile_from_geojson(geojson)
    145     ogrgmt_kwargs["schema"] = schema
    146     # Using geopandas.to_file to directly export to OGR_GMT format
--> 147     geojson.to_file(**ogrgmt_kwargs)
    148 except AttributeError:
    149     # Other 'geo' formats which implement __geo_interface__
    150     import json

File ~/opt/miniconda/lib/python3.10/site-packages/geopandas/geodataframe.py:1263, in GeoDataFrame.to_file(self, filename, driver, schema, index, **kwargs)
   1172 """Write the ``GeoDataFrame`` to a file.
   1173 
   1174 By default, an ESRI shapefile is written, but any OGR data source
   (...)
   1259 
   1260 """
   1261 from geopandas.io.file import _to_file
-> 1263 _to_file(self, filename, driver, schema, index, **kwargs)

File ~/opt/miniconda/lib/python3.10/site-packages/geopandas/io/file.py:572, in _to_file(df, filename, driver, schema, index, mode, crs, engine, **kwargs)
    569     raise ValueError(f"'mode' should be one of 'w' or 'a', got '{mode}' instead")
    571 if engine == "fiona":
--> 572     _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs)
    573 elif engine == "pyogrio":
    574     _to_file_pyogrio(df, filename, driver, schema, crs, mode, **kwargs)

File ~/opt/miniconda/lib/python3.10/site-packages/geopandas/io/file.py:601, in _to_file_fiona(df, filename, driver, schema, crs, mode, **kwargs)
    597     crs_wkt = crs.to_wkt("WKT1_GDAL")
    598 with fiona.open(
    599     filename, mode=mode, driver=driver, crs_wkt=crs_wkt, schema=schema, **kwargs
    600 ) as colxn:
--> 601     colxn.writerecords(df.iterfeatures())

File ~/opt/miniconda/lib/python3.10/site-packages/fiona/collection.py:558, in Collection.writerecords(self, records)
    556 if self.mode not in ("a", "w"):
    557     raise OSError("collection not open for writing")
--> 558 self.session.writerecs(records, self)
    559 self._len = self.session.get_length()
    560 self._bounds = None

File fiona/ogrext.pyx:1409, in fiona.ogrext.WritingSession.writerecs()

File fiona/ogrext.pyx:441, in fiona.ogrext.OGRFeatureBuilder.build()

OverflowError: value too large to convert to int

@weiji14 weiji14 Aug 21, 2023

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, good to know that it errors. We should probably cast int64 numbers to float64 or string dtype, xref geopandas/geopandas#2950 (comment) and https://github.com/OSGeo/gdal/blob/6e9103bd5acb1ff8da305c4e77fa30335ffa3b70/ogr/ogrsf_frmts/gmt/ogrgmtlayer.cpp#L754-L773. We can do this in a follow-up PR since this one is getting a bit long.

fig.plot(
data=gdf,
frame=True,
pen="1p,black",
close=True,
fill="+z",
cmap=True,
aspatial="Z=NPOINTS",
)
fig.colorbar()
return fig