New answers tagged

3

The principle The solution consists of a few simple steps: reproject your polygon to a Pacific-centered CRS and create a line for the antimeridian in the same CRS. Create a small buffer around this line and substract it with the Difference too from the polygon. See below for screenshots. How to do that in detail, step by step: Re-project your layer (polygon)...


2

A graph is said to be connected if there is a path between every pair of vertex. Therefore create a new full connected Graph with itertools import itertools Road = nx.read_shp('stac_graph.shp') # control Road2 = Road.to_undirected() nx.is_connected(Road2) False # new graph with path between every pair of vertex G = nx.Graph() # add original nodes G....


5

Geopandas 0.7 added a new rows parameter to read_file. You can use it to read the first n rows, or a specific slice of rows. import geopandas as gpd # Read the first 100 rows gdf = gpd.read_file("/path/to/my/shapefile.shp", rows=100) # Read the 5 rows from the 100000th gdf = gpd.read_file("/path/to/my/shapefile.shp", rows=slice(100000, ...


0

It helps if you explain your process and what you have tried. Otherwise, we don't have much to go off of. Have you seen the docs on how to do proportional symbols? Or you could just follow this very detailed answer to a similar question. I followed it and had no problems displaying the symbols in the legend. Ultimately, after completing the initial steps ...


0

I would implement the solution that is describe here Geometry based automatic fill of fields in attribute table The solution is essentially a trigger that will update the field automatically when ever a transaction occurs. This way, the data will not only populate the field upon data creation, but also if geometry changes. I would also move off of using a ...


0

You can also use package shapefile, which can be found here: https://github.com/GeospatialPython/pyshp import shapefile as shp sf = shp.Reader("percelen.shp") sf.bbox [41988.03170000017, 399987.0051999986, 47011.34809999913, 402276.9079999998]


2

There is no need to list maps etc. and I dont think you can clip a lyrx. Use ListFeatureClasses: import os, arcpy arcpy.env.workspace = r'C:\path\to\shapefilefolder' featureclass_to_be_clipped = r'C:\data.gdb\features123' out_folder= r'C:\outfolder' for shapefile in arcpy.ListFeatureClasses(): #List all shapefiles in env.workspace folder arcpy.analysis....


0

I have found the solution: with shapefile.Reader(shapefile_path) as shp: total_rows = shp.numRecords for row_num, row in enumerate(shp.iterRecords()): if row_num < 0: pass elif row_num > 100: break


0

A quick update to the answer. You may find that your data haven't been saved in UTF-8. You will need to use a different code in the cpg file. For example, Latin-9 requires the code ISO 88591 instead. Just that, one small line of text at the top of a text file. The encoding of this text file (the cpg file) can be UTF-8 or whatever the OS prefers. QGIS doesn't ...


4

shapefile.Reader doesn't support indexing. It loads all data to the memory. After using shapefile.Reader("file_path") you can get a specific range of features in different ways. Using itertools as in Thomas's answer. import shapefile import itertools shapefile_path = "C:/path/to/shapefile.shp" shp = shapefile.Reader(shapefile_path) ...


4

You can use Python itertools for this purpose import json import itertools import shapefile with shapefile.Reader("ne_10m_admin_0_countries.shp") as shp: # print(shp.fields) # Fields info # print(len(shp)) # To count to help decide range features = [] for s in itertools.islice(shp, 10000, 20000): # Loop only on a subset of ...


7

A generator is not subscriptable and iterRecords() returns a generator. Instead, use shapeRecords() (or records()). It gives you a list. rows = shapefile.Reader(shapefile_path).shapeRecords()[0:100] for row_num, row in enumerate(rows): print(row_num, row)


2

In Source code for cartopy.io.shapereader if reader.shp is None or reader.shx is None or reader.dbf is None: raise ValueError("Incomplete shapefile definition " "in '%s'." % filename) So one of these files is missing


1

You can use the join field tool to join your DBF to the raster attribute table, then use the tabulate area tool with shapefile / country as the in_zone_data / zone_field and raster / type as in_class_data / class_field


3

You have to specify the layer (feature class) name in the path string. layer_name = "countries" layer = QgsVectorLayer("path.gdb|layername=" + layer_name, "Countries", "ogr") QgsProject.instance().addMapLayer(layer)


3

Maybe this will help: library(sf) eco_l2 <- st_read("na_cec_eco_l2/NA_CEC_Eco_Level2.shp") Reading layer `NA_CEC_Eco_Level2' from data source `/home/micha/work/tmp/na_cec_eco_l2/NA_CEC_Eco_Level2.shp' using driver `ESRI Shapefile' Simple feature collection with 2261 features and 8 fields Geometry type: POLYGON Dimension: XY Bounding box: ...


0

Try this video: https://www.youtube.com/watch?v=N6WOzM00FLk Not mine but helped me a lot.


1

You are printing the NO statement inside the loop for every file that isn't a .shp. That should be left to the end, after you have checked. You are also printing the name of the file if it is a .shp, but that wasn't in your question. Get a list of files Filter the list by type (lower case in case they end with 'SHP') Output results files = os.listdir(r'E:/...


2

Move the check for no files out of the loop items = os.listdir(r'E:/folder/test') print('checking directory...') print('your files are: ') print() filenames = [] for names in items: if names.endswith(".shp"): filenames.append(names) print(filenames) if not filenames: print("there are no files ending with .shp in ...


9

Here my solution, entirely in QGIS (version >= 3.14, for versions >= 3.0 and < 3.14, maybe few adjustments are needed). The solution is to create a points grid in QGIS and retrieve data from JavaScript functions and update each point. The entry parameters of the function are : QGIS interface iface The solar eclipse elements, you can find it in the ...


3

It's quite easy if you have the same columns for all your shapefiles and GDAL/OGR installed and know how to use shell commands. In this case, you can use ogrmerge.py with -single option. Complete documentation with some examples is available at http://gdal.org/programs/ogrmerge.html


6

You can achieve this using ogr2ogr with # SHP from Natural Earth Data ogr2ogr out.shp ne_110m_admin_0_countries.shp -dialect SQLite -sql "SELECT * FROM ne_110m_admin_0_countries WHERE ST_intersects(ST_Buffer(ST_GeomFromText('POINT (1 43)', 4326), 10), ne_110m_admin_0_countries.geometry);" To sort out execution with index, you may want to add after ...


15

Demo: SOLUTION: This is a two-step solution, Step 1: Browser (JavaScript). It copies the popup text to the clipboard. Step 2: QGIS (Python). It gets the text from the clipboard, parses to obtain latitude and longitude, then adds the point to the active layer. In Browser: Open the website. Open DevTools Ctrl+Shift+C. Source > Snippets > New snippet &...


0

It turns out, when you create a a feature dataset within a file geodatabase, you have to define the coordinate system of the feature dataset to which you will import data at the time of the creation of the feature dataset (you cannot assign later), regardless of whether or not the data that you will import has a defined coordinate system.


3

There are two sets of classes for spatial data in R, and you are mixing them. They do not mix. The sp package provides the classes that start Spatial...., such as your SpatialPointsDataFrame. Then the sf package provides different classes for the same sort of data, and these are usually called "sf classes" or "spatial data frames". You ...


0

Please have a look at zonal statistics and zonal statistics as table tools. Given the objective, you might want to use Zonal Statistics as Table tool. It will generate a table as output (in contrast to the Zonal Statistics which returns a raster as output). Once you get the table, you may simply join it to your input shapefile to get population for each ...


0

okay, DatabaseUpdater Transformer is the answer of the question, but it my case it worked only with FME version 2021 and NOT 2020. In DatabaseUpdater you connect shapefile which in this case is updater, and for the matching columns you choose x and y and for the column to be updated you choose z.


3

Try this: import arcpy, os arcpy.env.workspace = r'C:\GIS\data\testdata\clip' #Your folder with the shapefiles to be clipped outfolder = r'C:\GIS\data\testdata\clipout' cutter = r'C:\GIS\data\testdata\clippolygon.shp' for fc in arcpy.ListFeatureClasses(): #List shapes in workspace print 'Clipping: ', fc arcpy.Clip_analysis(in_features=fc, ...


1

The FeatureWriter is a good transformer for this type of work. Here is an example: There is a shapefile of pois's and atms are cashpoints - so in this case used a tester to pull atms and renamed them. These are then written back to the original geodatabase but this time it has been set to update on just the 108 records changed by the id. You can control if ...


1

There are a few ways you can achieve this. Probably, what I would do is to add another reader to your workspace to read the file geodatabase data. Then, connect both inputs to GeometryExtractors. This will create an attribute which embeds the geometry. You can then use a FeatureMerger to merge the attributes of the shapefile to the file geodatabase feature ...


0

The issue got resolved. Since the database has layers from different sources, the CRS of one particular layer was not appropriate for web display. After reprojecting, the application is working perfectly now.


4

How about using zip: import geopandas as gpd import os frames = [gpd.read_file(r'C:\GIS\data\testdata\ak_riks.shp'), gpd.read_file(r"C:\GIS\data\testdata\ak_riks_2.shp")] names = ['r1.shp','r2.shp'] outfolder = r'C:\GIS\data\testdata' for frame, name in zip(frames, names): frame.to_file(os.path.join(outfolder, name))


2

A list comprehension is very Python. clipped = [gpd.clip(s, boundary) for s in shapefiles] But often it's easier and more flexible to loop, especially if there are clip polygons for each shapefile. Let's say clip_bounds contains a list of clip bounds in the same order as the shapefile list: shapefiles = glob.iglob('E:/folder/shapefiles/*.shp') for file, ...


1

Following code you can use to get what you want: for shapefiles in gdfs: clipped = gpd.clip(shapefiles, boundary) clipped.to_file("your local system path") Here I assume that your all shapefiles and boundary polygon are having same projections. If not then you can change the projection of your boundary polygon (if it's also in a ...


1

I'm seven years late, but Statistics Canada's postal code data (for FSAs and, if you use the PCCF, LDUs) are not wrong (generally). Postal codes in Canada can correspond to disjunct polygons, a bit like Russia proper and its territory of Kaliningrad. Both are Russian territory, but Kaliningrad Oblast (once known as East Prussia) is Russian soil. So goes it ...


4

I don't think it's possible. This is because shapefiles cannot store field names of more than 10 characters. (Shapefile is an old format from the 1990s) If you convert a data table into shapefile format, the field names are truncated (if longer than the limit). So, if you had longer field names, the information is already lost when it's converted to a ...


0

I encountered a similar issue and searching for a solution! How can I round the coordinates of the "SHAPE" column in an SDF created with pandas.DataFrame.spatial? The "round function would not work against "SHAPE" because it is not a number. import pandas as pd df = pd.DataFrame.spatial.from_featureclass("c:\\test.gdb\...


5

You can store each GeoDataFrame in a list or a generator to avoid any memory issue. import glob import geopandas as gpd files = glob.iglob('E:/folder/*.shp') gdfs = (gpd.read_file(file) for file in files) # generator # A list is an option for small files # gdfs = [gpd.read_file(file) for file in files] for gdf in gdfs: # clip stuffs Or import glob ...


0

I think you need to design your analysis to fit into memory and optimise so that you're only loading data from the Sentinel-5 data once. Here's a suggestion: Load all your polygons Use the find_datasets command in the ODC to identify the extents of each dataset Do an intersection with polygons and datasets, to group the polygons by dataset Load one dataset, ...


1

arcpy.env.workspace = r"C:\\" using a raw string, the interpreter will escape backslashes, so this will result in 4 backslashes. Remove the "r" as you have to escape to form a valid string. >>> dir = r"c:\\" >>> dir 'c:\\\\' here you are also prepending a backslash to the file name, so this will end up with 5 ...


0

I think you should look at this line of code: arcpy.env.workspace = r"C:\\" It does not look like a valid workspace so as a test try using: arcpy.env.workspace = r"C:\temp"


0

Assuming that the correct way to proceed in this case is simply to do mean <- extract(S2_stack, polygons2, fun=mean) as said in the comments, however, if you want to pass your stack in lapply you have to unstack it first. mean<-lapply(unstack(S2_stack), FUN=function (S2_stack) {data.frame(mean=extract(S2_stack, polygons2, fun=mean))})


2

You already have the geometry as a shapely multipolygon/list of polygons. Add it directly to the geodataframe. gdf3 = gpd.GeoDataFrame(geometry=[finalpol]) # Note GeoDataFrame geometry requires a list gdf3.to_file(filename='dPolygons.shp', driver='ESRI Shapefile') # Or gdf3 = gpd.GeoDataFrame(geometry=outmulti) # outmulti is already a list gdf3.to_file(...


3

The max limit for a shapefile it's 2GB because offsets in the files are specified with a 32bit integer, in other words, the limit is structural due to the file format. When downloading shape-zip larger than 2GB GeoServer should split the data into multiple shapefiles, there is logic ensuring that, but in your case, there might be an issue with how offsets ...


1

I've had a similar issue, the only difference was that I wanted the Sentinel-3 footprint that I extracted from the xfdumanifest.xml. Fixed it by including a srsName tag with a CRS that was GML3+ compliant to my GML file, I'll paste it here for you if it can be of any help: <gml:Polygon xmlns:gml="http://www.opengis.net/gml" srsName="http://...


Top 50 recent answers are included