texts: ogr2ogr fails to convert GeoJson file > 200MB

texts:
As a followup to my previous question (Converting SimpleGeo Places GeoJSON file to shapefile), I am now struggling with what I perceive might be an ogr2ogr "memory" limitation, but I cannot be sure.  I have successfully converted the SimpleGeo dump files into proper geojson files (verified that are directly readable in QGIS).

Now I want to convert this data to SHP or other formats (SQLite, PostGIS).  I am running the following command with success on smaller files (where "places_fr.gejson" has a file size of ~118MB):

ogr2ogr -f "ESRI Shapefile" places_fr.shp places_fr.geojson

When I attempt to run this same exact command with a geojson file that is 340MB, ogr2ogr crashes in one of two ways:

"lh_table_new: calloc failed" (first attempt)
Windows throws an "ogr2ogr.exe has stopped working.  Please close the program." (subsequent attempts)

I am running the ogr command via FWTools in a Windows 7 environment dual core with 8GB of installed RAM (4GB available, 12GB virtual memory available).  Is there a memory limit to ogr2ogr converting these datasets?  I have also tried converting into SQLite with the same results, so I don't think it is an SHP conversion issue.  Thanks.


texts: Georeferencing a GeoTIFF using with RasterIO

texts:
I am trying to Co-register a raster using Ground Control Points w.r.t. another raster using RasterIO. I have written a PoC, but having a bit of difficulty in measuring the quality of the final transformed raster. I do not have any ground truth transformed raster to verify it against but using QGIS tool and GDAL APIs are giving completely different and seemingly off results
I have manually identified few common points in both the rasters and now trying to setup the whole pipeline for the co-registration.
Terminology:
Query Raster: Raster to be co-registered
Reference Raster: Reference raster to co-register to
Following are my GCP-pairs (pixel values):
[[query_x, query_y], [reference_x, reference_y]]
[[4891, 2922], [4701, 2731]],
[[4902, 2898], [4710, 2703]],
[[4826, 1998], [4631, 1795]],
[[4551, 1293], [4360, 1102]],
[[4605, 733], [4413, 537]]

My rasters are in WSG84 co-ordinate system (crs = {'init': 'epsg:4326'})
So for doing the transformation, I am projecting my reference pixels to the WSG84 co-ordinates
reference_longitude, reference_latitude = reference_raster_dataset.xy(reference_y, reference_x)

and initialising the GCPs
rasterio.control.GroundControlPoint(row=query_y, col=query_x, x=<reference_longitude>, y=<reference_latitude>)

After this, am using
 transform = rasterio.transform.from_gcps(<list of GCPs>)

to get the Affine transform and re-projecting each band using the transform computed above.
So, my final doubts are:

Is my Ground Control Point format and initialisation correct?
Is there any metrics that I can use to verify the quality of the co-registration
Is the flow correct, because I am not able to find any working code flow etc, so just followed and implemented the theoretical guides to the co-registration process


texts:
I have a vector layer where some features are found entirely within other features and therefore I can't select that feature because the other feature is on top of it. Is there a way to reverse the order of the features of a vector layer after they have already been loaded? This will likely solve about 99% of issues that may occur with my mouseover popup.

Edit to add code:

I'm trying to do this on demand after the layer is already loaded.

var refreshStrategy = new OpenLayers.Strategy.Refresh( {
    interval : 60 * 1000
  });
var vLayer = new OpenLayers.Layer.Vector("Hazards with Population",
      {
        projection : "EPSG:4326",
        strategies : [new OpenLayers.Strategy.Fixed(), refreshStrategy],
        protocol : new OpenLayers.Protocol.HTTP(
        {
          url : "./getFeatures.php",
          format : new OpenLayers.Format.GeoJSON(),
          rendererOptions : {
            zIndexing : true
          }
        })
      });


texts:
I'm working in a project and i need to load some layers at startup.
I do it adding this code to startup file localted in the bin folder.

capa = QgsVectorLayer('Polygon?crs=epsg:32614&field=yeah:string(30)&index=yes', 'capita', "memory")
QgsProject.instance().addMapLayer(capa)

This step is right, however, when i start qgis, it asks if i want to save, discard o cancel changes.
If i select save or discard, it creates a new project without my layers.
but if i chose to cancel changes, nothing else happens and my layers are there ok (and this is what i want)

Does anybody knows how to hide that prompt and let my layers to be loaded without ask me to save?


texts: Calculating areas by overlay analysis of two tables with several information

texts:
I have to calculate the areas (m2) of the blue, green and yellow polygons within different vegetationtypes (also polygons - the table is overlaying the colormap). I have to calculate how large areas of the blue, green etc. polygons are within the different black symbols.

Can I calculate this in one SQL for all the color polygons?


texts: MapBox GL-JS vector tile - line rendering issue

texts:
I'm developing a MapBox GL-JS application and I faced with a problem regarding line styles. I want to display railway lines in high zoom level (20-22). I have one source, with two layers. The first is a black line, the second one is a white, dashed line. It's important that the line geometries are long, so on the pictures, there are no connecting points.
The problem is the rendering quality. See the images below. The first is a view with pitch:60, the pitch of the second one is 0.
The source is vector tile published by GeoServer v2.18.0 but I tried with pg_tileserv too. The rendering quality is the same.
What do you think? Do you have any suggestions how to fix it? Or is it a bug?
This is my code:
map.addLayer({
  'id': 'railways',
  'type': 'line',
  'source': 'railways','source-layer': 'railways_osm_3857',
  "paint": {
  'line-width':20,
  'line-color':   '#000000'
    }
});

map.addLayer({
  'id': 'railways2',
  'type': 'line',
  'source': 'railways','source-layer': 'railways_osm_3857',
  "paint": {
  'line-width':15,
  'line-dasharray':[5,10],
  'line-color':   '#ffffff'
    }
});


texts: Converting FeatureLayer to graphics using ArcObjects?

texts:
How do I convert a Map or FeatureLayer to graphics on ArcObjects 10.1?

Problem Solved

this is the working code:

    private void Convert_to_graphics( IMapDocument doc)
    {

        IMap pmap;
        IFeature pFeature;
        IElement pElement = null;
        IMarkerElement pp = new MarkerElementClass();
        IGraphicsContainer pGraphicsContainer;
        IActiveView pActiveView;

        pmap = doc.get_Map(0);
        pActiveView = doc.ActiveView;
        //pGraphicsContainer = doc.ActiveView.GraphicsContainer;

        int[] array = {8 , 7 , 6 , 5 , 1}; //layers indexes ,

        pGraphicsContainer = (IGraphicsContainer)pmap;
        foreach (var i in array)
        {

            IGeoFeatureLayer fl = (IGeoFeatureLayer)pmap.get_Layer(i);
            IFeatureClass fc = fl.FeatureClass;

            IFeatureCursor fcur =   fc.Search(null, false);

            IFeatureRenderer fr = fl.Renderer;

            pFeature = fcur.NextFeature()  ;

            while (pFeature  != null)
            {

                if (pFeature.Shape.GeometryType == ESRI.ArcGIS.Geometry.esriGeometryType.esriGeometryPoint)
                {
                    IMarkerElement pmkElement = new MarkerElementClass();

                        ISymbol sym = fr.get_SymbolByFeature(pFeature);

                        pmkElement.Symbol = (IMarkerSymbol)sym;

                        pElement = (IElement)pmkElement;

                }
                if (pFeature.Shape.GeometryType == ESRI.ArcGIS.Geometry.esriGeometryType.esriGeometryPolyline)
                {
                    ILineElement pmkElement = new LineElementClass();
                    pmkElement.Symbol = (ILineSymbol)fr.get_SymbolByFeature(pFeature);
                    pElement = (IElement)pmkElement;
                }
                if (pFeature.Shape.GeometryType == ESRI.ArcGIS.Geometry.esriGeometryType.esriGeometryPolygon)
                {
                    IFillShapeElement pmkElement = new PolygonElementClass();
                    pmkElement.Symbol = (IFillSymbol)fr.get_SymbolByFeature(pFeature);
                    pElement = (IElement)pmkElement;
                }

                if (pElement != null)
                {
                    pElement.Geometry = pFeature.Shape;

                    pGraphicsContainer.AddElement(pElement, 0);
                }

                pFeature = fcur.NextFeature();
            }

      }

            pActiveView.Refresh();

    }


texts: I'm having a mental block and was wondering if there was any way to click on an individual raster grid cell in ArcGIS and create a sample point readout of the data (other than in the identify tool).  I want to place sample points of numerous grid points similar to this (this is a different program).


texts: Exporting PNG image with background transparency from ArcGIS Pro

texts:
I am attempting to export a PNG image file from ArcGIS Pro 1.2.0 with the background set to transparent. In ArcMap 10.3 this is simple, there was a button for setting the transparent color to white during the export.

When I go to "share" the image in ArcGIS Pro, I do not see the same options to set my background or transparent color. Has this button been moved elsewhere? Or is it not an option anymore?

The goal is just to find similar settings to the screenshot above.
For this project, I have a basemap as a standalone PDF file. I need to overlay a callout box over top of that file--this requires that the callout box has transparancy all around it, and not a white box (like what happens when you clip to extent).


texts:
I'm attempting to calculate vegetation indices in ENVI 5.2 using level-2 AVIRIS data in band math.

The indices I'm working with are Vogelmann Red Edge 2, normalized difference red edge, and red edge NDVI (reNDVI). All these indices have the same format as for NDVI.

In band math I am using: float(B1-B2)/float(B1+B2)

These indices are supposed to all be in the range of -1 to 1, however they are typically outside of the -1 to 1 range.

For example, I tested this band math equation for NDVI and got a range of -1.163 to 0.946.

I've tried other iterations of the band math equation (e.g., altering which variables are designated "float"), but they also produce indices that are outside of the valid range.

What am I doing wrong?

EDIT:

Re-doing the indices in ENVI, most values are in the proper range. The table below is for reNDVI using (float(b1)-float(b2))/(float(b1)+float(b2)).
The second column, "count" indicates there are almost 19 million pixels with a value near -0.24, with another 10 million pixels in the two bins next to it. However, a couple hundred pixels are still out of range. I'm not sure why I get this. Does this make the whole band invalid (for the index), or are these really low values to be expected from things like bad pixels?


texts: Where is the drawn feature when the drawend is triggered in OL3

texts:
I have written a piece of code to check if the drawn line is intersecting other lines in the vector layer. If so, then the line should be split and new lines should be added to the layer. here is my code:

The events:

lineStringdraw.on('drawend', function(e2) {
var featureEnd = e2.feature;
var featureGeom = featureEnd.getGeometry();
var endNodeCoordinates = featureGeom.getLastCoordinate();
snapIndex= endNodeCoordinates.length-1;
checkSnapping(nodes, endNodeCoordinates);
if(snap){
    var featureCoordinates = featureGeom.getCoordinates();
    featureCoordinates[featureCoordinates.length-1]= snapToNode;
    featureGeom.setCoordinates(featureCoordinates);
    featureEnd.setGeometry(featureGeom);
    featureEnd.changed();
}
checkIntersection(featureEnd);
if (tempList.length!==0){
    //lineStringdraw.finishDrawing();    //my first solutuion
    //vector2.getSource().getFeatures().removeFeature([vector2.getSource().getFeatures().length-1]);
    //delete e2.feature;   //my second solotion.
    for (i=0;i<tempList.length;i++){vector2.getSource().addFeature(tempList[i]);}
}

drawing = false;
drawing_feature = null;
});

checkIntersection():

this function checks if the currently drawn line is intersecting other lines. if so then it calls the split function which will produce final lines to be added to the layer.

There is only one problem. The currently drawn line is also added to the vector layer and I cannot stop it as this feature is not yet added to the source so I cannot remove it from source. And also I cannot remove it from the draw interaction (at least I dont know how to remove it). It seems that the line is added after the drawend event is handled but how I can prevent it from happening. Any Idea?

EDIT

I now use the addfeature event for the source. when the feature is added, it will be checked for the intersection. but there rises another problem: when I add the new features from split function they will be checked again and again. I tried .unByKey() to remove the listener so that I can add the new features but then I cannot reactivate it.

var addFeatureEventListener = vector2.getSource().on('addfeature', function(e){
var currentFeture = vector2.getSource().getFeatures()[vector2.getSource().getFeatures().length-1];
var newLines =  checkIntersection(currentFeture);
vector2.getSource().unByKey(addFeatureEventListener);
console.info("the listener is deactivated now....");
vector2.getSource().addFeatures(newLines);
vector2.getSource().changed();
});

I am kinda fed up. I checked this answer but I seem cannot understand it.


texts: How to make Field Value constant size (length) by appending zeros towards end?

texts:
How do I append my field values in my .dbf to correspond with the values of my feature's attribute table.

I need to append 0's to the field value if the field value is less than 11 characters in length.

It needs to be done in the field calculator through VB Script or Python.


texts: Matching to symbols in style using ArcPy?

texts:
I have made a style with all my values and colours for my floodmap. I can match them manually at the properties, but is there a way (maybe with Python) to do this automatically in ArcMap?

When they add some new floods to the map it would be great if the colour based on the total floods will change automatically, because these floods are in despite a non-stop problem.

Does somebody know something for this?


texts: ogr2ogr KML to SHP generates invalid file names from layer names

texts:
I have a collection of very large KML files that I'm trying to convert to .shp. I'm iterating over 200+ KML files with an average size of around 20 MB. Each KML contains multiple layers, some of which have names that contain special characters including /.

When I run ogr2ogr on these KMLs the layers with names containing special characters do not get created, as the layer's .shp files are named by each layer and filenames cannot include things like /.

I wrote a script to replace all special characters using sed but due to the huge volume of text this looks like it will take several days to complete, even on a fast machine with SSD.

I was hoping to use ogr2ogr's -nln parameter to map outgoing layer names to character-safe names, but the parameter seems to only accommodate a single layer name.

Can anyone suggest a way to change the KML layer names before conversion that may be quicker than sed (Google Earth struggles and eventually crashes if I try to open the KML directly), or a way to explicitly tell ogr2ogr what to call the new .shp files?


texts: Error when setting "permanent" tooltip option with bindtooltip in Leaflet

texts:
I've got a geoJSON that I load into a map and add tooltips to. Pretty simple, figured that out fine, and now I want the tooltips to be permanent (displayed on map load instead of on marker hover), and I'm getting an error when I try to set that option in bindTooltip.
var newLayer= new L.geoJSON().addTo(map).bindTooltip(function(layer) {
   return layer.feature.properties["Property"];
}, {permanent: true, opacity: 0.5});

...later...

$.ajax({
   dataType: "json",
   url: "somewebaddress",
   success: function(data) {
   $(data.features).each(function(key, data) {
   newLayer.addData(data);
   });
   }
}).error(function() {});

If I remove the permanent: true option, I have no issues, but of course I do not get the desired permanent tooltip. But when I have that option in there, I get the following error message:

My theory was first that permanent: true makes it want to assign the tooltips instantaneously and that since the geoJSON didn't have any data assigned to it yet that was causing the error. But I tried binding the tooltip AFTER adding the data and again, if I have permanent: true in there it errors but without it I have no issues, so I guess that's not it.
I'm new to Leaflet, but I feel like this should be simple and I'm at a loss for a solution and struggling to find documentation.


texts: Raster to postgis - choosing right projection

texts:
I have a raster with the projection:

+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

I'm trying to insert a raster file into the postgis with a SRID 4326, but I don't know if it is right.

raster2pgsql -s <srid> -I -C -M <pathRaster> -F  -f rast public.rasterTable

How can I figure out the right SRID to every projection?


texts:
When I'm tryong to merge two line segments in QGIS, I get this error:

Invalid result: Could not store value 'Skipped' in field of type
Integer64.

What does this mean?
One example is showned in the picture below. I select two features, go to Edit > Merge selected features, and then the error.


texts: Difference of Assign projection, Reproject Layer, and Define layer projection in QGIS?

texts:
I have a point layer and a polygon layer that are all in UTM51N projection but do not align. Just want to know when should I use Assign, Reproject, or Define.

Experimented with a polygon in EPSG:4326 - WGS 84 - Geographic, which is the project CRS since this polygon is the first layer I loaded. Target projection is UTM51N.

Reproject - aligned

Assign - not aligned

Define - not aligned

Save layer as - aligned

Using 3.4 so there is no "on the fly projection" feature.


texts: QGIS Wien 2.8. Openlayers plugin doesn't work

texts:
I have installed the latest version of QGIS wien 2.8. Windows 8.1 downloaded the stand alone version from official QGIS website.

I had to uninstall the latest one I had because Openlayers plugin was downloaded and installed but one day it stopped appearing in my toolbar.

So I uninstalled and removed everything and installed again, but the problem remains, just as with the previous one.

Can anybody help me with this issue?


texts: Assigning id from a polygon to new field in intersected points using the Graphical Modeler in QGIS

texts:
I’m working with the Field Calculator in the Graphical Modeler and I’m stuck. I need to assign the feature id from a polygon layer (track) to a point layer (StopPoint) so that I can group together those points to later be connected by a line. The suggested expression works flawlessly in the Field Calculator from the toolbox but fails when run inside the Modeler. I get the error: Layer 'track' could not be loaded.
array_to_string (
    overlay_nearest(
        'track',
        id
    )
)

Is there a solution to this problem, or a workaround using other algorithms in the modeler?


texts: How to make appear values above the mean on rule based labeling QGIS

texts:
I've been looking at the other questions relating to rule based labeling but I haven't found any answer to my query.
I have polygons with average literacy rate for each canton of a country and therefore have a field named "LitRate". I'd like to show cantons with a rate higher than 5.
Getting inspiration from other questions I've tried:
"LitRate"= > "5" but it doesn't work.


texts: How to export an image from Print Composer with transparent background QGIS

texts: I am trying to export a map with a transparent background in QGIS (Using a slightly older version). I looked all around in the print composer but I couldn't find a place to set the background to transparent. I have already unselected the "background" box so I don't think that's it. If someone could direct me to the menu to make the background transparent, that would be great!


texts: QGIS crashes when changing projection to South Pole Lambert Azimuthal Equal Area

texts:
I can change the projection to every projection I have tried successfully except for South Pole Lambert Azimuthal Equal Area. When I click apply the loading bar appears in the bottom of QGIS and says 'loading tiles' it does this for several minutes before finally crashing. I went ahead and upgraded to the most recent version of QGIS to see if that fixed it and it did not.

Does anyone have experience with a similar problem or know what might be going wrong?


texts: What are some options for embedding Google Maps into a WordPress site?

texts:
Background
A non-technical user would like me to embed a Google Map into a WordPress-powered site. She needs the ability to add markers using a simple mechanism, and to link those markers to pages within the site.
Ideally I (as administrator) would have the ability to

use Styles on the basemap (eg, don't show roads)
write JavaScript which controlled map interactions (eg, restrict zooming to just 2 levels)
access advanced Google Maps functionality (eg, marker clustering)

There is no shortage of Google Maps plugins for WordPress, so my questions are:

Can anyone on GIS.SE recommend any particular plugin?
Do any of them allow customised JavaScript to run? (they appear to be aimed a non-GIS
endusers)

Thanks


texts: Creating map with same-area hexagonal polygons using ArcGIS for Desktop?

texts:
I'm trying to figure out

(1) what this type of map linked here is called (I've been unable to find the name through search), and

(2) whether there is an Arc tool that will transform a shapefile or feature class with geospatially referenced polygons into a map like this where the relative spatial positions of the polygons are preserved, but the areas are equal and the shapes are transformed?


texts:
I have a taxlot that is giving me troubles in a Union.  When I try to run the union I get the tried and true ''Geometry is not M-Aware'' message. (there are no Polygon-M features).
Based on readings on this site and the ESRI forums, I'm now trying to explode my taxlots prior to the union.
However I seem to have come across a mystical indestructible shape.
I've tried check and repair geometry with no results.
I've tried the Multipart-to-Singlepart tool and receive this message:

Cannot process feature with OID ...
Failed to execute

So I tried the explode function in the Advanced Editing Toolbar and go this lovely message:

I calculated (!SHAPE!.partCount) the number of parts into a field and came up with 22.

Here is what the shape looks like:

Help?


texts: Point features disappear when imported to geodatabase

texts:
I am trying to organize my files by creating geodatabase and importing shapefiles into it. I have created a feature dataset for my vectors which includes a few polygon layers, a line layer and points collected with a GPS. Before importing the files I have projected all of them to Ky state plane. Except point data everything looks fine. point data was collected originally with WGS84. I do not know if it was the reason.  My point data missing table and when I try to load features on it I get this error. Any idea why?

Warning 594 says

Description
These specific input features require attention. This often occurs when features are positioned outside the current XY Domain environment. The XY Domain environment sets the allowable coordinate range for x,y coordinates. If a feature falls outside this range, it cannot be added to the output. In addition, with some tools such as Append, which have a target dataset, if the features of the input fall outside the XY Domain of the target, they also will not be added. Similarly, the M Domain and Z Domain limit the coordinate range for m- and z-coordinates.

Solution
Check the XY Domain, M Domain, and Z Domain environments, and reset them to Same as Input' if they are too restrictive for your inputs. When using a tool like Append, check the coordinate system of the target dataset, and if unknown, define a coordinate system for the feature class. Also, consider creating a new feature class with a broader x,y domain if the current x,y domain is too restrictive.

I have checked domain it looks like the points within the boundary.


texts: What is the syntax for Python 2.7 to permanent display in percent in a field column?

texts:
What is the syntax to build to calculate field, for example,

Acres to Percent?

841173.00 to 82 %

What I want the field column in Percent to be display numbers with  percent with the alignment to the right of 2 characters? How is that possible?

I can run from the Calculate field toolbox and then use the Field Calculator and insert the expression with the python syntax.

I had to manually right click on the Percent field column's properties to get to the Percentage under the Number Format and change it from there. I feel this is a time-consumed for me because I will be doing this for about 80 of them. I want to get the expression and let it run from the Calculate Field.


texts: Shape area is not updated when importing shapefile from QGIS to Esri Geodatabase

texts:
I made some holes (difference) in a shapefile with Q-GIS 2.2.0, than I imported the shapefile in ArcGIS 10.0 and I exported the shapefile to a feature class in a File Geodatabase.

In the feature class the "shape area" doesn't show the right value.
I added a new field and I used "calculate geometry". ArcGIS calculated the right surface and it is different from the old one in "shape area".

Why ArcGIS 10 is not updating the value?

Is there a way to recalculate or edit the standard field "shape area"?


texts: Loading Shapefile in a standalone PyQGis App

texts:
I'd like to open a shapefile in a standalone PyQGis App. The App starts and runs without a problem but a get an false-Return in the isValid()-Function on the Layer and the Layer doesn't show up.

I tried this on two different shapefiles and did make sure that the filepath is correct.

Can anyone tell me is there is something wrong or missing in the code or how I can debug this problem any further?

import sys
from PyQt4 import QtGui
from qgis.core import *
from qgis.gui import *
from PyQt4.QtCore import *

app = QtGui.QApplication(sys.argv)

window = QtGui.QMainWindow()
window.setGeometry(550, 550, 500, 300)
window.setWindowTitle("Test Window")

window_frame = QtGui.QFrame(window)
window.setCentralWidget(window_frame)
frame_layout = QtGui.QGridLayout(window_frame)

canvas = QgsMapCanvas()
frame_layout.addWidget(canvas)

filepath = "/home/stefan/testlayer/test.shp"

layer = QgsVectorLayer(filepath, "testlayer", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(layer)
if not layer.isValid():
    print "Layer not valid"
canvas_layer = QgsMapCanvasLayer(layer)
canvas.setLayerSet([canvas_layer])
canvas.zoomToFullExtent()

window.show()

sys.exit(app.exec_())


texts:
We have a ArcGIS Network Analyst service which uses enterprise geodatabase (Oracle).

Frequently we receive updated data in shp format and we want to;

Overwrite existing data which we keep in database. (Delete old rows, update existing rows and add non-existing rows)
Reload - build network analyst service with minimum down time.

The method we use is really painful. We delete dataset, import data create network dataset.


texts: EPSG code lost when writing shapefile with sf

texts:
I've found that the EPSG can be recorded in the .prj file of a shapefile using the "AUTHORITY["EPSG","xxxx"] syntax. However, writing a shapefile with sf::st_write() doesn't seem to add that entry, and the EPSG code is lost when re-reading that shapefile. Patching the .prj file manually seems to work. Minimal reprex below.

What's the best way to record the EPSG in the .prj file?

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.1, proj.4 4.9.3
berne <- st_point(c(26e5, 12e5))
shp <- st_sf(name = "Berne", geom = st_sfc(list(berne), crs = 2056))
st_crs(shp)
#> Coordinate Reference System:
#>   EPSG: 2056
#>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"

st_write(shp, "berne.shp")
#> Writing layer `berne' to data source `berne.shp' using driver `ESRI Shapefile'
#> features:       1
#> fields:         1
#> geometry type:  Point

shp_read <- st_read("berne.shp", quiet = TRUE)
st_crs(shp_read)
#> Coordinate Reference System:
#>   No EPSG code
#>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs"

prj <- readLines("berne.prj", warn = FALSE)
new_prj <- gsub("]]$", '],AUTHORITY["EPSG","2056"]]', prj)
writeLines(new_prj, "berne.prj")

shp_read <- st_read("berne.shp", quiet = TRUE)
st_crs(shp_read)
#> Coordinate Reference System:
#>   EPSG: 2056
#>   proj4string: "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"


texts:
I'm quite new in either QGis-Plugin- and Qt-development.

I started to work on a Plugin, which displays a dialog (I used pluginBuilder and Qt-Designer) in which the user should select some options.
After the user has clicked Ok the action is performed. Everything works well so far.
Now I'd like to check the options, and if important information is missig I display a QMessageBox which tells the user what was wrong. So now to my Question: In this case the plugin-Dialog should not be closed but stay open for the user to fix the issue.
Can anybody give an example how to accomplish this?


texts:
I have a Leaflet application that uses circles. It works fine when I don't put the script in my HTML template (I use SingApp template), but when I integrate the script with the template, the circle does not appear on top of my map.

This is how it looks when I use it without the HTML template (SingApp template):

And here I put the script into the template:

I've checked the HTML element, there is a circle SVG element in the HTML document.

<svg class="leaflet-zoom-animated" width="1902" height="406" viewBox="-301 -64 1902 406" style="transform: translate3d(-301px, -64px, 0px);"><g><path stroke-linejoin="round" stroke-linecap="round" fill-rule="evenodd" stroke="red" stroke-opacity="0.5" stroke-width="5" fill="#f03" fill-opacity="0.3" class="leaflet-clickable" d="M434,140A52,52,0,1,1,433.9,140 z"></path></g></svg>

but still, circle does not appear on top of my map. I think there is conflict with template, but i dont know how to deal with it.


texts: Docker PostGIS .sql files not loading into database at build

texts:
I’ve recently started using docker and I’m trying to create a postgis docker image with raster and vector data inserted to a database at initialisation. I have no need to change or add to the data after this. I just want this to be available for requests from a website and geoserver.

I have made the image using the following Dockerfile:

    FROM mdillon/postgis

COPY ./nga_bth.sql /docker-entrypoint-initdb.d/
COPY ./nga_preg.sql /docker-entrypoint-initdb.d/
COPY ./nga_adm1.sql /docker-entrypoint-initdb.d/

EXPOSE 5432

The first 2 sql files are raster, the last is a vector.

I build the image using:

docker build -t dillon .

And I start a container with the following:

docker run -d -p 4000:5432 --name postgis_test -e POSTGRES_DB=gis -e POSTGRES_USER=postgres -it -e POSTGRES_PASSWORD=postgres dillon:latest

The container starts and then immediately disappears. I’ve tried this with different combinations of the above sql files. I checked the log on the container and it showed the following when I tried to load shapefiles:

server started
CREATE DATABASE

ALTER ROLE

/usr/local/bin/docker-entrypoint.sh: running /docker-entrypoint-initdb.d/nga_adm1.sql
SET
SET
BEGIN
CREATE TABLE
ALTER TABLE
2017-11-29 19:07:35.043 UTC [73] ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist at character 8
2017-11-29 19:07:35.043 UTC [73] HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
2017-11-29 19:07:35.043 UTC [73] STATEMENT:  SELECT AddGeometryColumn('public','nga_adm1','geom','4326','MULTIPOLYGON',2);
psql:/docker-entrypoint-initdb.d/nga_adm1.sql:18: ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist
LINE 1: SELECT AddGeometryColumn('public','nga_adm1','geom','4326','...
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

And the following when trying to load rasters:

server started
CREATE DATABASE

ALTER ROLE

/usr/local/bin/docker-entrypoint.sh: running /docker-entrypoint-initdb.d/nga_bth.sql
BEGIN
psql:/docker-entrypoint-initdb.d/nga_bth.sql:2: NOTICE:  table "nga_bth" does not exist, skipping
DROP TABLE
2017-11-29 19:03:39.043 UTC [73] ERROR:  type "raster" does not exist at character 66
2017-11-29 19:03:39.043 UTC [73] STATEMENT:  CREATE TABLE "public"."nga_bth" ("rid" serial PRIMARY KEY,"rast" raster,"filename" text);
psql:/docker-entrypoint-initdb.d/nga_bth.sql:3: ERROR:  type "raster" does not exist
LINE 1: ...ublic"."nga_bth" ("rid" serial PRIMARY KEY,"rast" raster,"fi…

I made the sql files from shp2pgsql and raster2pgsql for the vectors and rasters respectively. The above errors suggest to me that PostGIS isn’t enabled, which I don’t think is the case.

Could someone please tell me how I might create a postgis container with pre-inserted data?


texts: ArcCatalog says spatial references don't match between feature dataset and feature class with same coordinate system?

texts:
I have a feature dataset that is using the GCS_WGS_1984 as the geographic coordinate system. It includes a few feature classes.

This feature dataset is itself in a file geodatabase that contains a few more feature classes in the root of the gdb. One of these is a feature classes called "points" that is also using the GCS_WGS_1984 as its coordinate system.

I thought it would be trivial to use ArcCatalog to drag the feature class from the root of the gdb into the feature dataset since they have the same coordinate system. Unfortunately ArcGIS gives an error dialog that says:

Failed to paste points
The spatial references do not match

Why does this fail? If the GCS of the feature class and the feature dataset are the same, shouldn't it just work? I've verified this is the case via the dialogs as well as exporting the prj file for both the dataset and the class and using a diff tool to compare the two. They are identical.

Is the spatial reference of a feature class different than the coordinate system/projection?

Attemping to copy the feature class into the feature dataset by holding ctrl while dragging results in ArcCatalog crashing each and every time (shame on ESRI).

I guess I could try alternate means of moving feature classes. Use CopyFeatures into the dataset.  Project from the feature class into a new feature class within the feature dataset?


texts: Help converting a WKT projection to proj4

texts:
I have shapefile (from the local government gis office of a town in Massachusetts) that contains, in the corresponding .shp.xml file, the following WKT projection:

PROJCS["NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001",
    GEOGCS["GCS_North_American_1983",
        DATUM["D_North_American_1983",
            SPHEROID["GRS_1980",6378137.0,298.257222101]],
        PRIMEM["Greenwich",0.0],
        UNIT["Degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic"],
    PARAMETER["False_Easting",200000.0],
    PARAMETER["False_Northing",750000.0],
    PARAMETER["Central_Meridian",-71.5],
    PARAMETER["Standard_Parallel_1",41.71666666666667],
    PARAMETER["Standard_Parallel_2",42.68333333333333],
    PARAMETER["Latitude_Of_Origin",41.0],
    UNIT["Meter",1.0]]

Loading the shapefile into QGIS results in a null CRS ("This layer appears to have no projection specification...").  I've tried converting the above projection into proj4 syntax, but I'm running into some problems caused by the fact that I have no idea what I am doing.

The above projection is awfully close to the standard projection used by MassGIS (epsg:26986), although "Standard_Parallel_1" and "Standard_Parallel_2" are swapped. Trying to convert this to proj4 syntax using the Python osgeo.osr module results in:

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

I can replace Lambert_Conformal_Conic with Lambert_Conformal_Conic_2SP (as suggested elsewhere), which results in a successful export in proj4 syntax, but (a) it's not clear to me that's the right thing to do, and in any case (b) the resulting projection still doesn't place the features where they belong.

How do I get this data in the right place?

Edit: I have tried just using epsg:29686 as well, which didn't work either. The shapefile in question is here.


texts: Missing toolbars from ArcGIS Desktop on Citrix?

texts:
I was going to georeference some images using ArcGIS Desktop 9.3 on Citrix.

I turned on my "Georeferencing" toolbar and it was nowhere to be found. At first I thought I hadn't had enough coffee and it really was right in front of me. But it wasn't.

I've looked underneath other toolbars. I've changed the screen resolution to make sure that I didn't have it moved out off screen. I restarted my machine. I tried resetting the toolbar (custom-->toolbars-->reset). I deleted my normal.mxt, I created new mxds, I tried other mxd's. I deleted my whole ESRI user profile. I also have found that no toolbars, other than the defaults, show up when I click on them.

I've had other people look over my shoulder just to make sure I haven't gone crazy, that they really don't show up.

Has anyone else experienced this? What the heck is going on?


texts: Add and Calculate Field in Table with ArcPy?

texts:
I have written the following script to iterate it as a tool for all my tables with ModelBuilder, but i have a problem wth my codeblock. I want to add a field FID to all my tables and to add the values 0,1,2,3,4... for all rows of table. Each table has different number of rows. I am using ArcGIS 10.2.2 for Desktop.

The code i have written is this.

# Import system modules
import arcpy
import os
from arcpy import env
arcpy.overwriteoutput = True

# Set environment settings
inWorkspace = arcpy.GetParameterAsText(0)
inTable = arcpy.GetParameterAsText(1)
outTable = arcpy.GetParameterAsText(2)
data_type = ""

# Set local variables
# fieldName1
# fieldName2
fieldName1 = "FID_1"
fieldName2 = "Same"

arcpy.AddField_management(inTable, fieldName1, "DOUBLE", "", "", "", "", "NON_NULLABLE")
arcpy.AddField_management(inTable, fieldName2, "DOUBLE", "", "", "", "", "NON_NULLABLE")

arcpy.Copy_management(inTable, outTable, data_type)

# Calculate Field FID_1
fieldname = fieldName1
expression = "getCalc(!FID_1!)"
codeblock = """def getCalc(FID_1):
    FID_1 = 1
    for row in FID_1:
       FID_1 = FID_1 + 1
       return FID_1"""
# Execute CalculateField
arcpy.CalculateField_management(inTable, fieldname, expression, "PYTHON_9.3", codeblock)

I am newbie in Python.


texts: Time Manager in QGIS to perform spatio temporal analysis

texts:
I am working on identification and monitoring of pipeline oil spill using remote sensing and QGIS. (I am a new user.) I wish to know how to use the time manager in QGIS to monitor the incident events.

The objectives of my project are application of geo-informatics and remote sensing capabilities to identify incident cases in the hotspots where oil spillage is prevalent to perform spatiotemporal analysis in monitoring events.

My data was imported from excel sheet (raw data) and saved as CSV. Essentially the attributes of my data reflect the following features oil spill ID, Incident date, report date, estimated quantity, causes of spill and site location. I also have the coordinates of all the events. However, there is no information about the time.  Is this not a limitation to time manager plugin?

The layers of the spill events with potential IDs, dates and different quantities can be critically examined by zooming in to see the damage done by the oil spills to vegetation, water bodies and lands. This gives opportunity to look at the antecedents how the affected areas have been in the past years. It is like using the time slider on Google earth which give insight into the areas affected in the preceding years. This is what I think I may use time manager to perform on QGIS. But the time of events were not recorded in the data which may be difficult to perform using time manager. But as you mentioned that Crayfish focuses on visualizing which may explain my monitoring strategy, perhaps I can watch the video and see how it works. So essentially, it is to monitor the spill events from the perspective of looking at the same spot/affected area over a given period of time/year. Maybe I succeed in making my question clearer.


texts:
I have 120 rasters( jpg format ) and i want to convert them to tiff format to folder called
 "OtherFormat" in env.workspace with arcpy. In my code i try to convert all jpg file I read http://resources.arcgis.com/en/help/main/10.2/index.html#//001200000032000000 and understand that i can do it if i will write the exact name of the raster. But i have 120 rasters so i tried to use asterisk:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.workspace = r"C:\Project\out"
arcpy.RasterToOtherFormat_conversion("*.jpg","OtherFormat","TIFF")
print 'converted'

but get en error:

ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000732: Input Rasters: Dataset *.jpg does not exist or is not supported Failed to execute (RasterToOtherFormat).


texts: How to Add a single feature ( features[0] ) from Geojson ( without jquery )

texts:
var features = new ol.format.GeoJSON().readFeatures(geojsonObject,{featureProjection:'EPSG:4326'});
   console.log("I_got_features[0]_as_<point>_:"+features[0].getGeometry().getType());
 console.log("features.length_as_<132>:"+features.length);
/*
var featurething = new ol.Feature({name: "point",geometry: features[0]});
*/
 var vectorSource = new ol.source.Vector({
 //for_all_feature-----But i want Only 1 feature
 //features: features
 //so i tried all this below but failed
 //features: featurething
 //features: features[0].getGeometry()
   features: features[0]
 });
 var vectorLayer = new ol.layer.Vector({
 source: vectorSource,
 });
 map.addLayer(vectorLayer);

How can i convert my features[0] as a features......


texts: Is there a Postgis SQL Coding Style Guidelines?

texts:
Most programming languages have Coding styles. Some have defined standards, other have just some code of conduct that the majority of programmers follow.

Standard coding styles are have obvious benefits, in gis.stackexchange itself, it would help people better understand both questions and answers.

From some "googling", there are some reference to Coding Styles for SQL, although they are not consensual:

https://stackoverflow.com/questions/118288/sql-coding-style-guide
https://stackoverflow.com/questions/522356/what-sql-coding-standard-do-you-follow
http://www.simple-talk.com/sql/t-sql-programming/sql-code-layout-and-beautification/

While this SQL coding styles can be used for Postgis, my question is if there is a Postgis\SQL coding style to handle its functions, or at least some common practices.


texts: How to convert self intersecting polygon to multipolygon?

texts:
Working with the union and polygonize methods in GDAL/OGR, I often end up with self intersecting polygons. The case is one polygon which essentially consists of two polygons:

In the above example, the polygon has two vertices at the intersection point, but whether the lines meet and goes seperate ways at the point or crosses each other, is unknown.

Are there methods for splitting a self intersecting polygon like above, into two polygons or a multipolygon?
Alternatively are there any existing algorithms I could implement my self (I would rather not invent another wheel)?


texts: Switch line direction for PostGIS-lines in QGIS or Pgadmin

texts:
I am trying to switch the line direction in a PostGIS-layer. The postgis-geometry column is e.g. "LINESTRING(0 0, 1 1, 2 1, 2 2)".

The reason is that my symbology is direction-dependent, but I drew some lines in the wrong direction.

I have tried using my usual approach that seems to only work for shapefiles. It is the code by Nathan W. that is in the answer that is linked below.

Link to similar question:

How can I switch line direction in QGIS?


texts: Generating dartboard shapefiles using ArcGIS Desktop

texts:
I am attempting to generate a polygon 'dartboard' shapefile. Radially, there should be a circle from 0 to 20 km, an annulus from 20 km to 30 km, and an annulus from 30 km to 80 km. These three concentric shapes are divided into 8 segments to yield 24 sectors.

Within ArcGIS Desktop , within a new shapefile, I have generated three overlapping circles, with the same origin. However, I do not know how to take these multiple features and convert them into a single feature set (e.g. separate outer ring, middle ring, and inner circle).

I have attempted to use the Overlay tools in ArcToolbox but to no effect.

Equally, I have no idea, once that is done, how to divide this shape up into 45 degree sectors.


texts: How to connect points to only nearby points

texts: I have points for addresses and points for facades of the buildings in which the address lies. I dont have a buildings polygon shapefile. How can I connect each address point with only the facade points of that building? An example is attached in the image. Any ideas?


texts: Preserve column order of GeoPandas file read

texts:
I have a shapefile I'm reading in using GeoPandas and then I do some geometry conversion of the data but in the process it unfortunately sorts the columns alphabetically.  I need to preserve original order of the columns which has over 200 columns.

Here is the start of the code I'm using.  I imagine I have to somehow store the original column order and apply it to the GeoPandas dataframe.

import geopandas as gpd

shpfile = r"C:\data\myfile.shp"
poly = gpd.read_file(shpfile)
poly.head()


texts:
I'm trying to use the "apply" function with the R "eleavatr" package. I'm getting an error message but I'm not sure if the problem is in my understanding of the apply function (perhaps the root of the problem?) or a wholly different problem. The function works if you just use it alone ie get_elev(temp_df) but using it in apply causes problem. This is my sample data & code:

> temp_df
           x        y
41 -74.67020 41.36181
42 -74.66993 41.36208
43 -74.66966 41.36235
44 -74.66939 41.36262
45 -74.66912 41.36289
46 -74.66885 41.36316
47 -74.66858 41.36343
48 -74.66831 41.36370
49 -74.66804 41.36397
50 -74.66777 41.36424

get_elev<-function(ll){
  temp_el<-elevatr::get_elev_point(locations = ll, src = "epqs", units="meters",
         prj = ll_prj)
  return(temp_el)
}

sapply(temp_df, FUN=get_elev) #this should take the lat/lon coordinates of each row to search for elevation.

Error message:

 Error in if (attributes(class(locations)) %in% c("raster", "sp")) { :
argument is of length zero


texts: How do I change my "no data" values to get "true" values in raster calculator in Qgis

texts:
I have two rasters and I want to subtract the values from one of the layers.
The first raster is a altitude raster and the second raster is depth values for a lake.
I would like to subtract the depth values to get the altitude of the bedrock under the lake.
When I do this in the raster calculator (raster 1 - raster 2) the resulting raster is all wrong and have extreme negative values (-1,7 x10^35).
It has been pointed out to me that the "no data" values for one of the rasters are within the size of my  negative output.

How do I change the "no data" values to get the result I'm after?


texts: QGIS Raster Calculator returning only nodata

texts:
I have two perfectly aligned raster layers, both with a single band with pixels composed solely of 0s and 1s. I want to create a new raster layer also composed solely of 0s and 1s, where the 1s only appear where both of the input layers are 1. (I.e. I want to create the intersection of the two layers.)
I have tried using the QGIS Raster Calculator in various ways that I think should be equivalent, e.g.
("layerA@1" > 0) AND ("layerB@1" >0)

or
(("layerA@1" + "layerB@1") > 1) * 1

Whatever I do, I get an output layer where every pixel is nodata.
I also tested just:
"layerA@1" + "layerB@1"

And that also gave me an output of entirely nodatas (rather than, as I might have expected, a mix of 0s, 1s, and 2s).
What am I doing wrong? It's my first time using the raster calculator, so I imagine I'm getting the syntax wrong somehow. Thanks in advance.
EDIT: The two layer files in question can be downloaded here: https://gofile.io/d/XfXfop. In case you need to know, the coordinate system is EPSG:3857.


texts:
My map file is approximately one pixel = 10 meters.

I want to keep that relationship when I export the map as an image.

So I want to export it with original resolution.

How do I achieve that in QGIS 3.2.3?


texts:
Can anybody guide me with this error message from the FME workbench when I try to read a geotiff file and write it as a ECW file

  ERROR MESSAGE:  ECW writer: Unsupported band interpretation 'REAL32'.
  Please convert the interpretation to one of the following valid types:
  UINT8, GRAY8, GREEN8, RED8, ALPHA8, BLUE8 A fatal error has occurred.
  Check the logfile above for details A fatal error has occurred. Check
  the logfile above for details

How do I move on, can anybody link to some documentation eventually


texts: Check if layer is a marker when using Leaflet's eachLayer

texts:
I would like to use Leaflet's map.eachLayer(...) to do something to every marker in the map. I am not sure how to detect markers, however.

Is there a way to find out of what class the current layer is an instance of?

I guess I could add a custom property to each marker that I create, such as:

L.marker([40.49083 , -108.95966], {layerType: 'marker'}).addTo(map);

And then check if that property exists... but maybe there is a more elegant way?


texts: Unable to import arcpy into Spyder with ArcGIS Pro 2.7

texts:
I am trying to import arcpy into Spyder. I saw  ArcPy import arcgiscripting  relating to this, but I'm using ArcGIS Pro 2.7, not ArcMap. My original Python is 3.8 but for the environment that I created to run arcpy in Spyder, I installed Python 3.6 (following an error I got previously). This is similar to what was suggested at this answer to Importing arcpy in Spyder, but when I install Spyder on the newly-created anaconda environment, I try to import arcpy into Spyder and get the following error:
 File "C:\Users\anner\AppData\Local\Programs\ArcGIS\Pro\Resources\ArcPy\arcpy\geoprocessing\_base.py", line 14, in <module>
    import arcgisscripting

ImportError: DLL load failed: The specified procedure could not be found.

I've tried updating the Python interpreter in Spyder to the python.exe that was automatically installed with ArcGIS but got this message

Any ideas?


texts: Why do quickmapservices basemap labels shrink when printed?

texts:
I am using QGIS with quickmapservices for my basemaps. In the composer the image looks good, with labels appropriately sized.

When I go to print the PDF however, the scales on everything is much smaller, yielding unreadable labels for street names, etc.

Is there any way at all I can mitigate this, without taking a screenshot of each atlassed image and pasting it into my final result?


texts:
E.g. in SAGA GIS it is possible to set with one click the number of CPU cores to use during processing.

Is this also possible in Grass GIS?


texts: Simplifying TopoJSON causes leak of Natural Earth in d3

texts:
I'm working with topoJSON and D3 to display Natural Earth data on a globe. Obviously, we need to simplify our geometries while still keeping borders consistent. We've created a script as based on Mbostock's Map Tutorial. Here's the link to our gist.

Once we run the script, it outputs all the countries as .json files.

However, when we load the USA, we encounter a problem. When we mouse over the USA, the entire globe is highlighted.

United States without mouse hover.

United states on mouse hover

Is there something in our script that would cause this? We need to simplify our geometries for performance purposes, and we need to preserve boarders. TopoJSON was the clear choice, but we could be calling it incorrectly.


texts: Handling multiple map frames with different layers in one print layout?

texts:
Often I need two (or more) map frames in map layout to present different data next to each other. It is easy to make multiple map frames in QGIS Print Composer. BUT...

It is in the nature of my work that very often I need to update some map frames because data I am presenting in map changed or I need to add/remove layer, change scale, projection etc.

So if I create one map frame, lock it (Lock layers for map item), create second map frame with different layers, scale,projection etc. and than realize I need to change something in first frame, I can't find a simple way to do it...

To simply do that in Mapinfo, there are more Map Windows, in ArcGIS there are multiple Data Frames. What is in QGIS to achieve that?


texts: QGIS Spatial join using subset of features?

texts:
I have a road network shapefile, as well as a postgis layer of different areas in a city.  The areas file contains one big polygon containing everything (the service area), and then a series of smaller polygons that split up the city by neighborhood.  I want to assign info from the smaller neighborhood polygons to the road network based on location.  If I do a simple spatial join, then the road network will get info from the larger service area, which I do not want.  Is there an easy way to tell QGIS to exclude the service area from the join?

I realize I could export only the neighborhood polygons to a new file, then do a spatial join, but I would like to avoid creating a new file every time I have to do this (which will be weekly).

I am using QGIS 2.8.6 and Python 2.7.  A script would be ideal, but I can use the graphical interface if necessary

Edit: I tried using the intersect tool, but it is going really slowly for my network layer of about 180,000 features (on track to take like 20 minutes).  In the future I will need to do this for layers with over a million features, so this won't be acceptable


texts: Sorting feature class to calculate sequential ID field using ArcGIS Field Calculator?

texts:
I am in ArcMap and I want to sort the records based on a field and then do a Calculate Field with an auto-incrementing value with the sort in-place. Any ideas?

It looks like I could get an update cursor on the feature class then do an auto-increment, but I would still like to see if there is a way to do this in the Field Calculator in ArcMap.


texts: How do I successfully install LeCos plug in for QGIS2.0 on a Mac?

texts:
I have been trying without success to install the LeCos plug in on QGIS2.0 on my Mac. I originally tried it on QGIS1.8 but was unable to install the dependent packages (PIL,scipy and numpy). I updated to QGIS 2.20 and successfully installed the packages, however I am still unable to install LeCos. THis is the message that comes up:

The packages are definitely installed because I checked by doing this in python:

Has anyone else installed it successfully?


texts: Errors when creating a TextElement from ArcObjects

texts:
I am attempting to add some text in the layout view of ArcMap.   I’m using the following line for the TextElement that I’ve found online from numerous places including ArcObjects Library Reference.  Being new to ArcObjects, I'm not sure if I'm doing something wrong or if there is something else going on that I'm unaware of.

ITextElement textElement = new TextElementClass();

I’m receiving the following errors:

The type TextElementClass has no constructors defined
Interop type TextElementClass cannot be embedded.

Full Code:

protected override void OnClick()
    {
        IMxDocument mxDoc = ArcMap.Application.Document as IMxDocument;
        IMap map = mxDoc.FocusMap;
        IActiveView activeView = map as IActiveView;
        IPageLayout pageLayout = activeView as IPageLayout;

        if (activeView is IPageLayout)
        {
            IGraphicsContainer graphContainer = pageLayout as IGraphicsContainer;
            ITextElement textElement = new TextElementClass();
            textElement.Text = "Test";
            IElement element = textElement as IElement;
        }
        else
        {
            MessageBox.Show("This tool only functions in layout view");
        }
        activeView.Refresh();
    }


texts: Using QGIS2Leaf plugin to export into leaflet but output results into browser crash?

texts:
Amateur GIS technician here, Just tried to export my converted cad cities pipelines into a shape files and then export it into leaflet using QGIS2Leaf but the browser crashes before the displaying the output.

Is it possible that what I'm trying to do is too much for the broswer? The size of the shape files is 110mb.


texts: SQLite layer from AutoCAD MAP in QGIS editing problem

texts:
SQLite files generated with FDO (AutoCAD MAP) have FGF geometry instead of WKB.  FGF is an extended version of the Well Known Binary (WKB) format. QGIS can read these layers but they can't be edited, so they must be converted in this way:

create SQLite in MAP;
open SQLite in QGIS;
export (convert to WKB) SQLite in QGIS;
after that this SQLite works in both programs but its schema must not be changed in MAP or conversion must be done again- not very practical;

My question: is it possible to add to QGIS ability to edit SQLite layers with FGF geometry?


texts:
I have some code here that reprojects raster files. However, once in a while, I get these errors and I am not sure how to interpret the error messages.

Here is my code:

library(rgdal)
library(raster)

##PROJECT RASTER
# Create RasterLayer object
r <- stack('Fil1_0.tif')
# Define the Proj.4 spatial reference
# http://spatialreference.org/ref/epsg/26915/proj4/
sr <- "+proj=utm +zone=31 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

# Project Raster
projected_raster <- projectRaster(r, crs = sr)

And here is the error message I get:

Error in projectExtent(from, projto) : cannot do this transformation
In addition: Warning message:
In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[,  :
  25028 projected point(s) not finite

This is a single band raster so I don't understand how there can be a problem with the extent.


texts: Open multiple QGIS projects via PyQGIS simultaneously in different threads

texts:
I'm writing a standalone PyQGIS application. When the method qgis() is called in a new thread, it's supposed to open a QGIS project (with one layer in it) and close it again:

from concurrent.futures.thread import ThreadPoolExecutor

from qgis.core import *

qgs = QgsApplication([], False)
QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.4", True)

def qgis(number):
    qgs.initQgis()
    print('{} started'.format(number))

    QgsProject.instance().read('D:/aGIS/pyqgis/multitest/testprojekt.qgz')

    qgs.exitQgis()
    print('{} stopped'.format(number))

with ThreadPoolExecutor(max_workers=3) as executor:
    executor.submit(qgis, '1')
    executor.submit(qgis, '2')
    executor.submit(qgis, '3')

Unfortunately, this is what the console says:

  Type conversion already registered from type QPair to type
  QtMetaTypePrivate::QPairVariantInterfaceImpl

  3 started

  QObject::setParent: Cannot set parent, new parent is in a different thread

  QObject::connect: Cannot queue arguments of type 'QList'

  (Make sure 'QList' is registered using qRegisterMetaType().)

  QObject::connect: Cannot queue arguments of type 'QList'

  (Make sure 'QList' is registered using qRegisterMetaType().)

  QObject::connect: Cannot queue arguments of type 'QDomDocument'

  (Make sure 'QDomDocument' is registered using qRegisterMetaType().)

  QObject::connect: Cannot queue arguments of type 'QList'

  (Make sure 'QList' is registered using qRegisterMetaType().)

  1 stopped

And then Python crashes. I made some experiments with moving

qgs = QgsApplication([], False)
QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.4", True)

into the function, resulting in slighty different errors but nevertheless it's not working. Am I making a mistake or is it just not possible loading multiple project instances at once?


texts:
I have a list with 4 raster files and another list with 4 vector files. I want to relate each of the raster file 1 file list 1 vector list.
I try to do and I related file 1 raster to vector file 1,2,3,4. How I can do to relate only 1 with 1, 2 with two and so on?

ListaRaster = arcpy.ListRasters("", "img")
1_23165_1997204.img
2_23065_1997197.img
3_23166_1997204.img
4_23066_1997197.img

ListaVector = arcpy.ListFeatureClasses("r*")
rec_1.shp
rec_2.shp
rec_3.shp
rec_4.shp

for r in (ListaRaster):
  for v in (ListaVector):
    print ('Raster', r, 'se relaciona con vector', v)
    ('Raster', u'1_23165_1997204.img', 'se relaciona con vector', u'rec_1.shp')
    ('Raster', u'1_23165_1997204.img', 'se relaciona con vector', u'rec_2.shp')
    ('Raster', u'1_23165_1997204.img', 'se relaciona con vector', u'rec_3.shp')
    ('Raster', u'1_23165_1997204.img', 'se relaciona con vector', u'rec_4.shp')
    ('Raster', u'2_23065_1997197.img', 'se relaciona con vector', u'rec_1.shp')
    ('Raster', u'2_23065_1997197.img', 'se relaciona con vector', u'rec_2.shp')
    ('Raster', u'2_23065_1997197.img', 'se relaciona con vector', u'rec_3.shp')
    ('Raster', u'2_23065_1997197.img', 'se relaciona con vector', u'rec_4.shp')
    ('Raster', u'3_23166_1997204.img', 'se relaciona con vector', u'rec_1.shp')
    ('Raster', u'3_23166_1997204.img', 'se relaciona con vector', u'rec_2.shp')
    ('Raster', u'3_23166_1997204.img', 'se relaciona con vector', u'rec_3.shp')
    ('Raster', u'3_23166_1997204.img', 'se relaciona con vector', u'rec_4.shp')
      .......

And I want the result to be:

('Raster', u'1_23165_1997204.img', 'se relaciona con vector', u'rec_1.shp')
('Raster', u'2_23065_1997197.img', 'se relaciona con vector', u'rec_2.shp')
('Raster', u'3_23166_1997204.img', 'se relaciona con vector', u'rec_3.shp')
('Raster', u'4_23066_1997197.img', 'se relaciona con vector', u'rec_4.shp')

How can I do it?


texts:
I'm tryin' to use the mask plugin associated with atlas composer (QGIS 2.7 dev).
I'd like to have in the legend of map composer only the features that belong to mask instance (inside the mask).
if I select the "green filter" for the Legend elements, there are visible in map composer also the features outside the mask.

Is it possible?


texts: Assigning unique Values to adjacent raster cells

texts:
I have a raster file like the one below. All green cells have the same value. The white space between the green cells is no data. I want to combine adjacent raster cells to one "count" and give them a unique value. Is this possible?  I´m using ArcMap.

This is what I want to do:


texts: How do I remove the orthoimage layer from the new USGS GeoPDF topo maps?

texts: I would like to remove the orthoimage layer from a USGS Topo GeoPDF map. I have Adobe Acrobat Pro 9, but can't find a way to delete a layer. I've tried setting the layer to "Never export" then exporting the pdf, but the image is still there (and in the process I have somehow deleted the geo-referencing info).


texts: Query for numeric field against text field that contains a numeric value

texts: In my database, I have one field, Service_Day, formatted as numeric, which contains a single digit number 1 through 7 (for days of week). I want to compare it against another field, called Region, formatted as Text/String which also contains a single alphanumeric character 1 through 7 (again representing days of the week). It seems that I should be able to query the records where "Service_Day" <> "Region", by using the proper syntax. I thought AsString or AsNumber or Convert might be appropriate, but can't figure it out. Any ideas, other than creating a new field and formatting it as either text or numeric and then copying whichever data to that field?


texts: Access Sextante (processing) in standalone QGIS app?

texts:
Sextante is addressable via python using the QGIS python console (unfortunately not from outside QGIS). I´m writing a standalone app and want to include Sextante geoprocessing functionality. In QGIS python console this is working perfectly. Is it possible to address the console from an outside python script (standalone app)?

My code so far:

import sys, os
from PyQt4.QtGui import QApplication, QAction, QMainWindow
from PyQt4.QtCore import SIGNAL, Qt, QString
app = QApplication(sys.argv)
from qgis.core import *
import qgis.utils
import qgis.gui

QgsApplication.setPrefixPath("C:/OSGeo4W/apps/qgis")
QgsApplication.initQgis()

sys.path.append("C:/Users/.../.qgis/python/plugins")
sys.path.append("C:/OSGeo4W/apps/qgis/python/plugins")
sys.path.append("C:/Users/.../.qgis/python/plugins/sextante")

import sextante
sextante.core.Sextante.Sextante.initialize()

map = qgis.gui.QgsMapCanvas()
layer = QgsVectorLayer("D:/Python_Test/a.shp",'a','ogr')
QgsMapLayerRegistry.instance().addMapLayer(layer)
map.setExtent(layer.extent())
map.setLayerSet( [ qgis.gui.QgsMapCanvasLayer(layer) ] )

output = "D:/Python_Test/b.shp"

sextante.runalg("qgis:convexhull",layer,None,None,output)

QgsApplication.exitQgis()


texts: How to use Time as a parameter in Network Analyst Service Area and Closest Facility

texts:
I have a road network downloaded from Ordnance Survey (Meridian 2, with A roads, B roads and Minor roads Merged into one shapefile and then built into a Network Dataset for use in ArcGIS' Network Analyst) and would like to use it to determine Service Areas for 4 hospitals (preferably with drive times of 30 minutes and 60 minutes) and Closest Facilities (same 4 hospitals and directions from 539 points).

I've used the Network Dataset to find the Service Area and the Closest Facility but the only parameter I can use is distance (metres).

Question: How do I set up the network to use drive time, not distance, and work out Service Area by, say, 30 minutes?

When I try to create the Network Dataset from the downloaded road network, I try to add Drive Time as an Attribute but I can't get the Network Dataset Wizard to accept the evaluators.

I've tried to follow the ESRI help pages and I've Googled around but I'm lost. I'm not sure what the time setup should be or if it's something I can just do or something I need to set up from the start.


texts: Model or script to calculate distance from each point in a multipoint file to the center of a matching polygon in ArcGIS Pro

texts:
I am working in ArcGIS Pro. I have a series of multipoint features that represent release locations for individual animals. I have a series of these, for multiple study sites. At each study site - I would like to have a method to calculate the distance between the release location of an individual animal, and the center of a polygon that represents it's home range during the study period.

I have a large number of individual polygon shapefiles calculated using R that represent the homeranges, their file names begin w a letter for the site and a two digit number for the animal ID. Because these were calculated in R - the attribute table for these polygons does not have a unique identifier that I can easily use to match this polygon to the corresponding release location for that animal.

Ideally - I would like a script to import the polygon shapefile, find it's center (convert polygon to point) and then calculate the distance from that point to the release location of the same animal - which is part of the release site shapefile. I know in R it is possible to use something like list.files() to get the needed data from the file name, but I'm not sure if that is possible here with ArcPY?

So, basic steps are select a release site - then measure distance to the polygon centroid that matches the identity of that individual, in this case C01.

If there's no elegant solution to this because my files are messy I suppose I can go the clunky way and measure all of these manually, but I would rather not and I am trying to learn more ways to make my workflow more efficient.

Here's a photo example of that workflow:


texts: QGIS modeler: Error when extracting by non existing attributes. How to avoid displaying of error?

texts:
I have a problem with the QGIS modeler:

The idea is to extract attributes ("DN" between 40 and 44) out of one layer.

I want to use the outcoming layer for further algorithms.

I calculated with the field calculator a new field named "selection2" using the following formula:

"DN"  >= 40 AND  "DN"  <= 44

Then I extracted by attribute for selection2=1.

The attribute table is empty, as in the original layer exists no DN between 40 and 44.

The problem is that the fied calculator is not able to use empty attribute tables for further algorithms (here: "extract by location"). Therefore this error is displayed:

Error executing algorithm Extract by location float division by zero See log for more details

It would be good to avoid the error and have a layer algorithms can work with.

Is there a algorithm in the modeler which functions like an if/else-function?

If there is none I have the following approach:
Fill in the attribute table any value, e.g. 9999, in order to not have an empty one. The question is then how to create this command?

As I do not have any experience with python, I cannot use it.


texts: "Unknown datum" opening GRS80 Transverse Mercator in R

texts:
I am working with a shapefile (available here as a 15MB ZIP) with a projection defined as follows:
GEOGCS["GCS_GRS_1980",DATUM["D_GRS_1980",SPHEROID["GRS_1980",6378137,298.2572221]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]]
As far as I can tell this is GRS 1980 Transverse Mercator.
I am able to open this shapefile using st_read("srprec_p04_state_shp") but the output includes the following:
Geodetic CRS:  Unknown datum based upon the GRS 1980 ellipsoid
Is this "unknown datum" a red flag, or is it OK?
I just want to avoid erroneous output in subsequent spatial operations.
Edit: I've been informed by the data provider that the correct datum is NAD83, so I'll be setting that manually after reading in the shapefile. I'm still curious, though: What is "D_GRS_1980"? My understanding is that GRS80 is just an ellipsoid/spheroid, not a complete CRS, as discussed here. And I can see that a NAD83 .prj file:
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
looks quite similar to the one that was packaged with the shapefile I'm working with, aside from the first argument in DATUM.


texts: Orfeo mshift segmentation asks for condition

texts:
I'm trying to run a segmentation using the mean shift. I leave all the parameters as default, but when I run it, it says "Wrong or missing parameter value: Condition".
Condition is not a parameter for the meanshift, but for the cc algorithm.

When running meanshift as batch processing, you can see that a parameter called "Condition" is requested.
The image is a small Landsat composite.


texts: Shapefile layer not showing in Windows/QGIS

texts:
I am trying to load a US county shapefile in QGIS (for Windows). It appears to load OK and a new layer appears in the project, but I cannot see the shapefile outlines. When I do the same thing with the same file on Mac, it works just fine. Any thoughts?
I am on QGIS 3.4.14(Madeira) on Mac and 3.14.1(Pi) on Windows. Link to shapefile:
https://www2.census.gov/geo/tiger/GENZ2018/shp/cb_2018_us_county_500k.zip


texts: Autocomplete polygon or autocomplete freehand in QGIS

texts:
Is there any similar function or plugins for QGIS as in ARCMAP for polygons called autocomplete polygon or autocomplete freehand. I am digitising polygons on a DEM as geomorphological research.  I need to digitise polygons next to each other with 0 overlap.

I'm a searching for this function in QGIS 3.4.9-Madeira on a Mac.


texts: Creating a dropdown widget in QGIS for QField project

texts: I want to collect data on crop species affected by diseases using QField. I want to create a widget so that I can have a dropdown list of the different crop species. In theory, if I were to export this project onto QField and collect data, I would select the crop species from the drop down list. Most of the answers I've found suggest changing the widget type to relation reference and link it to a CSV, but my data is already in the attribute table. Has anyone done this before? Which widget type should I use? Or any tutorials?


texts: Handle Multiple Layer Max Zoom Levels in Leaflet

texts:
Along with Mapbox layers, we are also adding layers from WMS and ArcGIS servers hosted by our customers to display map data relevant to them. Is there any way to figure out what the max zoom for a layer is through code so the maxNativeZoom can be set to that so they always display even at zoom levels higher than their max zoom?

Another thing I've noticed is that certain areas of Canada have satellite imagery up to zoom level 17 while others have it up to 18.  When moving to a specific location on the map, is it possible to figure out what the max zoom is in that area and then set the maxNativeZoom appropriately?


texts:
I am trying to create a database connection from ArcGIS desktop 10.1 to SDE Geodatabase on Oracle server with direct connection using Oracle easy connect syntax
So, from ArcCatalog I use "Add database connection"
then I chose Oracle in database platform

then in the instance I write down

myserver: port/serviceid

so I write

dbserer:1539/devDatabase

but it gives me ora-12154 TNS:could not resolve the connect identifier

I thought by using the easy connect I don't have to worry about TNS names


texts: Changing background color to transparent in GeoServer?

texts: How to set the background color to transparent. I want use style which describe in *.sld file for all my photo with tif format. How to do that?


texts: Exporting used shapefiles and saving them in separate folder using QGIS?

texts:
In my project, I have a lot of shapefiles with different information. I won’t need all of those shapefiles, but just a few.

I would like to know if there is a way to export all the used shapefiles in my project and save it in a separate folder? That way, I would just have to copy one folder, without having to select every shapefile one by one. Another advantage would be, that in case of a change in my folder, I would just have to change it one time and I could use my origin information in different projects as well.

I heard that this is possible in ArcGIS. I am not sure if I could express myself clearly, but I hope so.

I use QGIS 12.2.3.


texts: Finding the index of a field with its name using arcpy

texts:
I want to extract the index of a field with its name using arcpy. It seems what I have written is so dumb. Is there any simpler method with out using for-loop?

    import arcpy

    fieldList = arcpy.ListFields(self.address_ref_layer1)
    self.field_name = []
    for field in fieldList:
        self.field_name.append(field.name)
    #self.remove_pseudo_nodes()
    print self.field_name
    print self.field_name.index(u"NAMN1")

It could also be written:

    fieldList = arcpy.ListFields(self.address_ref_layer1
    name_index = [f.name for f in fieldList].index(u"NAMN1")

The ListFiled function returns list of field objects with many properties as well as name, but no function to return the index of the field which has a specific name?


texts: primary key in sql server uneditable in mapinfo attribute?

texts:
I trying create table relationship in sqlserver 2008 R2. Tables will have spatial table and non-spatial (tabular) .  I having trouble while selecting varchar field as primary key, this field cannon be editted on mapinfo attribute while I create new shape.
I tried use MI_PRINX as primary key but still success to get primary key field editable in mapinfo attribute editing.


texts: Selecting line features by topology in QGIS (select dangles)

texts:
I want to select line features by their relationship to their neighbouring line features.
To perform generalisation on the data set I need to select lines that:

Do not intersect with their end point with other features
OR
Do not intersect with their start point with other features
AND
are below a certain length (let's say 100)

I need to select dangles (line segments with one unconnected end (or start) point).

So far, I have tried (updated):

(intersects(start_point($geometry), difference(collect($geometry),$geometry)) AND  disjoint(end_point($geometry), difference(collect($geometry),$geometry))) OR (intersects( end_point( $geometry ), difference(collect($geometry),$geometry)) AND disjoint(start_point($geometry), difference(collect($geometry),$geometry)))

The output is not what I am looking for though.
Some dangles don't get selected, and some line features that aren't dangles are selected.

How can I do this?

Update regarding topology tools: The tools "Topology Checker" and "Check Geometries" seem to have problems finding real dangles. Both tools target end or start points of line features that are definitly connected to other features (vertexes and segments).


texts: How to use ModelBuilder to run tool extract by mask on many rasters in ArcGIS for Desktop?

texts: I'm trying to extract the same part of multiple rasters (Landsat scenes) using a polygon so I can do a time series.  I've been trying to do it in ArcMap 10.2 ModelBuilder but the following error "ERROR 000581: Invalid parameters. Failed to execute (Extract by Mask)" keeps occuring. I've attached a screen shot of my model.  Any ideas on what I'm doing wrong?


texts: Getting Projection Code (EPSG) from Coordinates and known location

texts: I have Coordinates (x=3381438.573, y=5794542.432) and a location (somewhere south of Berlin). Is there a way to find out the Projection code this data was defined with? A website would be nice.


texts: Most efficient way to check graph edge or node overlap with circle border on euclidean plane (2D space)

texts:
Suppose that there is a given graph, in terms of nodes with specific placement on the euclidean plane (2d space) (for ex. node A is at (7,14) 7 is x-axis and 14 is y-axis) and edges that specify which nodes are connected with which (weights or directed or not, are irrelevant).

Now suppose that we are given a circle (by specifying the center and the radius). What is the most efficient way in terms of complexity to find, if any, all parts of the graph (nodes or edges) that overlap with the circle border?

Is there a specific data structure that the graph must be saved as to enable quick indexing ?? A kind of binary search for 2d points for example (if a similar thing even exists) ???

Keep in mind that the reason i am pinpointing algorithmic efficiency is the fact that the input data sets are huge. Millions or maybe billions of nodes and edges. Subsequently if the approach is parallelizable too, it would really be of great help.

Update: I've found that 2 approaches may be feasible, r-trees and spatial hashing. I do not know though if r-tress support(and perform well on) equality queries (for the circle border), and if it is possible to incorporated the edges in a spatial hashing implementation.


texts:
I am creating a shapefile based on three columns in my CSV File - Latitude, Longitude, and Laps. For the most part, I understand how to loop through the data and create the many points based on the latitude and longitude column, but I am having issues with my Laps column. I need to be able to loop through the data but exclude blank columns. Below is a sample of my data from the CSV File. The row that needs to be excluded is the row with the total Lap time for that Lap (i.e. #Lap 0...)

Below is my code that I am using. When I run it, I get an Error 99999: Something unexpected has caused the tool to fail.
# This script reads GPS tracks in CSV format and
# writes the geometries (polyline) from the list of coordinate pairs

import csv
import arcpy

def addPolyline(fc, array, sr):
    polyline = arcpy.Polyline(array, sr)
    with arcpy.da.InsertCursor(fc, ("SHAPE@",)) as cursor:
      cursor.insertRow((polyline,))
      array.removeAll()

# setup workpace to save the file
arcpy.env.workspace = r"C:\\"

# setup the feature class and spatial reference
polylineFC = "carpath.shp"
sr = arcpy.SpatialReference("GCS_WGS_1984")
arcpy.CreateFeatureclass_management(arcpy.env.workspace, polylineFC, "POLYLINE", spatial_reference = sr)
spatialref = arcpy.Describe(polylineFC).spatialReference

# open the input file
wakefieldRaceWay = open (arcpy.env.workspace + "\WakefieldParkRaceway_20160421.csv", "r")
# setup csv reader and process the header
csvReader = csv.reader(wakefieldRaceWay)
header = next(csvReader)

latIndex = header.index("Latitude")
lonIndex = header.index("Longitude")
lapIndex = header.index("Lap")

# create an array
vertices = arcpy.Array()

# create a loop to iterate over lines but
# exclude lines that have zero data in them
for row in csvReader:
    isLap = row[lapIndex]
    if isLap == '\n':
        if vertices.count> 0 :
            addPolyline(polylineFC, vertices, spatialref)

    # Loop through the lines in the file and get each coordinate
    # put the coordinates into a tuple and add it to the list
    Latitude = row[latIndex]
    Longitude = row[lonIndex]
    vertex = arcpy.Point(Latitude, Longitude)
    vertices.add(vertex)

# add the last segment manually
addPolyline(polylineFC, vertices, spatialref)


texts: Measuring length from one polyline vertex to other polyline

texts: I have two polyline which represent the river bank of a particular section of a river drawn based on Satellite imagery at 1984 and 2019. As you can see in the picture, it is drawn in red line with yellow vertices and another is orange line with blue vertices. I want to measure perpendicular length from one polyline's vertices to other polyline i.e yellow vertices perpendicularly on orange polyline and blue vertices perpendicularly on red polyline and resuting length will be added to the attribute table of vertices layer showing the perpendicular line also. How can I proceed to do with QGIS 3.X?


texts: it is possible to use GeOnetwork java library for exposing a spatial resource(a server side spatial resource in a web app)?

texts:
we have a web app using spatial data, saved inside .
We would expose our spatial data to a csw service . A indipendent ensri geoportal data should map the resources those are present in our web app.
So i thing we would expose the spatial resources in a specific protocol.
GeoNetwork library can publish a csw resource using federation or is it better esposing wfs interface for then bindint it to csw service?


texts:
I'm trying to check if several gdbs exists, but I keep getting an error. I really don't understand why and I have tried to check the separately and that works.

sti = r"C:\\DTM"
#DB
if arcpy.Exists(sti + "DB.gdb"):
    print ("DB.gdb finnes")
else:
    print ("DB.gdb finnes ikke, oppretter GDB")
    arcpy.CreateFileGDB_management(sti, "DB.gdb")

#DTM
if arcpy.Exists(sti + "DTM.gdb"):
    print ("DTM.gdb finnes")
else:
    print ("DTM.gdb finnes ikke, oppretter GDB")
    arcpy.CreateFileGDB_management(sti, "DTM.gdb")

I get this error:

  File "<string>", line 8
if arcpy.Exists(sti + "DTM.gdb"):
 ^

SyntaxError: invalid syntax


texts: Creating SpatiaLite view (for use in QGIS) from non-spatial SQLite database?

texts:
I have a plain SQLite database that contains X, Y coordinates in a table.  Data are edited via an MS Access front-end using the ODBC driver and so I want to avoid storing them in a geometry blob that can only be edited using Spatialite.

I know that QGIS can plot tables of .csv files by specifying an X and Y column and that Arc can do a similar thing with plain SQLite tables.

How can I create a view that turns my X, Y columns into Spatialite geometry objects so that the resulting data can be viewed on a map in QGIS?

Example table structure

-- Location table (non-spatial)
CREATE TABLE location (
  id INTEGER PRIMARY KEY,
  name TEXT,
  x INT,
  y INT,
  epsg_code INT
);

INSERT INTO location VALUES
  (1, 'North', -3, 55, 4326),
  (2, 'South', 0, 52, 4326);

-- Sample data to be plotted on map
CREATE TABLE sample (
  sample_id INTEGER PRIMARY KEY,
  description TEXT,
  location INT
);

INSERT INTO sample VALUES
  (1, 'big', 1),
  (2, 'bigger', 2),
  (3, 'biggest', 2);

Example view

I would like to create a view that will let me plot the following attributes on a map:

SELECT
  s.sample_id,
  s.description,
  loc.name
FROM
  sample AS s LEFT JOIN location as loc
  ON s.location = loc.id;

Similar questions

I am posting this question to answer it myself. I found parts of the answer in the following questions, but wanted it to be documented in a single place.

How to make a Spatialite View (with geometries) for QGIS based on two table
Can QGIS Read Spatialite Views?
Spatialite view appears 'filtered' in attribute table, won't load in canvas
Creating SpatiaLite View in QGIS?


texts: Finding the highest value contiguous pixels in a raster

texts:
I have a raster of ~60,000 pixels, and each pixel has an integer value. I would like to find the ten contiguous pixels which sum to the greatest value. Ideally this would be achievable in QGIS.
I suspect that the number of possibilities means finding an exact solution may be impossible, and this will end up an exercise in image segmentation/object recognition to group high value pixels together
Can someone with more experience suggest whether this is possible or not?


texts: Columns and their types for shapefile in OGR/python

texts:
I am trying to list all the columns (or features) and their types for a shapefile.
shp2pgsql shows following schema:

CREATE TABLE "liberia" (gid serial,
"local_name" varchar(100),
"eng_name" varchar(100),
"cntry_name" varchar(75),
"type" varchar(75),
"cntry_code" varchar(10),
"continent" varchar(50),
"iso_code" varchar(10),
"volume" int4,
"id_7_1" varchar(100),
"id" numeric(10,0));
ALTER TABLE "liberia" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('','liberia','geom','5578','MULTIPOLYGON',2);

The following code snippet tries to achieve that:

 from osgeo import ogr
 reader = ogr.Open("file.shp")
 x = reader.GetLayer(0)
 fields  = x.schema
 for l in ll:
    print l.GetName()
    print l.GetType()

Which prints output as follows:

LOCAL_NAME
4
ENG_NAME
4
CNTRY_NAME
4
...

I am facing two problems. One the geom column/feature is missing from the list and second deciphering type id to python or postgres type.


texts:
I'm trying to create random points using slpancs library in R. But I have this problem:

random_points<- csr(poly = "zona", 1000)
Error in poly[, 1] : incorrect number of dimensions


texts: Arcpy if statement to act on layers at the top level of hierarchy

texts:
I am about to save copies of all of the layers in my mxd. I only want to save the top level group layers though. I have some group layers within group layers - I dont want to save them individually. So far I can export all group layers like this:

mxd = arcpy.mapping.MapDocument("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
 if lyr.isGroupLayer == True:
     lyr.saveACopy(lyr.name.replace(":"," -").replace("&","and").replace(">","") + ".lyr")

What I'm looking for is either a property that says what level the layer is on so I could say if layer level = top level...  (instread of if lyr.isGroupLayer == True:) or a property that says that the layer is within a group layer so I could say if lyr.isInGroupLayer == False:.

I suppose this could be in my for loop or a change to the mapping.ListLayers

Any ideas?


texts: Line still appearing but not existing in the metadata of a shapefile

texts:
For a project, I had to merge the cells of a grid to upscale it. It was possible to merge all I wanted but there is kind of a displaying problem: there is this line which appears visualy but that is not existing in the metadata or anywhere (see below, as you can see there are 2 undesired line, the first on the column of the selected cell and the other one more on the right)

What is weird is that when I select the cell, it is selected as one cell as if this line doesn't exist. Same in the metadata, there is no information about this line.
It is not a problem for my analysis as this line doesn't really exist but I would like to display it properly so would you have any suggestion on how to erase (even manually), this line?
Cheers


texts: How to identify grid shift file used by osgeo in Python?

texts:
I wrote a Python script reprojecting points from a source coordinate reference system (SCRS) to a target one (TCRS) both identified by an EPSG code.
I would like my python script to retrieve details about the grid shift file it uses when doing reprojections from NAD83 to NAD27.

Such information is given by CS2CS in a terminal window when the variable PROJ_DEBUG is set to ON.

Can my python script retrieve the same information?
Any ideas?

Thanks

My script is:

from osgeo import osr

#Input coordinates
xin = 436898.0
yin = 5577554.0

# Transformation setup from NAD83 UTM Z15 to NAD27 UTM Z15
srs = osr.SpatialReference()
srs.ImportFromEPSG(26915)
trs = osr.SpatialReference()
trs.ImportFromEPSG(26715)
rsTrans = osr.CoordinateTransformation(srs,trs)

# Doing transformation
(xout, yout, zout) = rsTrans.TransformPoint(xin,yin)


texts: Convert Vector field marker to line shapefile

texts:
I am using QGIS. I have a set of data with bearings from set points. I have used vector field markers to show them as lines from the point in the direction of the bearing but I want to have this as a line shapefile so I can manipulate the data and identify intersections with other files.

I am quite new to QGIS so it may be easy to sort this in a way I haven't noticed, but I couldn't see it and I didn't find a question on here where it has been answered.


texts: database-backed WMS service runs as transaction

texts:
I'm trying to tame a geoserver 2.3.2 installation (backed by Postgres 9.1 and Postgis 1.5) by determining which of the many, many layers are actually being used.  To this end, I've set up an audit table that logs SELECT activity on the tables I think might be deprecated.

However, when I make a WMS request to a layer that should be audited, I don't get any records in my audit table.  Of course, auditing SELECT in PostgreSQL isn't foolproof because INSERTs into the audit table are susceptible to ROLLBACKs.  Indeed, when logging all queries from geoserver to postgres, I see that geoserver is querying data in a transaction and then rolling back

bind S_1: BEGIN
execute S_1: BEGIN
parse <unnamed>: SELECT ...
bind <unnamed>: SELECT ...
execute <unnamed>: SELECT ...
parse S_2: ROLLBACK
bind S_2: ROLLBACK
execute S_2: ROLLBACK

There doesn't appear to be any reason to run a SELECT inside a transaction, and even less of a reason to ROLLBACK the transaction as there are no errors.  Is there a setting to turn off this transactional behavior when responding to a WMS request?

I like the database audit solution because it is easily searchable and reportable.  As a work-around I could certainly just turn on verbose logging at either the database or geoserver level and grep the logs for layer names.


texts: Creating a trigger containing St_intesects in Postgres

texts:
I am new to SQL. I need a help creating a simple trigger function that can contain records from a intersection.
I have points layer(vnp14imgtdl_nrt_global_24h) and a polygon boundary(states).
The trigger needs to fetch all new records that are inserted into jharkhandfire layer that overlap the polygon boundary into the new table.
Is it possible?

The trigger function is shown below:

    CREATE OR REPLACE FUNCTION public.firefunc()
    RETURNS trigger
    LANGUAGE 'plpgsql'
    VOLATILE
    COST 100
AS $BODY$   BEGIN
      select(new.gid, new.latitude, new.longitude, new.bright_ti4, new.scan, new.track, new.acq_date, new.acq_time, new.satellite, new.confidence, new.version, new.bright_ti5, new.frp, new.daynight, new.geom);
      into jharkhandfire
      from public.vnp14imgtdl_nrt_global_24h
      where st_intersects(new.geom, states.geom);
      RETURN NEW;
   END;
$BODY$;

But this function gives error: error near into


texts:
How can a WMTS layer be added to QGIS?
I'm trying to add World Imagery (Wayback 2014-02-20) to QGIS 3.16.
I'm using the following WMTS link...
https://wayback.maptiles.arcgis.com/arcgis/rest/services/World_Imagery/WMTS/1.0.0/default028mm/MapServer/tile/10/%7Blevel%7D/%7Brow%7D/%7Bcol%7D

Here's the error I'm receiving.


texts: Setting input parameter to be float for Python script tool?

texts:
I made a Python (2.7.8) script, where one of the parameters is height (float) and I wan to make tool in ArcGIS 10.3 now. However in data type list there is only double. If I input integer value into this parameter field (e.g. -480), the tool works fine. If i input e.g. - 480,5 i get error

  Traceback (most recent call last):   File "D:\HETMAN\GIS\Projekt Model
  Złoża\Sktypy.py", line 27, in  ValueError: invalid literal for
  float(): -480,1

  Failed to execute (SSKv1).

This is a fragment of code, where i try to convert value to float:

kontakt_dbl = sys.argv[3]
kontakt = float(kontakt_dbl)

How to convert double value into float and not to get error?


texts: binding nested arrays as geojson popups in leaflet

texts:
I am struggling with geojson popups in leaflet.js
I have nested arrays as properties for my features -
an example would be counts of fish caught at a location (not my actual project data!)

{"type":"FeatureCollection",
"features":[{"geometry":{"type":"Point",
"coordinates":[-122.679416,
45.516632]},
"type":"Feature",
"properties":{"fisher_id":"401",
"river":"Willamette",
"fishcaught":[{"type":"Salmon",
"count":5},
{"type":"Sturgen",
"count":1},
{"type":"Bass",
"count":2}]}},
{"geometry":{"type":"Point",
"coordinates":[-121.318584,
44.050165]},
"type":"Feature",
"properties":{"fisher_id":"401",
"river":"Deschutes",
"fishcaught":[{"type":"Trout",
"count":10},
{"type":"Whitefish",
"count":1}]}},
{"geometry":{"type":"Point",
"coordinates":[-124.102962,
43.983057]},
"type":"Feature",
"properties":{"fisher_id":"401",
"river":"Siuslaw",
"fishcaught":[{"type":"Steelhead",
"count":2},
{"type":"Salmon",
"count":6}]}},
{"geometry":{"type":"Point",
"coordinates":[-121.776289,
42.220382]},
"type":"Feature",
"properties":{"fisher_id":"401",
"river":"Klamath",
"fishcaught":[{"type":"Trout",
"count":3},
{"type":"Steelhead",
"count":2},
{"type":"Salmon",
"count":3}]}}]}

I am trying to put a popup on each feature (onEachFeature) that states river and a count of each type caught.

I'm binding using a function like this

function onEachFeature(feature, layer) {
        if (feature.properties.fishcaught) {
            var popupcontent;
            popupcontent = "River: "+feature.properties.river
            for (i=0;i<feature.properties.fishcaught.length;i=i+1){
                var fishstring = "
"+feature.properties.fishcaught[i].type+": "
                +feature.properties.fishcaught[i].count
                popupcontent = popupcontent+fishstring;
            }
            layer.bindPopup(popupcontent);
        }
    }

geojson = L.geoJson(inData, {
    onEachFeature: onEachFeature
}).addTo(map);

The problem is that the onEachFeature is caught in an endless loop. If I remove the for loop from the function, everything works fine.

Any thoughts on this?


texts: How can I sort table by clicking  one by one feature, now I am selecting the features but arc map sorting it automatically  according to OBJECTID.  I want, it should be sorted by clicking then I should give a sequence number to the features.


texts:
I downloaded the shape file https://dati.comune.milano.it/media/rdp/comunemilano/dati/MM_FERMATE.7z, it represents the Milan's subway stations. I used ogrinfo to output the shapefile's summary and guess the SRID:

$ ogrinfo -summary MM_FERMATE.shp MM_FERMATE

INFO: Open of `MM_FERMATE.shp'
  using driver `ESRI Shapefile' successful.

Layer name: MM_FERMATE
Geometry: Point
Feature Count: 100
Extent: (1506722.964600, 5027236.460200) - (1534186.769400, 5043728.155000)
Layer SRS WKT:
PROJCS["Monte_Mario_Italy_zone_1",
GEOGCS["GCS_Monte Mario",
    DATUM["Monte_Mario",
    SPHEROID["International_1924",6378388,297]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",9],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",1500000],
PARAMETER["false_northing",0],
UNIT["Meter",1]]
ID: Integer (10.0)
FERMATA: String (100.0)
N_ARCHI: Integer (10.0)

Then I loaded the data in PostgreSQL:

shp2pgsql -s 3003 MM_FERMATE.shp | psql -d gis_sample

Finally I tried to convert "the_geom" field in degrees but I always got the same results.

Output using st_astext function:

gis_sample=# select st_astext(the_geom) from mm_fermate where n_archi = 4;
    st_astext
----------------------------------
POINT(1515086.3214 5037713.221)
POINT(1516846.3194 5036826.4702)
POINT(1515864.6252 5036806.3372)
POINT(1514766.6816 5034612.3064)
POINT(1513860.5034 5034993.1312)
(5 rows)

Output using st_asewkt:

gis_sample=# select st_asewkt(the_geom) from mm_fermate where n_archi = 4;
        st_asewkt
--------------------------------------------
SRID=3003;POINT(1515086.3214 5037713.221)
SRID=3003;POINT(1516846.3194 5036826.4702)
SRID=3003;POINT(1515864.6252 5036806.3372)
SRID=3003;POINT(1514766.6816 5034612.3064)
SRID=3003;POINT(1513860.5034 5034993.1312)
(5 rows)

How can I get the lat/lon in degrees?

Below there's my PostGIS version:

gis_sample=# select PostGIS_Full_Version ();
                    postgis_full_version
-------------------------------------------------------------------------------------------------------
POSTGIS="1.5.3" GEOS="3.3.3-CAPI-1.7.4" PROJ="Rel. 4.7.1, 23 September 2009" LIBXML="2.8.0" USE_STATS
(1 row)