texts: Choosing technology to start Web GIS development?

texts:
I am very new to GIS development, and to be be frank I have no background about it at all. I searched the web but the tutorials I found seemed to assume the reader has some background information.

the thing is that I am confused about what to read or learn, there seems to be lots of technologies, and I feel lost since some speak about openlayers, geoserver, mapserver, google maps, and open street maps.

So here is what I am supposed to develop, and I hove you could give me an advice about which technology to use, and where should I start reading - given that I know almost nothing -.

Case 1: a closed system for about 20 users only, who can specify locations on the map, and the web application will store the latitude and longitude of the locations and show the markers. I wanted to use google maps api, but I cancelled that since there license requires you to purchase the service if the system is a closed one. so what technology should I use in such case? I need a free option, also I will be only using web server, so if the solution includes using my own geoserver, or something like that I won't be able to do it.

Case 2: I am supposed to display the roads and routes between two given points, and probably add some notes on the map. For this I case I can use my own map server/geo server, but again I want your suggestions.

The solution needs to be open source.

Finally, I hope you could tell me what to start reading first.


texts: How to use the JS GEE API to allow users to download vector data from the user interface when pushing a button?

texts: I have created a GEE app where users can generate training data for a simple classification. I would like for them to be able to export their training data so that it can be reused. I have already tried to create a button that would use Export.table.toDrive. This approach works when I can enter my tasks tab and click run. However, this does not work in the UI.


texts: Dividing shapefile with many polygons into equal parts

texts:
I have a shapefile with hundreds of caribou movement corridors.
using this data, I am trying to create a map displaying the frequency of use to determine important movement and migration pathways for different herds.

The original shapefile has several years of data of individual collared caribou movements derived using Brownian Bridge Movement Modelling.
From this file, I have successfully created my desired outcome by:
1. Using the union tool
2. Multipart to singlepart
3. Spatial join to itself
4. Change shapefile symbology.

The problem I am now encountering is that I have extended my study dates and included several years into a single model. When I try to follow the above steps with the large dataset, I a thrown an error code 'insufficient memory' at the spatial join step (3.).
Clearly, my computer does not have enough ram to run this calculation, so I have been trying to divide my shapefile into 4 equal parts, using a fishnet grid, which I plan to merge back together after I run the spatial join on each parcel. I have successfully created the grid and again using the Union tool, created 4 separate shapefiles of data. Now, when I try to 'split' the caribou polygons from the grid attribute, the operation fails.

I understand that ArcGIS pro and perhaps a patch for arcmap could help my computer access more of its ram to run such operation, however, I do not currently have access to these, so am looking for some kind of alternative method of achieving my desired results.

I will try to add some images shortly.


texts: Calculating the percentage of polygons overlapping within the same layer in QGIS

texts: I have a layer with over 400 polygons. They are all really close within the same area of a city. I would like to calculate the percentage of overlap of every polygon with every other polygon. I unfortunatley only find information about how to do this with two different layers. But I want to calculate the overlap of the polygons with each other in only one layer. My goal is for example to know that Polygon1 overlaps 40% with Polygon2 within this layer. Maybe you have some ideas.


texts:
I'm working on exporting projects from QGIS to QField for data collection.  I have read the documentation that requires the CRS to be the same for the layers and the overall project.  The project is in WGS84-EPSG:4326 as are my layers.  The basemaps are in EPSG:3857 - WGS 84 / Pseudo-Mercator.  If I change this in the Layer Properties for the base map my project ends up in the ocean off the coast of Africa.  If I change the Project CRS to EPSG:3857 - WGS 84 / Pseudo-Mercator, the layers end up in some unknown place.
How can I do this?
I usually work in GA83-WF (USA-Georgia-West) in AutoCAD Map3D and then move in and out of QGIS.

Add your data and basemap and adjust symbology. See the Resources section of this document for QGIS tutorials.
a. PLEASE NOTE: You will want to make sure your data layers are saved with the same CRS as your project in order for them to show up properly in QField. If you want to use a basemap, it is best to set your project CRS to EPSG: 3857 (WGS 84/Pseudo-Mercator) since popular basemaps you can access through QuickMapServices in QGIS use this CRS.


texts: Filtering value layer according to condition of other layer using QGIS?

texts:
I have 2 point layer: layer A have 500 point; layer B have 50 point and these point have coordinate duplicate with some point in layer A.
I want to find these 50 point out of the 500 points in layer A.

Is there any way to filter ?

I don't see any way to filter this layer based on conditions other layer.


texts:
I'm using San Francisco OSM database, which contains the name of each edge, the id source and id target of each vertex.
I need to use the information of the most congested routes in San Francisco, but the problem that i didn't found the information that corresponds the data that i have on my data base.
For example:
Jackson Street to 16th Avenue Overcrossing is a congested road.
On my dataset:

All these vertex have the same name "Jackson Street" and I don't know which one is the vertex that compose the road to "16th Avenue"


texts: Functions disappear with layout view on map QGIS 3.16

texts:
I'm using QGIS 3.16 LTR and MacOS Big Sur.
When i'm opening a new print layout it is opening in a new tab and the original map is disappearing (turning black). I always need to close the print layout in order to change something in my layers.
Also some default shortcutes are not working in the print layout tab. It is not recognizing ESC for example.
Can anyone help me with that issue? I think it has something to do with the mac update but i can't find anything online. This issue started yesterday. Before everything worked fine.
I tried also an old version of QGIS 3.14 Pi but there I have the same problem.


texts: Distinction between pixel-based and object based classification?

texts:
I am struggling to clearly understand the distinction between pixel-based and object-based classification in the remote sensing domain and am hoping someone from this community can provide insight.

Based on the information I have so far, my current understanding is along these lines:

Pixel-based classification:
Classification is done on a per pixel level, using only the spectral information available for that individual pixel (i.e. values of pixels within the locality are ignored). In this sense each pixel would represent a training example for a classification algorithm, and this training example would be in the form of an n-dimensional vector, where n was the number of spectral bands in the image data. Accordingly the trained classification algorithm would output a class prediction for each individual pixel in an image.

Object-based classification:
Classification is done on a localized group of pixels, taking into account the spatial properties of each pixel as they relate to each other. In this sense a training example for a classification algorithm would consist of a group of pixels, and the trained classification algorithm would accordingly output a class prediction for pixels on a group basis. For a crude example, an image might be partitioned into n segments of equal size, and each segment would then be given a class (i.e. contains object / does not contain object).

Is this thinking accurate regarding the meaning of these terms, or is there something that I have missed?


texts: Is it possible to export a WMTS layer using ArcGIS Pro?

texts: I have a WMTS basemap layer I loaded via a server connection. Are there any ArcGIS-specific tools that allow me to export the images from that layer so I can work with them directly?


texts:
I'm using QGis 2.18.3 on Windows 10.

I'm trying to join a point layer with a polygon layer. However, since it is a city grid map, points that do not fall within the polygons do not get joined.

My expectation was that QGis could join those points to the nearest polygon.

Can somebody help me with that?

Here's a picture of my map.


texts: How to use pgRouting's dijkstra_sp_delta_directed with WHERE clause?

texts:
Right now I am using pgRouting for vehicle routing and I have problem with Dijkstra shortest path. My  query is:

SELECT * FROM shortest_path('SELECT gid  AS id,
      source::int4 AS source,
      target::int4 AS target,
      cost::float8 AS cost,
      reverse_cost::float8 AS reverse_cost
from network ',
257027,
276521,
true,
true)

But, my table network is very large, so I should use "dijkstra_sp_delta_directed" wrapper function with bounding box to make query faster.

Select gid from dijkstra_sp_delta_directed('network', 5700, 6733, 0.1,true,true);

But I don't know how to write some where clause like "from network where road_class=1"

Also, I don't know where and how to write (road_class * cost)::float8 AS cost


texts: Is it possible to disable Infowindow in cartodb editor, not in cartodb.js?

texts:
I am using cartodb Editor( to create maps, is it possible to disable the infowindow? I've found a few things for cartodb.js, but can't seem to figure out how to do this in the Editor on their website.

Or perhaps something like this? but within the editor itself?


texts: Cannot publish mxd file using ArcGIS publisher

texts:
I am having trouble publishing a map file using ArcGIS publisher. This mxd is intended for use by ArcReader but I keep getting a "file access error message asking me to check that the file permissions are correct and that it is not already open by another application".
I have been granted permission to publish on a shared drive and do not believe this mxd is being used by another application.
Does anybody have an idea what else I can try to make it work?


texts:
In Google Earth, the Rio Grande Gorge Bridge is not there. They've got the road going in and out of the canyon. This would be a blast to drive, but there really is a bridge there. They had to have purposely taken it out. Why?
Lat/long: 36.476176° -105.732763°

Edit: Same applies for:
Royal Gorge Bridge in Canon City, Colorado.
Lat/long: 38.461088°-105.325118°

New River Gorge Bridge in Victor, WV
38.066808°, -81.082553°


texts:
I recently upgraded, desktop, server, engine and SDE to 10.3, and of course had to upgrade some other features and extensions. It seems that now that I've upgraded to 10.3 my new 10.3 geocoders will not find or identify addresses. Strangely, when using the Address Inspector on the Geocoding Tool bar it will identify address with the 10.3 GC, but won't find any when using the search bar.

Has anyone had any similar issues?


texts: Esri.ArcGISRuntime.Toolkit - Map legend with more than 2 levels

texts:
I'm having trouble creating a multi-level legend for my map. I have a ArcGISDynamicMapServiceLayer that has 3 visible layers beneath it. The first layer then has 7 sub-layers beneath it. In my legend I want to show the following (I'll user my layers' actual names):

Contaminated Soil
-----> Former Cleanup Areas
       -----> Soil Remediation Complete

For the last layer Soil Remediation Complete I want to actually display the color coding associated with the layer.

This is what my code is so far:

<esri:MapView x:Name="MyMapView">
    <esri:Map x:Name="Map1" >
        <esri:ArcGISTiledMapServiceLayer ID="Base Layer" DisplayName="Base Layer" ServiceUri=https://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer  />
        <esri:ArcGISDynamicMapServiceLayer ID="Contaminated Areas" DisplayName="Contaminated Areas"  ServiceUri="https://xxxxxxxx" />
    </esri:Map>
</esri:MapView>

Can someone point me in the right direction to creating a legend for this. I have been surfing the web for 2 days now looking for something that can help but am coming up empty.


texts: Which is better, DEM or Triangles (TIN/Terrain)

texts:
Which is better DEM or TIN?

What are the benefits/negatives of each?


texts: Create a new BasicOverposterLayerProperties

texts:
ArcGIS 10.1 with VB.NET

I have been trying to figure this out for days but can't seem to make the connection. For a layer I am trying to create a label. I create a new  BasicOverposterLayerProperties

Dim pBasicOverposterLayerProps As IBasicOverposterLayerProperties4 = New BasicOverposterLayerProperties
        <setup properties>
        Dim pAnnoLayerProps As IAnnotateLayerProperties = Nothing
        Dim pAnnoLayerPropsColl As IAnnotateLayerPropertiesCollection = pGeoFeatLayer.AnnotationProperties
        pAnnoLayerPropsColl.QueryItem(0, pAnnoLayerProps, Nothing, Nothing)
        Dim pLabelEngineLayerProps As ILabelEngineLayerProperties2 = CType(pAnnoLayerProps, ILabelEngineLayerProperties2)

pLabelEngineLayerProps.BasicOverposterLayerProperties = CType(pBasicOverposterLayerProps, IBasicOverposterLayerProperties)

Setup all my label properties and apply it and it works great for the ESRI Standard Labeling Engine. Now with the widespread use of Maplex I want to convert the labeling properties before applying.

Dim pUidOld As UID = New UID
Dim pUidNew As UID = New UID
pUidOld.Value = "{01004145-0D1C-11D2-A26F-080009B6F22B}" 'ESRI Standard Label Engine
pUidNew.Value = "{20664808-0D1C-11D2-A26F-080009B6F22B}" 'ESRI Maplex Label Engine
Dim pAnnotationPropertiesConverter As MaplexAnnotationPropertiesConverter = New MaplexAnnotationPropertiesConverter

Dim pGeomType As esriGeometryType = pFeatCls.ShapeType
Dim pDisplayTrans As IDisplayTransformation = MyDoc.ActiveView.ScreenDisplay.DisplayTransformation
Dim pAnnotateLayerPropertiesCollectionOld As IAnnotateLayerPropertiesCollection = ???

Dim pAnnotateLayerPropertiesCollectionNew As IAnnotateLayerPropertiesCollection = pAnnotationPropertiesConverter.Convert(pGeomType, pDisplayTrans, pUidOld, pAnnotateLayerPropertiesCollectionOld, pUidNew)

I can't figure out how to get the pAnnotateLayerPropertiesCollectionOld object from the pBasicOverposterLayerProps.


texts: Creating points along a line and keep the Z value

texts:
I want to create points along many line features and keep the Z value from the line in the points. I'm trying to do this using Densify and Feature Vertices To Points tool, but this takes a lot of time (almost 1.5 hours for 600 features).

Is there any method that I can use that is faster than the one I'm using today?

I'm working in ArcMap 10.6 with Advanced license.


texts: Difference between osgeo gdal/ogr python module and GDAL C++ library

texts:
I am developing application in python and I use the osgeo gdal/ogr python bindings for gdal, but its not clear to me what exactly is the difference between the gdal/ogr in python and the gdal library in c++.

What I would like to know is, if I upgrade to the last version of c++ gdal (1.12), does my python applications using osgeo gdal/ogr will automatically use this new library or do I have to upgrade osgeo as well?

As fas as I know, gdal c++ library is being upgraded regularly but not the osgeo gdal/ogr python library.


texts: QGIS Model Builder Script - Possible to merge layers in a folder?

texts:
I would like the user to select a folder via the interface in which all the vector shapefiles inside that folder would be merged.

I have been using the Model Builder and looking at example scripts and tried to incorporate simple commands into my own script. The following code is incorrect but any guidance would be appreciative:

##Test=name
##input=folder
##merge=output vector

outputs_0=processing.runalg("qgis:mergevectorlayers", input, None)

I have looked at Merge scripts from the Processing and SAGA plugins, but I was wondering if it's possible to do something similar in just one or two lines of code?


texts: High resolution is RGB image is not coming clear and sharp?

texts:
I am woking with LISS IV image of Indian Satellite Resousesat-2. This image has 3 bands i.e., NIR, Red, and Green band and spatial resolution of five meter.

When I view the image in ArcGIS the details are not sharp as it should be. For example, below is the screenshot of the image I am using i.e., five meter spatial resolution, and a screenshot from a published article which uses similar image at Eight meter spatial resolution.

5 m spatial resolution

8 m spatial resolution (screenshot from a publish article)

Why does my image not have the fine details despite having the finer resolution than later? Why the colors are blurry for different land cover?

How to present it like the published image shown here?


texts: Exporting GeoTIffs from Data Driven Pages for use on Tablet?

texts:
I'm a beginner at ArcGIS 10.2 and have encountered a small problem that I couldn't solve by googling. I'm trying to export a set of raster and vector data for use on a tablet. The aim is to find LiDAR features in the field using the GPS from the tablet and I want to be able to switch between various visualizations and the polygons and lines of the interpretation.

For this purpose I try to batch create GeoTiffs with various visualizations of the data and export them using Data Driven Pages, since it lets you very conveniently tile a map and I wanted to create a hard copy mapbook of the data anyway. I found an ArcPy script that batch exports Data Driven Pages to tiffs and this is working fine, with the major problem that it's not possible to include spatial information when exporting from a Layout. Since I'm generating a very lage amount of tiles (over 600) and would like to do so with various visualizations, manually georeferencing them is not an option.


texts: How to select points that are within a polygon when both layers are FeatureLayers using a Geometry Service?

texts:
I have a map service that has 2 layers. One is a polygon layer that contains a single large polygon, and the other is a point layer that contains many points located both inside and outside of the large polygon.

I've turned both layers into FeatureLayers:

var polygonFeatureLayer= new FeatureLayer(myPolygonLayer, {
    outFields: ["*"]
});

var pointFeatureLayer= new FeatureLayer(myPointLayer, {
    outFields: ["*"]
});

And I want to use a Geometry Service to find out how many points are located within the polygon.

geometryService = new GeometryService(myGeometryService);

var relationParams = new RelationParameters();
relationParams.geometries1 = myPolygonGeometry;
relationParams.geometries2 = myPointGeometries;
relationParams.relation = RelationParameters.SPATIAL_REL_WITHIN;

geometryService.relation(relationParams).then(addRelateResultsToMap);

But I can't figure out how to get the geometry of the features from the FeatureLayer.

I thought the most obiously way was to use the FeatureLayers's graphics property which supposedly returns an array of Graphics for each feature in the layer.

I then figured after that I'd use graphicsUtils method getGeometries() which supposedly returns an array of Geometries for each graphic in the array.

Hey! That sounds great! Exactly what I need, but I can't get it to work.

I did a simple test:

var polygonFeatureLayer= new FeatureLayer(myPolygonLayer, {
    outFields: ["*"]
});
map.addLayers([polygonFeatureLayer]);

polygonFeatureLayer.on("load", getPolygonGeometry);

function getPolygonGeometry(evt) {

   console.log(graphicsUtils.getGeometries(polygonFeatureLayer.graphics));

    }
}

The console just logs an empty array [].

In my mind this should work. Any idea why it doesn't? What do I need to do to get this to work?

EDIT

If you have a different approach, then I'll consider it, but it has to work knowing that both layers are part of a map service.


texts: Batch importing same layer into multiple mxd files using ArcPy?

texts:
I am new to python.

Will someone help me with the code to import the same layer into all the arcmap files in a folder on my desktop?

I used the code below to import the layer into one map but I need to import it into 20 maps.

import arcpy.mapping as mapping

mxd = mapping.MapDocument("CURRENT")
df = mapping.ListDataFrames(mxd)[0]
layer= mapping.Layer(r"c: PATH TO LAYER FILE")
mapping.AddLayer(df, layer, "TOP")


texts: Add attribute information from polyline to a point shapefile (QGIS or ArcGIS)

texts:
Is it possible to add attribute information from a polyline to a point dataset. I have a network layer with information per street how many scooters went over this street. I also have a point dataset of people with information on the level of noise annoyance. I want to add the number of scooters passed by to be added in the attribute table of the point layer. Not all points are exactly on the lines of the streets, so maybe there is a functions which links the closest line with its attributes to the point layer.

Is this possible and if so, how do I manage this?

I have ArcGIS and QGIS software installed.


texts: QGIS - print composer legend default settings

texts:
In the QGIS user guide, it says

The legend will be updated automatically if checkbox Auto-update is checked. When Auto-update is unchecked this will give you more control over the legend items.

Can I change any settings so that Auto-update is not checked when the legend is created? And can I build the legend from scratch by default, rather than dismantle it by having to remove unwanted layers (many of which may be invisible on the map canvas).

The reason is: when I have a layout that includes raster layers with bands (in my case Ordnance Survey tif), the legend takes a long time to load when first created then unchecking the box takes ages, then deleting the unwanted raster layers... If I have a lot of these layers, it can take over 5 minutes just to load the legend, so I can be looking at nearly 10 minutes just to get to the point that I can edit the legend.

I'm currently using QGIS 3.8.2

The similar question has not been answered in a way that solves the question I have asked, although it is a loose work around:
QGIS- how to prevent layers from appearing in legend on composer


texts: How to create a data package for arcgis 9.3 in arcgis 10?

texts:
I need to create a data package in ArcGIS 10 for a company that is using ArcGIS 9.3.
Is there any way to downgrade the versions like GDBs. I'm grateful if someone has a solution for this and help me.

Thanks.


texts: Is this transformation using square metres or kilometres (3857/32633)?

texts:
I am using geopandas to compute population per km2. However, I don't know if this transformation is using meters or kilometers.

gacs_df = gpd.GeoDataFrame(acs_df.merge(census_gdf2, on="geo_12", how="inner").drop_duplicates(subset='geo_12'), crs={'init': 'epsg:4326'})
gacs_df = gacs_df.to_crs({'init': 'epsg:3857'})
gacs_df["area"] = gacs_df['geometry'].area/ 10**6
gacs_df= gacs_df.to_crs({'init': 'epsg:32633'})
gacs_df["area"] = gacs_df['geometry'].area/ 10**6

I am replicating the code from Getting polygon areas using geopandas?. Thanks!!


texts: Any way to speed up an append of data to SDE database while it is still versioned?

texts: My group works on a database and about once a month I receive data that must be appended. This process has required me to reconcile and post all versions to default. I then pull down the database locally and delete all relationship classes and geometric network so I can append to the new data. Then I re-upload database, register data, recreate versions and grant privileges. Is there any way around this process?


texts: Can I plot points in ArcGIS if I only have distances from a set coordinate?

texts:
I understand there is a way to plot lat,long data from a table (dbf) in ArcMap, but is there any way to plot x,y data if you only have distances away from a set coordinate instead of lat, long?

To clarify and specify my task: Square plots were set up and a GPS coordinate was recorded for each corner. Then the distances (m) from the north corner were recorded in the x and y direction for every tree in the plot. I was hoping to plot these trees in ArcMap, but don't want to plot each tree individually by way of creating a shapefile and adding each feature in (too many trees!!).So, I was wondering if there was any way to upload all these points with one or a few steps?  I already have a .csv file with all of the tree numbers/attributes and their respective x and y distances away from the corners.


texts: Determining if point is surrounded using raster processing

texts:
I am attempting to improve a currently extremely cumbersome vector/python process for a natural hazard model. At the moment we have a lengthy script which generates distance/bearing lines from a given point to determine:

the type of polygon that it intersects (e.g.. Forest, grass, swamp, etc)
the distance to that polygon
how many of these lines intersect polygons, to determine how 'surrounded' it is.

There's a lot more involved but that's the gist of it. I am trying to find a way to improve this and am currently stumped on part 3. The idea is to determine if a point is completely surrounded by polygons, within say 200m

So in my attached image, I would want point A to be marked as being at higher risk than point B as it is completely surrounded by my polygons. This is repeated for about 13 million points so is not a small task and I would rather have a surface to derive values from, rather than running our script. I'm thinking there's got to be a variation of hydrology tools or cost-path to do this but I can't seem to get my head around it.

How could I go about this?


texts:
I need to clip all of the feature classes in the .gdb to a layer then add the results to a results dataset. However it returns the error.

  Error 000210:Cannot create output C:...gdb\FC_Results\tranterm Failed
  to execute.

How do I set the output parameter of clip to create a featurename_clip output for each feature?

import arcpy
import os
arcpy.env.workspace = "C:/Student/WashingtonDC.gdb"
out_workspace = "C:/Student/WashingtonDC.gdb/FC_Results"
ClipFC = "C:/Student/WashingtonDC.gdb/citylim"
arcpy.env.overwriteOutput = True
ListFC = arcpy.ListFeatureClasses()
for FC in ListFC:
    output = os.path.join(out_workspace, FC)
    arcpy.Clip_analysis(FC,ClipFC,output,"",)


texts: I am looking to create a customized version of the find tool where the layer and field tools are pre-configured.  Ideally I could add this to ArcMap 10.1 as a button for end users to who use same search using the find button repeatedly.  Any help would be appreciated.


texts: Layer is not valid error when connecting QGIS to PostGIS

texts:
I am new to QGIS.
I am trying to connect a view from PostGIS to QGIS (3.16.4, same with 3.12).
The PostGIS connection in QGIS browser have no issues and when "Test Connection" it says it can connect.
But when I am  trying to add the layer I am getting this error:

Layer is not valid: The layer dbname='dwa' host=localhost port=40802 user='user' password='user' key='"bin"' srid=4326 type=Polygon table="public"."height_grid_layer_view" (bin) is not a valid layer and can not be added to the map. Reason:

The pgAdmin 4 shows the view in Geometry Viewer without any problem with the SQL:
SELECT * FROM public.height_grid_layer_view
What can I do?


texts:
I download data from http://land.copernicus.eu/global/products/lai?qt-lai_characteristics=5  and when I visualize it in ArcGIS I set epsg 3857, but scale is wrong and measure show pixel with resolution 1x1m.

Can you tell me the correct system for this data?


texts: Displaying elevation profile and storing DEM elevation back to KML line using QGIS?

texts:
When I draw a track in Google Earth, for example along a road, I can right-click it to "Show elevation profile", and it displays a 2D plot of distance vs elevation. However, if I save the KML, the elevation coordinate is always zero.

So let's suppose I import the KML LineString in QGIS, and I have a SRTM DEM layer, how can I:

Display an elevation profile (distance along road vs DEM elevation at point);
"Replace" the (zeroed) elevation profiles of KML by DEM elevation values, so that I could export back to kml and have elevation "hardcoded" in the LineString.


texts: How to create a squared buffer around a point with Python GDAL?

texts:
I have a large collection of points with X and Y coordinates in a projected coordinate reference system (EPSG:28992). I would like to compute a square buffer around it (say, 500 meters around the point), so I can crop this polygon with a raster layer (in TIFF format) and proceed with a further analysis with the values within the cropped windows.

I am aware that this is possible in Python for QGIS, but for the sake of the workflow I am implementing, I prefer to have it entirely in the Python GDAL environment, so it becomes automatic. Since I know the X/Y coordinates and the size of the grid cells in the raster, I could add/subtract to the coordinates the desired buffer, but I wonder whether there is a less handcrafted version of this "square-buffer" operation, which is less prone to errors. Other Python libraries (e.g. shapely) would be fine for me as well.

Note that a band.GetRasterBand(1).ReadAsArray(X, Y, 5, 5) does not return the window around the point, but the window leaving the point in a corner.

For instance, the following code:

path_curr_lyr = r"path_to_tif_file"
tif = gdal.Open(path_curr_lyr)
dat = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 1, 1)
print(dat)
sq_buf = tif.GetRasterBand(1).ReadAsArray(9550, 9620, 5, 5)
print(sq_buf)

Accesses the raster grid cell X/Y (9550, 9620) and returning the value of forest coverage in that pixel, and the value of a window of width= 5 pixels.

>>> [[44]]
>>> [[ 44  87  99  86  43]
     [ 96  99  99  79  14]
     [ 99  99  98  34   9]
     [ 99  86  27  93  71]
     [ 98  98 100  97  65]]

As you can see, this approach returns a window in which the "44" is not a central value, but one in the corner. I would like to have a procedure in which a window around the given grid cell coordinate is created.

Any advice on how to implement this?


texts: Styling in OpenLayers2: unknown symbol name

texts:
I'm using Geoserver 2.8.1 and OpenLayers 2 to visualize cartographies, the geoserver documentation indicates I can use extra symbols in my SLD, the "shape://oarrow" for example :

 <Mark>
   <WellKnownName>shape://oarrow</WellKnownName>
     <Fill>
       <CssParameter name="fill">#000000</CssParameter>
     </Fill>
 </Mark>

But when I try to load this kind of style, I get this error from OpenLayers.js:

  Uncaught Error: shape://oarrow is not a valid symbol name

What am I missing here ?


texts: Google Map iframe KML not displaying correctly

texts:
Looking for help in why our embedded KML Google map stopped working today. This is the current code:

iframe width="100%" height="350" frameborder="0" scrolling="no" marginheight="0" marginwidth="0" src="http://maps.google.com/maps?q=docs:%2F%2F0B0CyyOTyQR_7T0NfQmhVakFndG8&t=p&output=svembed">

Here is a link to the Google Map with the KML that the site is pulling from:

https://docs.google.com/file/d/0B7SRdsTwGw_MTVJJMDVnZ3lhX2c/edit


texts:
To test something I have a very simple startup.py file in my .qgis/python folder:

import qgis
qgis.utils.iface.actionShowPythonDialog().trigger()

print "Welcome chump"

Python console launches fine, Welcome chump consistently prints out twice. Why is this?


texts:
I have tried without any success to install LASzip¹ ² (paper by Isenburg about LASzip) on Ubuntu 16.04.
The basic way gave me nothing relevant:
./configure && make & checkinstall

which seems to install some .a, .la, .so.6 and other .hpp files under /usr/local but no binaries (which are mandatory to install other "packages" like pdal)
I wonder if these tools were included in libLAS but it doesn't seem so: some of them does exist within libLAS and others not.


texts:
At the moment I am creating a webmap with PostGIS and GeoServer. I want to style the data using SLD-files, because my data keeps updating.

I have made SLD-files for Points and Polygons and that worked.

(I imported the data in QGIS, then styled the data and exported the SLD-file, which I then imported as a style in GeoServer. When I changed the data in the database, then (for example) the color changed of the polygon in the WMS-layer created by GeoServer.)

Now I did the same with a heatmap. I imported the data in QGIS, made a heatmap and tried to export a SLD-file. Unfortunately, this does not work. The file contains the following sentence:

FeatureRendererV2 heatmapRenderer not implemented yet-->

I cannot use this SLD-file now for GeoServer.

What do I need to do to make this work?


texts: What is the best or recommended way to store a GeoJSON FeatureCollection in MongoDB?

texts:
The comments in this question: Shapefile to GeoJson conversion using C# raised another one: what is the best way to store a GeoJSON FeatureCollection in MongoDB?  Is there any sort of convention when it comes to adding GeoJSON FeatureCollection to a MongoDB?

When I create a GeoJSON output with ogr2ogr:

ogr2ogr -f GeoJSON outshp.json c:\temp\us_counties_.shp

I get a json file and load that into MongoDB collection.  The result is collection with a single document of type FeatureCollection, which contains all of the features:

It seems like a mongo collection and a featurecollection are redundant.  Is there any sort of convention on how the Feature documents are stored?  Are they typically left in their GeoJSON FeatureCollection? Or are they typically moved out of their GeoJSON FeatureCollection container once in a MongoDB?

Would there be any advantage in doing so?

The scope of my questions is the group of features that would be analogous to a shapefile or featureclass, not for cases were you start nesting geometries or feature json inside other documents.


texts:
I have been struggling for quite a long time now with keeping my project files organized.

What are your tips for keeping your datasets, images, shapefiles, etc.. organized?


texts: How to filter out a geometry from another geometry in Google Earth Engine?

texts:
I have a shapefile out of which I need to filter out another shapefile. The .filterBounds() function filters to the second shapefile, which is actually what I need to "delete". How do I do this? Here's the code I'm working on:

// sampling pixels for control group: (MA) - (PES)

var MApes = MA.filterBounds(PES);

var MAnot = MApes.not(); // this was a solution I found in another thread, but GEE didn't recognize the function

// lossyear Atlantic Forest without the PES area

var lossyearMA = ee.Image("UMD/hansen/global_forest_change_2020_v1_8")
  .select(['lossyear'])
  .clip(MAnot); //


texts: Which geo-format does this data represent


texts:
How to add a dynamic scale bar in (or close to) map extent while data view is turned on?

I'm panning a lot (raster images) and I need to get more precise "feel" of distances just in one glance since I'm panning in various scales.


texts: Clipping raster in Python with ArcPy and GDAL using coordinates

texts:
I am having trouble clipping a raster with GDAL.  I have the coordinates of the bounding box I would like to clip to. There is no error but also no output file is being saved. Could it have to do with my coordinate system? My coordinates are in a UTM zone (same as raster I am clipping)

import arcpy
from arcpy import env
from arcpy.sa import *
import os
import subprocess
import gdal

path = "R:\\Sentinel-1\\kevf\\"
image = "S1A_IW_GRDH_1SDH_20190303T100649_20190303T100714_026172_02EBE6_AD79_out"
segmentedRasterSave = "R:\\Sentinel-1\\kevf\\seg\\S1A_IW_GRDH_1SDH_20190303T100649_20190303T100714_026172_02EBE6_AD79_out_seg.tif"

def split(segmentedRasterSave, image):
    arcpy.CheckOutExtension("Spatial")

    elevRaster = arcpy.sa.Raster(segmentedRasterSave)
    Extent = elevRaster.extent
    #set values
    xL = float(Extent.XMin)
    xR = float(Extent.XMax)

    yB = float(Extent.YMin)
    yT = float(Extent.YMax)

    x_delta = xL - xR
    x1 = xL + float((x_delta/3))
    x2 = xR - float((x_delta/3))

    leftSave = path + "split\\" + "_left.tif"

    ds = gdal.Translate(leftSave, gdal.Open(segmentedRasterSave), projWin = [xL, yB, x1, yT])
    ds = None

split(segmentedRasterSave, image)

There is no error message, just no output file is being saved.


texts: Grass processing not configured in OSGEO4w QGIS 2.14.1

texts:
I get the below error when I try to use grass processing from within QGIS. The error is given on all commands and I have got this error on 2 machines. I am using windows.

  "This algorithm cannot be run :-(
  It seems that GRASS is not correctly installed and configured in your system. >Please install it before running GRASS algorithms."

I have tried downloading standalone GRASS 6.4.4 and changing the grass and mysys folder to that but no sucess.

Any help?


texts: Using Geodatabase Annotation in ArcView 3.2?

texts:
I've been asked to take a parcel layer (shapefile) that contains house numbers and create an annotation layer so the house numbers are displayed in a separate annotation layer rather than a label to the parcel layer.

They'd like the ability to move and rotate the house numbers for readability.

Is there a way to create an annotation layer for ArcView 3.2 that will do this?

I've experimented with exporting a couple different ways and it loses it's usability once it comes outside the geodatabase in which I originally created it.

Will ArcView 3.2 even read a geodatabase?


texts: Fill larger polygon empty space with closest subset of polygons

texts:
Looking for an ArcGIS or QGIS solution.

I need to divide up the space of the large single green polygon. Each new polygon would  represent "nearness" to the purple polygons and have the ID of the purple Polygon it is closest to. Need to fill all the empty green space.

This would be a Thiessen Polygons solution if the purple were point but my input is a polygon and I need the equivalent but "what polygon am I closest to" and the ID of the original polygon has to carry over to the new closest polygons. I am open to turning the purple polygons into raster and having a raster solution.

I know I can turn the polygon boundary into equally spaced point and do this but I am looking for something more elegant and fully accurate.


texts:
I am a beginner with QGIS. I am trying to learn it through lynda.com. I have tried downloading the nanaimo crime map in GeoJSON format, but I get a message saying that the data source is invalid.

I am bemused at to what I do as I am following the instructions to the letter.


texts: Publish an orthophoto without needing mxd file

texts: I wonder if anyone has a tool or a method that can help me to Publish an orthophoto without the need to include it in mxd file.


texts: Italy, road network shapefile with full speed information

texts:
I am perusing the internet since one day and still cannot find something satisfying.
The most important thing I am interested in is the speed limit (or an average speed) of the respective elements. I don't need information on small residential streets or city-level streets. Just the main road network would be enough.
The extracted OSM stuff I downloaded from geofabrik.de or OSM extract Italy has tons of missing speed values and is thus not useful for me.

The closest I could get was the NASA SEDAC Global Roads Open Access Data. I subsetted on Italy via a GADM0 shapefile, it's slim, only with main roads, but there is no speed information.
I found this page with information on speed limits, but a manual editing would be virtually impossible.

Even a file that has consistent and complete (in the files I have, half of the roads go as "unclassified" - even though some of them are highways) information on the category of the street would be kind of o.k. In this case I could just impute average values for each type of street; which would be enough for the type of rough estimates of travel times I need.

Does anyone have a shapefile or information that could potentially help me?


texts:
This question has been bugging me for a while now, and I still haven't managed to solve it, so I'm giving it another go, hopefully someone may help me out.

I've been having problems loading several QGIS plugins. I've just done a reinstall of QGIS 1.7.4., followed all the instructions, so just did a default installation (I've tried the osgeo4w installer already and it only made matters worse...).

Whenever I try to use, let's say the "1 band raster color table" plugin, I get error messages. First of all, when I load the plugin, I get the message "Cannot continue. Problems!". If I close this message, the plugin stays open anyway and let's me select the options I want. When I click "Create", I get yet another error message, that starts off with this:

"Traceback (most recent call last):
  File "C:/Users/Tiago/.qgis/python/plugins\bcccoltbl1\bcctbl.py", line 826, in applyTBL
    myStats = self.myGRD.getStats()
AttributeError: 'str' object has no attribute 'getStats'"

Now, I've been told that this might have to do with me having non-Latin characters in the path... I've been trying to understand this, but I'm getting really confused... When I used QGIS 1.6, I had no problems of this sort. But I didn't change a thing in my disk's organization when I upgraded to QGIS 1.7, so how can this be, that the plugins ran smoothly before but not anymore? If I had non-latin characters, they would be there before, anyway, so the plugins should have had these same errors...

This is getting me really frustrated, because it will make it very hard for me to keep using QGIS if I can't even change the looks of rasters!

Please, I really need help with this. If necessary, I can send some more information.

Thanx!!!


texts: PostGIS Join Queries results in Missing MultiPolygons

texts:
I'm trying to join two PostGIS table, one of DMAs and another of states.

After getting some help, I ran this query:

CREATE TABLE combined AS
     SELECT nextval('polyseq_1') AS id,
     b.name AS state_name,
a.dma_1 AS dma_1,
CASE WHEN ST_Within(a.geom, b.geom)
    THEN a.geom
     ELSE ST_MULTI(ST_INTERSECTION(a.geom,b.geom)) end AS geom
FROM us_states b
JOIN us_dma a ON ST_INTERSECTS(a.geom, b.geom)

However, the results I get have some gaps, for instance, it looks like the Dallas DMA in Texas is completely missing. I'm not familiar with PostGIS but I am unable to figure why the query will cause polygons to be missing. A picture is below:


texts: Whenever I open the one project I've ever done in QGIS, I get a warning at the top of the screen about "Used a ballpark transform from EPSG:4326 to EPSG:26916". I presume this means I've done something wrong, I'm thinking that at least one layer doesn't have the correct projection (should be NAD83 UTM 16N). However, under Project>Properties, it lists the correct projection. I've also selected all of the layers, Set CRS to the correct projection and still I get this warning when I open QGIS.


texts: Python toolbox documentation from within pyt file

texts:
Python toolboxes are a great way to create custom, GUI toolboxes for ArcGIS without actually interacting with the GUI (since everything can be defined within the pyt file). Regrettably (and incomprehensibly) documenting a python toolbox does require the GUI (see Documenting a tool in a Python toolbox). This help page also states:

  For Python toolboxes, the documentation for the toolbox and tools are stored in .xml files that are associated with the toolbox and tools by name.

In a question on this topic the user @dmahr stated in his answer:

  Fortunately, Python supports for reading and writing XML, so you should be able to dynamically edit the help text from a separate script.

How would something like this be achieved? Since my toolbox is constantly changing, I'd  like code and documentation to live in the same file.

EDIT:

As @user2856 commented, the first step would be to look at the xml file generated when documenting the tool via the GUI. This step should have been part of the question, so here's the xml for the tool named "Tool" (saved with the name New Python Toolbox.Tool.pyt.xml:

<?xml version="1.0"?>
  -<metadata xml:lang="de">
    -<Esri>
      <CreaDate>20190401</CreaDate>
      <CreaTime>10300800</CreaTime>
      <ArcGISFormat>1.0</ArcGISFormat>
      <SyncOnce>TRUE</SyncOnce>
      <ModDate>20190401</ModDate>
      <ModTime>14352000</ModTime>
      -<scaleRange>
      <minScale>150000000</minScale>
      <maxScale>5000</maxScale>
      </scaleRange>
    </Esri>
    -<tool xmlns="" toolboxalias="" displayname="Tool" name="Tool">
      <arcToolboxHelpPath>c:\program files\arcgis\pro\Resources\Help\gp</arcToolboxHelpPath>
      -<parameters>
      -<param displayname="Input Features" name="in_features" expression="in_features" datatype="Feature Layer" direction="Input" type="Required">
          <dialogReference>

this is the in_features dialog explenation

</dialogReference>
          <pythonReference>

and here is the scripting explenation

</pythonReference>
        </param>
      </parameters>

This is a nice summary

      <usage>

This is how this tool should be used

</usage>
    </tool>
    -<dataIdInfo>
      -<idCitation>
          <resTitle>Tool</resTitle>
        </idCitation>
        -<searchKeys>
          <keyword>Here are some tags</keyword>
        </searchKeys>
    </dataIdInfo>
    -<distInfo>
      -<distributor>
        -<distorFormat>
          <formatName>ArcToolbox Tool</formatName>
        </distorFormat>
      </distributor>
    </distInfo>
    -<mdHrLv>
      <ScopeCd value="005"/>
    </mdHrLv>
    <mdDateSt Sync="TRUE">20190401</mdDateSt>
  </metadata>


texts: How to adjust the granularity of a heat map in QGIS?

texts:
Not sure I am using the correct terminology here, so I have tried to describe my problem in detail. Please edit as needed.

I have a layer of points, a layer of roads and a heat map in a raster layer.

I would like the heat map to reflect density of points in the point layer in such a way that I can zoom into a part of the city and see the colours for that area.

To make the heat maps I am using the heat map plugin in QGIS 1.8.

Unfortunately, my heat maps' resolutions don't seem to capture the level of granularity I am looking for. I have been fiddling with the settings but don't really know what settings are relevant in this case.

Here are my points with the heat map overlaid:

Here is what I get when I zoom in:

As you can see, then replacing points with roads, the result is not very illustrative:

I would like to obtain a result which gives a better idea of the point density at the neighbourhood level.

These are my current settings:

Update to respond to Underdark's comment below:

I have tried lower values in the radius field like you suggest but that results in two problems:

maps are mostly "cold", it's almost as if the points don't pass a certain threshold to really register. I've thought increasing the radius makes it more likely QGIS considers them as forming clusters.
As you can see in the picture below, the result is highly pixelated when one zooms in. Far too much for what I am trying to do.


texts: How to control how points are displayed across different programs

texts:
I'm completely new to GIS, but I'm trying to control how Points are displayed (shape, color, size, etc) across different programs (ArcGis, qGis, etc).

Is there a way to control the rendering of Points in a shapefile consistently in different GIS programs?


texts: Problem with Swiss coordinate system overlay in WMS

texts:
I have administrative boundaries data for communities in Switzerland from SwissTopo. Original data came in shape files with 'Bessel 1841 Hotine Oblique Mercator Azimuth Natural Origin' projection.

Full info from ArcGIS:

Projection: Hotine_Oblique_Mercator_Azimuth_Natural_Origin
False_Easting: -9419820.590700
False_Northing: 200000.000000
Scale_Factor: 1.000000
Azimuth: 90.000000
Longitude_Of_Center: 7.439583
Latitude_Of_Center: 46.952406
Linear Unit: Meter (1.000000)

Geographic Coordinate System: GCS_Bessel_1841
Angular Unit: Degree (0.017453292519943299)
Prime Meridian: Greenwich (0.000000000000000000)
Datum: D_Bessel_1841
  Spheroid: Bessel_1841
    Semimajor Axis: 6377397.155000000300000000
    Semiminor Axis: 6356078.962818188600000000
    Inverse Flattening: 299.152812799999990000

I tried defining and projecting this data into 'CH1903 LV03' (as described here) in ArcGIS 9.3 and loaded into GIS Cloud service.

When I overlay my data onto Google Maps or Open Street Map, they seem to be shifted ~200/300m NE from correct position. This happens to both original data in CH1903 LV03 and data reprojected to oter CS.

Has anyone of you encountered solution for such problem?


texts: TableToTable maintaining table header in csv file as alias in ArcGIS

texts:
I have large excel tables with geochemical information that I am importing to ArcGIS. TableToTable works, but there are likely restrictions to the number of columns that can be imported directly from excel as the last columns of my sheet are not included in the table when I open it in ArcGIS. Converting the excel table to csv and importing fixes this, but....

My question is related to the resulting field names and aliases. The headers in my excel/csv table contain special characters that are converted to underscore on import, but I would like to keep these headers as aliases for each field. It is possible to set the resulting aliases manually using the Field Map in the TableToTable conversion tool, but since I have hundreds of columns, I would like to automate this process in a python script.

Can I add the Field Map function to my script such that the headers in the csv table are maintained as Aliases in the new table?

import arcpy
inTable = "Data.csv"
outDB = "Data.gdb"
outTable = "Data"
arcpy.TableToTable_conversion(inTable, outDB, outTable)


texts: Is there a way in QGIS to display vector files on a shared drive without others being able to edit the file?

texts:
So just like lyr-files in ArcGIS.
I have read that qlr-files are QGIS answer to lyr. But as far as I can see when I try, the receiver can edit in the main source of the data through the qlr-file.
So basically I would like to display some files on a shared drive, where I am the only one who can manage the source in postgreSQL, geopackage or just shape.


texts: Including multiple layouts in one ArcGIS Desktop project

texts:
Is it possible to have multiple different map layouts (like the one accessed by clicking the little page icon at the bottom of the viewing panel of ArcMap) within one ArcGIS Desktop project?

I often find myself wanting to produce a number of different map outputs from one project, and this requires doing each one separately, and then 'undoing' it when I want to create the next one, making it difficult to go back and edit any earlier ones.


texts: Is there any difference between a stream's Longitudinal Profile and its Elevation Profile?

texts: I have derived elevation profile of a river channel using USGS SRTM 1 ArcSec. I am interested to know whether Elevation Profile and Longitudinal Profile are synonymous?


texts: Trying to Count Number of Points within Specified Distance from Lines using Carto?

texts:
There were multiple points stacked on the same coordinates that I did not notice.

I am trying to count the number of points within 10km of lines (routes) so the routes can be styled according to the number of points within proximity. However, the results are not turning out as expected.  Using Carto Builder, I created a 10km distance buffer around the points. On the lines layer I added the analysis for Intersect to the buffer layer with the count operation. However, when I check the results the count is incorrect (or rather, now what I want). For example, there are 6 points show below but results for the counts on the lines vary from 11, 12, 21, and 23.

I then tried the reverse to see why the results were behaving in this manner.  I created 10km distance buffers around the lines (routes). On the points layer I added the analysis for Intersect to the buffer of the lines with the count operation. These results were accurate.

I welcome any suggestions on a different perspective of the problem, even if it requires a SQL query instead of the built in analyses.


texts: Monthly statistics based on raster file name

texts:
I have a raster of a national park in Norway showing the temperature every day between the 1st of November 1999 and the 31st March 2015. I have created booleon images showing which pixels are above zero and which are below zero. I now want to go through the folder and sum up how many times each pixel has a temperature above zero each month. The file format is rr_YYYY_MM_DD.tif I tried adapting a script I found here (https://geonet.esri.com/thread/70044) to do this based on the file name, but I get the error message

"[function] + Utils.flattenLists(in_rasters_or_constants))"

when I try to run it.

The code is below.

import arcpy
arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = r"G:\Hardangervidda\Suitable pixels\Precipitation_3mm"
outWks = r"G:\Hardangervidda\Suitable pixels\monthly_stats"
for year in range(1999,2015):  # 1999 through 2015
    for month in range(13):      # 1 to 12
        # wildcard - "rr_YYYY_MM*.tif" (format month 1 -> "01")
        wild = "rr_{0}_{1:02d}*.tif".format(year,month)
        rasList = arcpy.ListRasters(wild) # list all rasters for a month
        outRas = arcpy.sa.CellStatistics(rasList,"SUM") # monthly mean
        outPath = outWks + "/" + "m{0}m{1:02d}".format(year,month)
        outRas.save(outPath)


texts: How to use arcpy.mapping to batch export layer to KML?

texts:
I am trying to export map layers from an MXD to KML using arcpy with the following code:

import arcpy, os, sys, traceback

arcpy.env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument(r"C:\Layer_KML\MXD\Test.mxd")

try:
    for lyr in arcpy.mapping.ListLayers(mxd):
        i = 1 # forgive the wonky naming convention, I just wanted something quick for test purposes
        outKML = "C:\\Layer_KML\\KML\\" + "Test" + str(i) + ".kmz"

        # These lines are where the error/s is occuring
        newLayer = arcpy.SaveToLayerFile_management(lyr, "C:\\Layer_KML\\Layers\\" + str(lyr) + ".lyr")
        arcpy.LayerToKML_conversion(newLayer, outKML, "20000", "false", "DEFAULT", "1024",     "96")

        i = i + 1

except:
        ...

Though my traceback keeps telling me : {ERROR 000623: Invalid value type for parameter in_layer.}  for both (I commented out each indivudally to test) the .SaveToLayerFile, and LayerToKML routines.  I am assuming these routines want to except a full file path, while the arcpy.ListLayer() only yields a basename.  Maybe?  I am really not sure.  Any thoughts as to why this isn't working?

I have already seen this question : How to batch export layers to kml in Arcmap and this is not what I am looking for.  I want to do this in Python rather than through a toolbox tool.

Any help is appreciated.  Thanks!


texts: Symbolic machine learning with OpenStreetMap data

texts:
I want to apply symbolic machine learning (e.g. case-based reasoning, inductive logic programming) to OSM data. In detail, I want to use OSM Map Features to predict characteristics of residential buildings.

Does anyone know algorithms, applications or related work to that problem?


texts: Performing QGIS Datum Transformation from WGS 1984 to NAD 1983?

texts:
Is there a way to transform a shapefile's datum from WGS 1984 to NAD 1983 in QGIS 2.18.0?

In the picture I have attached, the green dots are points that were collected in the field using a Trimble R2 receiver, which by default, collects coordinates in WGS 1984. The red dot is a benchmark with known GPS coordinates and is projected, along with the rest of my shapefiles in NAD 1983 StatePlane Virginia North FIPS 4501 Feet. Optimally, those green points would fall directly on top of my benchmark, but those points fall about 3.5 ft north of my benchmark. I believe this is due to the fact that they are projected in WGS 1984, while the rest of my shapefiles are in NAD 1983.

So far, all I can find in QGIS is ways of re-projecting the CRS. I believe what I need is a datum transformation. I am aware of the "Default datum transformations" tab within the QGIS settings, but if this is the solution to my problem, I cannot figure out how to use it.


texts: GDAL processing of GRIB2 raster prints hundreds of unsuppressable warnings

texts:
I am processing some GRIB2 raster files from NCEP (weather model stuff).  An example file would be:
http://nomads.ncep.noaa.gov/pub/data/nccf/com/hiresw/prod/href.20180906/ensprod/href.t00z.conus.mean.f06.grib2

When I use various GDAL tools to process the file it prints this line for each of the dozens of bands in the raster:

Warning: Master table version == 0, was experimental
I don't have a copy, and don't know where to get one
Use meta data at your own risk

Ok, that warning is fine -- but it's a major problem when I process the file with raster2pgsql as it dumps that message dozens of times into the SQL, which causes an error when I attempt to load it into the db.  Using gdalwarp -q or gdal_translate -q still prints the message hundreds of times despite being in quiet mode.

Is there any way I can suppress this warning, either with GDAL or some sort of pre-processing of the file?


texts: Configurable map viewer cannot display non-cached basemap?

texts:
I download the sample from github. Then I tried to add a custom basemap at the mapOption section. I have followed the instruction given in here. Also, I followed this post as well. But none worked, no map comes out. I only see empty map grid with "esri" logo at the bottom right corner. What have I done wrong ? How can I show a custom (non-cached) basemap in the mapOptions section.

Update :
I find out that if i use cached basemap (like http://services.arcgisonline.com/arcgis/rest/services/World_Terrain_Base/MapServer), it works. I have no luck with non-cached basemap. Then I tried to an older version (1.3.4) too. And its exactly the same, it woks only on a cached basemap.

Here is what I have done using non-cached map service:

Edit viewer.js to add the "Define" keyword and add "mapOptions"

define([  ... 'esri/geometry/Point', 'esri/dijit/Basemap',
'esri/dijit/BasemapLayer' ],

function (units, Extent, esriConfig, /*urlUtils,*/ GeometryService,
ImageParameters,GoogleMapsLoader,i18n, topic, has, Point, Basemap,
BasemapLayer)

    mapOptions:{
        // tiled layer
        basemap: **new Basemap({
            id: 'soco',
            layers: [new BasemapLayer({
                url: 'http://gis-webpub.sonoma-county.org/arcgis/rest/services/BaseMap/SonomaCountyBaseMap/MapServer/'
            })]
        }),
        center: new Point({
            x: 6356023.515330248,
            y: 1917313.443425702,
            spatialReference: {
                wkid: 102642
            }
        }),
        zoom: 5
    }**,

Am I missing something ?


texts:
Today, Google Earth allows you to import images to overlay on existing terrain.

Is it possible to import custom terrain to replace default terrain in Google Earth? There doesn't appear to be any options for doing so?


texts:
The reason I'm asking this question here is that I've seen a thread about ODL already. Let me know if I'm in a wrong place.

I'm trying to crosscheck a route plan using ODL Studio, which uses JSpirit to route vehicles. I have a customer data that is 400k rows long; and I keep getting memory errors.

I presumed it was because my data was too big; so I picked a portion from it and created 100k rows worth of data. No luck there as well. Here's where I'm confused.

I went ahead and deleted like 95k rows. Now with 5k rows, it created a matrix of 2000x2000 and worked perfectly. So instead of deleting random data, I created a new data with only one week's worth of information, which was around 4700 rows. For some reason my matrix is 4700x4700, and I can't get it to work again.

Can someone at least point me in the right direction so that I can have the optimizer working?


texts: Dynamic data loading in openlayers 3 using Geojson

texts:
I am able to load the multiple wfs features using GeoJSON, but now i need to load the features dynamically instead of loading randomly
Here is my code

// First WFS Layer
var url = 'http://localhost:8080/geoserver/wfs?service=WFS&' +
        'version=1.1.0&request=GetFeature&typename=WFS:Layer2' +
        'outputFormat=text/javascript&format_options=callback:loadFeatures' +
        '&srsname=EPSG:4326';
$.ajax({url: url, dataType: 'jsonp', jsonp: false});
var geojsonFormat = new ol.format.GeoJSON();
var vectorSource = new ol.source.Vector();
var loadFeatures = function(response) {
    vectorSource.addFeatures(geojsonFormat.readFeatures(response));
};
var fibreEquipment = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'rgba(55, 155, 155, 0.5)'
        }),
        stroke: new ol.style.Stroke({
            color: 'rgba(55, 155, 155, 0.8)',
            width: 2
        }),
        image: new ol.style.Circle({
            radius: 5,
            fill: new ol.style.Fill({
                color: 'rgba(255, 255, 255, 0.5)',
            }),
            stroke: new ol.style.Stroke({
                color: 'rgba(155, 55, 155, 1)',
                width: 2
            })
        })
    })
});

// Second WFS Layer

var url = 'http://localhost:8080/geoserver/wfs?service=WFS&' +
    'version=1.1.0&request=GetFeature&typename=Test:Layer1' +
        'outputFormat=text/javascript&format_options=callback:loadFeatures' +
        '&srsname=EPSG:4326';
$.ajax({url: url, dataType: 'jsonp', jsonp: false});
var geojsonFormat = new ol.format.GeoJSON();
var vectorSource = new ol.source.Vector();
var loadFeatures = function(response) {
        vectorSource.addFeatures(geojsonFormat.readFeatures(response));
};
var Layer1 = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style({
        fill: new ol.style.Fill({
            color: 'rgba(55, 155, 155, 0.5)'
        }),
        stroke: new ol.style.Stroke({
            color: 'rgba(55, 155, 155, 0.8)',
            width: 2
        }),
        image: new ol.style.Circle({
            radius: 5,
            fill: new ol.style.Fill({
                color: 'rgba(255, 255, 255, 0.5)',
            }),
            stroke: new ol.style.Stroke({
                color: 'rgba(155, 55, 155, 1)',
                width: 2
            })
        })
    })
});


texts: Is there a PostGIS function for determining whether a linestring intersects itself?

texts:
I've been looking for a function that can tell me if a LineString intersects itself. I've tried calling st_intersects with the same LineString twice, but obviously identical linestrings will intersect. My plan if there is no function would be to get all the points in the line string, and create individual linestrings and then check each new linestring against each other one with st_intersects. I do not want it to come to this, but I'm afraid it will.

So are there any PostGIS functions for checking if a linestring is self-intersecting? It should be similar to finding out whether a polygon is complex I'd assume


texts: Is the GIS Repository In Ubuntu in a different stream than the main repository?

texts:
I was checking out this:
https://wiki.ubuntu.com/UbuntuGIS
and it outlines how to install PPA for GIS.

I performed a basic apt-cache search on a few things and QGis and mapserver are available.

Such as:

$ apt-cache search mapserver
cgi-mapserver - CGI executable for MapServer
libmapscript-perl - Perl MapServer module
libmapscript-ruby - Ruby MapServer library
libmapscript-ruby1.8 - Ruby MapServer library
libmapscript-ruby1.9.1 - Ruby MapServer library
mapserver-bin - MapServer utilities
mapserver-doc - documentation for MapServer
php5-mapscript - php5-cgi module for MapServer
python-mapscript - Python library for MapServer
qgis-mapserver - Quantum GIS mapserver

Running linux mint Nadia which works out to be quantal for Ubuntu version (10.10)

Is there any difference for adding this package?  I notice that there is a testing version as well.  Is the package difference similar to stable, testing, and unstable in debian or are the packages the same?


texts: arcpy absolute file path failing with AddLayer()

texts:
I have the following code:

import arcpy.mapping as mapping

mxd = mapping.MapDocument("current")

print mxd.title

df = mapping.ListDataFrames(mxd)[0]
print df.name
layer = mapping.Layer(r"C:\path\clippedCountries.lyr")
print layer.name
mapping.AddLayer(df, layer,"TOP")

# Refresh things
arcpy.RefreshActiveView()
arcpy.RefreshTOC()

And it works fine.  However, if I try to use the absolute filepath to set the mapping variable (instead of "current") it doesn't work.  Any ideas?  I''m running the script through the ArcGIS 10.2 python window.

* update *

Interestingly if I use the full path name and add mxd.save() it works, so long as the mxd is not currently open in ArcGIS.  Next time I open the mxd the layer is there.

If I try to do this with the mxd referenced in the filepath OPEN at the time I get

IOError: MapDocObject: Unable to save.  Check to make sure you have write access to the specified file and that there is enough space on the storage device to hold your document.

Which I believe is due to a lock.

If I switch it back to "current" the lock vanishes and it will save and show the layer added in real time.


texts: How do I construct a geometry point in SRID 4326 from lat and long?

texts:
I have a geometry(Point,4326) column on my table:

CREATE TABLE foo ( geom geometry(Point,4326) );

I want to insert into it. Currently' I am running

INSERT INTO foo (geom) VALUES (
  ST_GeomFromText('POINT(latitude longitude)', 4326))
);

To save a geometry value. I think ST_GeomFromText('POINT(latitude longitude)', 4326) is wrong method, I tried from here


texts: Resolve Building Conflicts - the open source way?

texts:
In ArcGIS there is a tool called Resolve Building Conflicts and I'd like to achieve the same using only open source tools (QGIS, GDAL etc), preferably on Mac OS instead.

In the long run I'm trying to do exactly what Anna Vetter tries to do in her thesis called The automatic generalisation of buildings whilst maintaining the settlement structure: A case study
based on the 1:50‘000 Swiss National Map Series but without using ESRI.

I've also seen this discussion about buildings.

Any ideas?


texts:
Maperitive, a desktop OSM renderer, stores its rulesets for rendering in a subfolder called "Rules" - at least when running under Windows.

I tried to generate my own rules (or to be more honest I changed a ruleset) and saved it under a new name like "myruleset.mrules" in this folder. Unfortunately, Maperitive does not seem to "recognize" new files in this subfolder "Rules". Restarting Maperitive does not help.

How can I "register" this new ruleset, so that it is available under the menu Map -> Switch to Rules?


texts: Rotate a table but not the text in a QGIS Layout

texts:
I'm working in a QGIS Layout and I have an special requirement: I have to rotate a table but I need that the text remains horizontal. When I check the rotation  settings, the "Expression" option just let me to make an expression for changing the orientation of the table but not the text.

Is there any way to solve this tricky problem?


texts: Leaflet: Find my "this" in GeoJSON.pointToLayer function

texts:
In leaflet, I have a class that extends L.GeoJSON:

L.MyLayer = L.GeoJSON.extend
({
    initialize: function ()
    {
        this.markers = [];
...

I redefined the pointToLayer method to create markers for points:

pointToLayer : function(feature, latlng)
{
   var marker = L.marker(lat...
   ...
   return marker;
},

I need to store all the markers locally, so I created the markers array.
My problem is that if I call this. inside pointToLayer function I obtain the Window object instead obtaining the current instance of MyLayer.

How can I get a reference to current MyLayer instance inside pointToLayer function?


texts: How to install Quantum GIS Lisboa and scipy

texts: I'm trying to install some plugins of the QGIS but an error occurs. The program asks the scipy program. I installed Python, the Numpy and scipy. Also installed the program R. However, QGIS continues asking for scipy. Could you help me?


texts: Planet asset metadata different from search endpoint

texts: I ingested some PlanetScope 4-band Analytic assets into Google Earth Engine. I'm noticing the metadata of each scene is different compared to what the search_endpoint.py script from the Planet docs page returns. Are there some auxiliary metadata files somewhere/am I missing something?


texts: How to Use Custom Basemaps with BasemapToggle widget?

texts:
I'm trying to upgrade my app from version 3.3 to version 3.7 of the JavaScript API and am switching out the basemap gallery for the new BasemapToggle widget.

I have two custom basemaps that I want to use in my app using BasemapToggle. Is this possible? The API homepage suggests it is but I haven't seen a sample or any documentation on how to do it properly.

Can anybody point me towards a solution or example of how to use my own basemaps with this widget?

EDIT

I'm trying to do this with two basemaps, loading one of them as the initial basemap and then being able to toggle between them. My attempt at Derek's answer is HERE


texts:
Hi Guys is it possible to use regexp_split_to_table to split my column on 2 different levels. AS in split by ';' and then split by ' '.
I have this "-29.617356578819454 30.394328366965055 0.0 0.0;-29.61735976394266 30.39434471167624 0.0 0.0;" a column called "RPT_FTR_TRACE". ; seperates each coordinates and a space seperate long from lat.
I want regexp_split_to_table to split it for me so my final query have "-29.617356578819454 30.394328366965055 0.0 0.0; in one column and long=-29.617356578819454 column and lat= 30.394328366965055 column

WITH a AS (
         SELECT DISTINCT "_URI",
            "_SUBMISSION_DATE","INFO_TEAM_NAME","INFO_SUPERVISOR_NAME","INFO_LOCATION","RPT_FTR_WIDTH",
            "RPT_FTR_TRACE" AS t
           FROM odk_duct.duct_monitoring_core
        )/* selection from odk_duct.duct_monitoring_core table,*/
    ,b AS (
    SELECT a."_URI",
           a."_SUBMISSION_DATE",initcap(a."INFO_TEAM_NAME"),"RPT_FTR_TRACE",
           b.p,
           b.id

           FROM odk_duct.duct_monitoring_core a,
           LATERAL regexp_split_to_table("RPT_FTR_TRACE", ';'::text) WITH ORDINALITY b(p, id)
           WHERE b.p <> ''::text
            ) select * from b


texts: How to select polygons by duplicated attribute value from a shapefile?

texts:
I have a shapefile that contains many polygons(>200,000). Some of these polygons are assigned the same ID number. None of the polygons overlap other polygons. I need to find a way to select out all of the polygons that have the same ID number as another polygon for further processing.

I have looked into finding all the duplicates in an Excel sheet and then joining this to the shapefile attribute table but this would take almost as long as manually selecting each shape. I would like a way to select all these out only using QGIS.


texts: VBScript - If/Then Statement - Field Calculator

texts:
I am trying to do a simple (or so I thought) if/then statement in the ArcGIS field calculator.

I am trying to add the HydroCode field - hydologic soil codes - to an existing data set. I know what the code should be based on the "MUSYM" number. For example, "MUSYM" 1667 has the hydrologic soil code A. The "MUSYM" field is a number and the "HydroCode" field is text.

I vaguely remember doing if/then statements years ago when I more actively used GIS. I have tried simple If "MUSYM" = 1667 Then "HydroCode" = A Endif along with some other variations but none are working.

What am I missing?? All of the help I have found online deals with very complex coding, I can't find anything for such a simple statement. I also can't find anything with just basic operational instructions as far as where the operators should go/punctuation/etc.


texts: Map resize displays google layer at incorrect zoom level

texts:
This is a screenshot of the icons in the wrong position due to the google layer showing the wrong zoom level.

This is how it should look, and how it does look after the user pans the map.

How does this happen? I have a menu which opens a panel on the right side of the map. When the panel is displayed, the map div is resized accordingly. Usually this does not produce any problems. However, at certain zoom levels, resizing the map changes the zoom level which the google layer is displayed at. As soon as the user pans the map, the google layer image changes to the correct image for the current zoom level.

This behaviour doesn't always occur, and it only seems to happen in IE (8 and 9 at least). Resizing the browser window also resizes the map, and sometimes this causes the same behaviour. However the most predictable way to produce this has been to resize the map.

I am using SmartGWT and OpenLayers for the UI, Geoserver for the overlayed WMS layer. The map div is setup in a SmartGWT Canvas. I am handling the canvas resize event and calling map.updateSize() from there. As I have mentioned, the base layer is a google hybrid layer.

Although this problem doesn't occur 100% of the time, it does happen often, and it happens at the kind of zoom levels which users are likely to view very often.


texts:
I would like to delete several values as they are equal to zero. In my example, 0 -5 - 6 - 7 - 8 ...

I couldn't find an answer on the internet.


texts: Compressions to support when Implementing GeoTiff Reader?

texts:
I'm writing a GeoTiff reader for GeoTrellis for Google Summer of Code. I have read in all tags and geokeys and now I'm working with the decoding process.

Currently the Tiff 6.0 spec recommends that baseline readers support the following decompressions:

Uncompressed (duh)
PackBits
Huffman

And the extension recommends the support of:

LZW
Jpeg (Old spec)

I have implemented all the baseline compression types and LZW but as I understand it the compression type 6 (Jpeg old spec) is deprecated and replaced by compression type 7 (new Jpeg). Except these compression types mentioned Deflate compression (compression type 8) is mentioned in an extension document by Adobe.

As I understand it nobody uses the old Jpeg compression, is there any common other compressions I should support? That isn't bound for a single reader/writer software brand? And I should implement the compression type 7 and 8?

The goal with the reader is to support as many GeoTiffs as possible.

For the interested the reader is written in the Scala programming language.


texts: Combining two queries to compare different columns

texts:
Suppose I have a table

create table positions
(
 location geography,
 property1 int,
 property2 varchar
);

I can do the following queries:

select property2 from positions
  where st_dwithin(POINT, location, dmax) and not st_dwithin(POINT, location, dmin);

select property2 from positions where st_dwithin(POINT, location, dIgnore);

I want to

select * from positions
  where [results of query1 except where results of query1 have property1 = property1 of the results of query2]

I have tried various things (which I will spare you the trouble of reading, because they all resulted in syntax errors). How would I do this query?


texts:
I use QGIS SCRIT 2.14

When I run the script I would like to see the progress percentage

##Input_raster1=raster
##Band_raster1=output raster
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *
import os
from osgeo.gdalconst import *
import numpy
import sys
import osr
driver = gdal.GetDriverByName("GTiff")
ds = gdal.Open(Input_raster1, GA_ReadOnly)
dst_ds_red = driver.Create(Band_raster1, ds.RasterXSize, ds.RasterYSize, 1, GDT_Float32)
dst_ds_red.SetProjection( ds.GetProjection())
geotransform = ds.GetGeoTransform()
if not geotransform is None:
  st_ds_red.SetGeoTransform( geotransform )
bandList = []
band = ds.GetRasterBand(2)
data = band.ReadAsArray()
bandList.append(data)
rasterlayer_red = numpy.zeros((ds.RasterYSize, ds.RasterXSize), dtype=numpy.float32)
for i in range (0, ds.RasterXSize):
  for j in range (0, ds.RasterYSize):
    RED = (float(bandList[0][j][i]))
    rasterlayer_red[j][i] = RED
dst_ds_red.GetRasterBand(1).WriteArray(rasterlayer_red)


texts: Create a custom Vector Source for OpenLayers for referenced tiles

texts:
I have a custom vector storage, using Flatbuffers (if you’re not familiar: similar to pbf). It would take me a bit longer to describe the format in details, though it is different to other vector formats like mvt in that it is independent of a projection: Where other vector formats have to adhere to a certain projection and have to alter geometries (e.g. cut along tile bounds which is then also depending on the projection / tile bounds), my format always returns complete geometries.

This works already in a product available for iOS and Android. Though the culprit is that this format requires an extra step when fetching the tiles: Each tile may have referenced tiles, which also need to be loaded. This is because I use full geometries, and therefore tiles where geometries laying in one tile and overlapping other(s) have to be loaded as well.

I should I implement such behaviour in OpenLayers, loading depending tiles and also showing geometries of referenced tiles which might nit be in the current view (as per commonly implemented vector / raster sources), but containing overlapping geometries which need to be shown?


texts: Existing features not in Attribute Table and Data View

texts:
In ArcMap 10.7.1 I have drawn some lines in a already existing feature class containing other lines. The preexisting lines are displayed correctly in data view and attribute table. The new drawn lines are not shown either in attribute table and data view. If I guess the Object_Id of a new line and select by attribute, the line is displayed in data view and the table. But disappears immediately if unselected. Symbology is by categorie, unique values. Unique values are enabled. There is no definition query. Any body an idea?

I added this to screenshots. I hope this will clarify the issue. My question is, why are there lines which are shown only when selected?


texts: I keep getting the "Paste Features Failed." dialog when trying to copy/paste point features in a layer in ArcGIS pro

texts:
I have a "plant_species" point feature class, and after a period of time working with the layer, the ability to copy/paste an individual plant_species feature from within the map frame no longer works. I have been working on symbology for the map so I think it may have something to do with that.

The layer's symbology is set to "Unique Values", and has proportional size effects applied to it via the "Vary symbology by attribute" >> size menu, where I apply a multiplier to a radius attribute stored in the geodatabase.

The steps to reproduce are very simple:

I add the feature class data to the map, edit the symbology using custom *.svg files, and it works fine for a while. I can copy/paste a tree and it works just as you would expect.
After some amount of time working with the layer (creating new plant point features, copy pasting plant features, and editing symbology), I get a simple "Paste Features Failed" dialog box when I try to, for example, copy/paste a single plant point feature by selecting it, right clicking it, and copy/pasting in the map frame. There is neither a detailed explanation nor a link to help provided.
A (sort of) workaround is to remove the layer from the map and re-add
it, but this requires me to set up the symbology again, which takes
forever. So I have had to resort to using a "print map" which
contains my symbology, and a less fancy working map that I use for
editing.  I'm not sure if the copy/paste problem will happen in that map, which would suggest it doesn't have anything to do with the custom symbology.

I would prefer to use the print map as my editing map, as I am working with a lot of different species and the symbols really help.

Edit -- A better workaround:

Copy the layer with the problem.
Rename the original layer (important).
From the symbology tool of the new layer, click the hamburger menu, and then "import symbology...". Import the symbology from the original layer to the new one.
Edit any fields in custom expressions as necessary(some of these needed tweaks for me).

This is still not ideal, but much more workable than redoing all the symbology by hand, and persists any custom expressions you are using.


texts: Creating multiple polygons based on attributes in field using ArcGIS Desktop?

texts:
I have a file that contains a field with multiple pieces of attribute data for one entry. The data in the field is separated by a semi-colon.
I need to create a separate polygon for each attribute in that field. Is this even a possibility, or am I going to have to copy each polygon multiple times and delete the attributes I don't want from each entry?

Basically I want the info in the second column, to be broken into each, individual code to have its own polygon.


texts: How to identify the right road using postgres and postgis

texts:
I want to identify a road at particular point of intrest using post gis.

Example

Input String : 16th cross near xyz mall.

Problem I am facing

As there will be many "16th cross roads" in the map , how do I identify the "16th cross near xyz mall"


texts:
I'd try to merge multipart feature line, but didn't work.
would you please teach me how to merge multipart feature line?

I work in version environmental and geometry network.

Regards,


texts: Turning on developer options of Collector for ArcGIS on Android (Early Release)?

texts:
I am trying to turn on logging in order to troubleshoot an issue.

And found this is on GeoNet:

  Open the settings > tap your login name at the top 3 times >
  "Developer" options should appear at the bottom of the settings > Open
  the Developer options and turn on logging > reproduce the download
  failure > back in the Developer options, click Logs > you will see
  more detailed errors that you can review or share out via email, etc.

I don't know what device or version is referenced in the above-mentioned post but I don't have Settings. I did go into the About screen and started trip tapping things.  And I did get a developer options enabled notification message, but no actual UI change or additional options appeared.  And I'm not really sure what was actually being tapped that triggered that message.

How do you turn on the developer options?


texts: Changing point cluster symbol color by expression (color ramp with manual intervals)

texts:
I am working on cluster map (QGIS 3.0).
My map display both "not clustered" markers  with color attribute based on manual brackets ( [1-1.9); [2-4.9]; (5+ ) and clustered markers (red color).

What I want : I want to apply the same color attributes for my clustered markers using an expression.

What I tried (based on assistant) :


texts: Data for atmospheric correction with 6s from Modis

texts:
I want to make an atmospheric correction with Py6S. So, I want to get ozone, water vapor, AOT from Modis - MYD08_D3.

I got water vapor from dataset "Atmospheric_Water_Vapor_Mean". The average value is 0.0410000019474 in "cm". Did I take a correct dataset? And how to convert "cm" to "g/cm2"
For AOT I took dataset "AOD_550_Dark_Target_Deep_Blue_Combined_Mean" and got average value 0.175000008312 without any unit of measurement. Did I take a correct dataset? Should I covert this value to something?


texts: Obtaining file path in QGIS Training Manual

texts:
6.2.5 Try Yourself: Convert Layers’ CRS
...
8. Repeat this process for each layer, creating a new layer in
the vector_analysis.gpkg GeoPackage file with _34S appended to the
original name and removing each of the old layers from the project.

link: https://docs.qgis.org/3.16/en/docs/training_manual/vector_analysis/basic_analysis.html
I appear to have successfully added roads_34S to the vector_analysis.gpkg file:

My goal now is to repeat the process with restaurants, schools and buildings. I open the Save Vector Layer as … dialog and enter the Layer name and find the correct CRS. Now I want to use the ellipse button by file name to find the file path:

I search for vector_analysis.gpkg to obtain the file path. When adding the roads layer, I put the new package into a folder, myExerciseData, reserved for files that I create. I click the button and get the following, showing that the filenames are inactive.

Since the file names are inactive, I cannot open a file, which seems to mean that I cannot save the restaurants_34S layer with the new CRS to the file vector_analysis.gpkg. What have I missed? Do I have to enter the file path manually?
I try entering vector_analysis.gpkg. I click on Save. I get this dialog:

What has gone wrong here?
The following was added after Ben W's response.
I filled the boxes in the order in which they appear, as in the following image. The filenames remained inactive as above in the Save Layer as... dialog.

When I selected OK, I got a new dialog box:

Does this offer a clue? Are there any OS X users on this forum?


texts: Extract X Y data from Mapinfo layer data files?

texts: I have been provided with layer data files from MapInfo (.TAB, .DAT, .ID, .IND, .MAP). These seem to consist of polygons of various sites, but was told they include centre points of each site. I would like to extract any useful data from this including X Y coordinates into a form I can use ideally text. I do not have MapInfo. I have QGIS. Is there a tool or anything I could use to do this please? I am pretty new to the GIS world so please be gentle!


texts: Using java with mapinfo for integrated mapping?

texts: Can we use integrated mapping with java? Are there any examples?


texts: Showing lines in shapefile with different size of border depends on feature value

texts:
I'm trying to write a PyQGIS script.

I have shapefiles with lines which looks something like this:

ID   name   weight
1     a       15000
2     b       100
3     c       450
4     d       50
5     e       10000

And I need that lines depend on the weight show with different size of the border.

I have the dataset with 100 lines and I need just first 10 to show with different weight of borders (first to have the biggest border, second little less than that, etc), and rest to set to be invisible (or to colour in white).


texts: WKT: What is the reasoning behind the concept of a POINT EMPTY?

texts:
In data formats like WKT and WKB, we have ways of representing "empty" geometries. That all pretty much makes sense: a "LINESTRING EMPTY" is a LINESTRING type with 0 vertices. Same with a MULTIPOINT, POLYGON, etc.

But what about POINT EMPTY? For starters, there is no WKB representation for a POINT EMPTY, only WKT. GeoJSON also doesn't define this case in its specification.

According to the Wikipedia definition of a geometric point, "points do not have any length, area, volume, or any other dimensional attribute". So how can it be empty, and why does the WKT specification allow this?


texts: Displaying the sum of values within a given polygon

texts:
I am using some sample data to try to figure out the best method to complete the desired task so that I can ensure all the proper data is transcribed in the coming weeks. I have hundreds of reports that include a proposed number of roads that may or may not have been constructed. There is little post-project reporting or much for ground-truthing on the part of compliance or enforcement agencies, so there is no linear spatial data associated with these proposed roads. All that we have is a general project location (point location), and the total proposed kilometers of road or trail. I would like to display the total proposed kilometers of roads, for each project that falls within a subsubdrainage (Polygon), or a given 10x10km grid (polygon).
I would like to calculate the sum of attributes from a point later that fall within an overlapping grid layer (polygons), to then display a gradient of proposed total kilometers of road within a given area.

I have tried working with spatial join in ArcMap to combine the grid layer with the point layer, but I was not able to find a means to calculate and display the sum roads.

After doing some research I found the Statistic by Category Tool in QGIS.
It looks like it could calculate what I want, and then I could potentially (?) join the results to the polygon grid and display as a gradient.

The issue is that when I tried to use the statistic by category tool the sum and range are not included in the results table. When I was doing my research it appeared as the two attributes were included in other people's results. I updated my QGIS from 3.6 to 3.12 but it did not make any difference. Am I missing something in order to obtain these results?


texts: Unknown proj4: how to create custom proj4 from wkt?

texts:
I have a shapefile with a .prj file I've never seen before. In fact, it seems that GDAL doesn't even know what it is (see below).

Therefore, when I read it into R with raster::shapefile() the CRS isn't set. When I set it to the closest thing I can find (ESRI:54030, or ESRI:53030), the projection is still off (even when I modify the Central Meridian with "+lon_0=11" in the proj4string).

Here's the WKT:

(base) Richs-MacBook-Pro:whymapdata_world richpauloo$ gdalsrsinfo wrb11hyg_7d84d.prj

ERROR 6: No translation for Robinson_ARC_INFO to PROJ.4 format is known.
PROJ.4 :

OGC WKT :
PROJCS["WGS_1984_Robinson_ARC_INFO",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137.0,298.257223563]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Robinson_ARC_INFO"],
    PARAMETER["False_Easting",0.0],
    PARAMETER["False_Northing",0.0],
    PARAMETER["Central_Meridian",11.0],
    UNIT["Meter",1.0]]

If I compare this to the WTK for ESRI:54030 and ESRI:53030 they're clearly not the same.

What projection do I have? How can I set the projection for this object in R for either an sp or sf object? Do I need to write a custom proj4 string from this WKT?


texts:
I wrote some Tools for CAD in QGIS with python. Therefore I add a special temporary layer with python, if this layer does not exist. That works fine. To add a new feature I start edit, add the feature, end edit and commit. The feature appears and is added. I can identify it and snap it. That works fine too. But then I cannot snap to every further feature I add to the same layer with the same code, until I open the project again. All these features are visible and can be identified. Only snapping fails.
I think it could have to do with the feature ID,  but I have no idea what to do.
Snapping is activated for all layers, before I start the plugin.

The code to create the layer:

def createMeasurePointLayer(iface, LayerName):
prCRS=iface.mapCanvas().mapSettings().destinationCrs().authid()
vl = QgsVectorLayer("Point?crs=" + prCRS, LayerName, "memory")
pr = vl.dataProvider()
vl.startEditing()
pr.addAttributes( [ QgsField("fid", QVariant.Int),
                    QgsField("Punktnummer", QVariant.String),
                    QgsField("Foto", QVariant.String),
                    QgsField("geändert von", QVariant.String),
                    QgsField("geändert am", QVariant.DateTime),
                    QgsField("erstellt von", QVariant.String),
                    QgsField("erstellt am", QVariant.DateTime),
                    QgsField("erstellt durch Methode", QVariant.String),
                    QgsField("geändert durch Methode", QVariant.String)])
vl.loadNamedStyle(u"C:\\OSGeo4W64\\apps\\qgis\\python\\plugins\\EingabeMesspunkt\\Messpunkte.qml",False)
vl.endEditCommand()
vl.commitChanges()

QgsProject.instance().addMapLayer(vl)

The code to add a feature:

...
    fet = QgsFeature()
    fet.setGeometry(pt_z)
    fet.setFields(fields)
    fet.setAttribute(0,idx)
    fet.setAttribute(1,PNR)
    if PointPhoto != "":
        fet.setAttribute(2,PointPhoto)
    fet.setAttribute(3,User)
    fet.setAttribute(4,strTime)
    fet.setAttribute(5,User)
    fet.setAttribute(6,strTime)
    fet.setAttribute(7,PointMethod)
    fet.setAttribute(8,PointMethod)

    layer.startEditing()
    layer.dataProvider().addFeatures([fet])
    layer.endEditCommand()
    layer.commitChanges()

...


texts: How to List Tables from current map document (mxd)?

texts:
We have many maps containing Feature Classes and Tables in a Personal Geodatabase that need their information updated.

Is there a function to List Tables from the current map document like

mxd = arcpy.mapping.MapDocument("Current")
arcpy.mapping.ListLayers(mxd)

in ArcGIS 10.1?


texts: QGIS expression to access page-number of atlas for certain features

texts:
This question is based on another one that found only a partial solution (workaround) - but there should be a better solution, I'm convinced. However, I'm stuck. Here is what is already clear:

In an atlas, you can get for each page of the atlas output a list of the elements (say: towns) on this particular atlas page - a list of the towns that appear on this particular atlas page. So for every page you get a list of the towns that appear on this particular atlas page.

Now here is the question:

How is it possible compose (let's say on a separate layout-page) all features (all towns) and list them together with the no. of the atlas page they appear on? In the most simple case you have one layer towns and one layer rectangular polygons, the last one is the one that creates the atlas pages.

As a starting point I thought about using an expression used on the towns-layer, based on intersects( $geometry , @atlas_geometry ), however I'm not sure how to integrate the page-no. The variables @atlas_featurenumber and @atlas_featureid can both be used to refer to the current page-no on each individual atlas-page, but i don't see how it could be used the refer to all the pages (not just the current one) in a separate list. I thought of using an array, but I'm stuck with how to set it up and use it accordingly. Final aim is to create a register of all towns that lists the page-no. of the atlas on which they appear, like:

town a  page 4

town b  page 7

town c  page 19

town d  page 4


texts: Why ST_Intersects() returns one row more than expected?

texts:
I am trying to understand the behaviour of ST_Intersects() in PostGIS.

My goal is to get the mean value of a raster area corresponding to each polygon in a vector layer.

To test my SQL I run initially the script with a vector layer containing only one geometry:

CREATE TABLE mean AS
SELECT id, (ST_SummaryStats(ST_Clip(rast,1,geom))).mean
FROM vector, raster
WHERE ST_Intersects(geom,rast)

What returns is a two rows table:

1;-146.258943781942
1;-115.587301587302

How it is possible that the same geometry intersects the same raster twice with two different values?

I don't know if it is a coincidence, but the mean value I got doing the same operation with the same two layers with the QGIS Zonal Statistics is -138.788659793814 (the mean between -146.258943781942 and -115.587301587302 is -130... pretty close)


texts: What is the best method for route cost estimation?

texts:
What would be the best method for calculating a rough adhoc route cost i.e. a potential road

(been some years since I have used Network Analysis in ArcGIS), alternatively would it be better using Google Map/ArcGIS Javascript API ..(I am okay with really rough cost)

Ideally freehand e.g. http://gmap-pedometer.com/ and not just between 2 specific places or abiding by an already set network like network analysis requires
It would be intelligent i.e above a certain relief the route would be greater in cost per Km (if a tunnel is required), then as it declines below X metres it would resume normal cost per Km


texts: Merge node results into one MultilineString

texts:
I'm trying to generate MultiLineStrings from street data nodes. Right now my workflow is as follows: I run something like the following through the Overpass Api, to get all the nodes:

[out:json][timeout:25];

{{geocodeArea:Berlin}}->.searchArea;

(
  node["highway"]["name"="Holteistraße"](area.searchArea);
  way["highway"]["name"="Holteistraße"](area.searchArea);
  relation["highway"]["name"="Holteistraße"](area.searchArea);

);
out body;
>;
out skel qt;

Then I export it as geojson and merge all features into one MultiLineString with turf.js. I then add this generated data as MultilineString into Mapbox.

I'm quite new to this, so this workflow probably isn't the best way to do it.

How could I directly call for a MultilineString via the overpass API?


texts:
I am unable to call shp2pgsql using python as I  have to pass various options like (-s -a -i)
and I am not getting the required output ..

Method I am using is

p1 = subprocess.Popen(["cmd", "/C", "shp2pgsql","-a","-s 4326","roads.shp" "roads", ">","roads.sql"],stdout=subprocess.PIPE)

>>>p1.communicate()[0]

as my intension is to add the features of the shapefile into the already existing table in the database

I am done with the environment settings..
roads.shp is the name of the shapefile
roads is the table name
and roads.sql is the sql file name..


texts: Merging .tif files downloaded from Google Earth Engine in correct order in Python

texts:
When you download a large image from Google Earth Engine, it is in multiple small images. Each sub-image name contains the x_min and y_min coordinates of that image.

It is described here - https://developers.google.com/earth-engine/image_upload#tfrecord

How to join all these .tif files into one .tif file in right order by using a Python script.


texts: How to obtain coordinates from points with similar symbol in an image?

texts:
I would like to obtain a table containing all the coordinates of the triangles in the following image.

I know it is possible to georeference my map using the georeferencing plugin and then I can access the coordinates of any point using the Copy_Coords plugin. But this takes much time and I would like to know if you know a way to do it faster ?