texts: Changing QGIS plugin's icon, why doesn't it change in the menu/toolbar?

texts:
I've been writing a QGIS plugin (using the "Plugin Builder" plugin), and decided I wanted to change the icon (icon.png).

After running make deploy to rebuild the plugin locally, the new icon appears in the Plugin Manager (both in the list of plugins, and the plugin description panel). So far, so good.

However, the old icon persists in the menu, and toolbar icons.

I've established that the icon.png is correct in the metadata.txt, and is getting copied into my ~/.qgis/python/plugins/myplugin directory.

This is especially puzzling - I'm not replacing the default plugin icon, but an icon I previously created to replace the default plugin icon. So this has worked before... :/

I've tried these...

restarting QGIS
uninstall, then reinstall plugin using plugin manager
destroy the plugin using make derase followed by make deploy
reload the plugin using the "Plugin reloader" plugin

Does QGIS cache icons somewhere? That's the only explanation I can think of. If it does, is there any way to flush the cache?

I'm using QGIS 2.10.1 Pisa on Ubuntu 14.04. I'm using "Plugin Builder" version 2.8.3

I suspect it won't be a problem for anyone who hasn't installed the plugin before, but it's a bit strange.


texts: mapserv calculate BBOX instead of requiring user to provide it

texts:
We have a requirement to email out a map showing a particular feature in our database.  When the query is made to generate the map all it has is a feature id.  All the queries I've seen to MapServer/WMS include a BBOX.  Is it possible to calculate the BBOX  based on the feature id?

The feature id will be stored in a PostGIS/PostgreSQL database.


texts: How to keep leading zeros in csv file when working with QGIS

texts:
I'm new to QGIS, and Geographic Information Systems in general, and also have no experience with informatics.

I a have shapefile which works fine, and then I have a file in CSV format, and I try to join these two files, using a column which should be equal between two files. The problem is, some of the numbers have leading zeros (they are some sort of postcode in my country), that disappear from the CSV file when I import it to QGIS...

Is there any simple way to keep leading zeros?

In the shapefile the leading zeros don't disappear.


texts: Google Earth Engine: Stack an ImageCollection

texts:
Does anyone know how can I stack all this collection?

I made my own collection from January Landsat images (L7 and L5) and I would like to stack all the images to check which are the areas with more data (without clouds).

var col = imageCollection
  .sort('system:time_start')
  .filterBounds(region)
  .map(function(image) {
    var date = ee.Date(image.get('system:time_start'));
    var savi = image.expression(
      '(1 + L) * float(nir - red)/ (nir + red + L)',
      {
        'nir': image.select('B4'),
        'red': image.select('B3'),
        'L': 0.5
      });
      var mask = image.select(['cfmask']).neq((1,2,3,4));
    return ee.Image(1).addBands(savi).updateMask(mask);
  });


texts: FME assigning sequential start and end for each line segment

texts:
To simplify, I have a line that is 500 meters long.  As part of its attributes it has a attribute called 'Start Chainage' of 0 and an attribute called 'end chainage' of 500.

In FME I have used snipper to cut the line into lengths of 20 meters.
What I need is to give each one of those segments their own sequential start and end changes.

So the first segment will have a start chainage of 0 and an end chainage of 20, the next segment will have a start chainage of 20 and end chainage of 40 and so on.....

The reason for this is so that I can identify at what point (chainage) along the road some other feature crosses or intersects it.


texts: Creating Metadata for new feature class using ArcPy with ArcGIS Pro

texts:
Using Jupyter Notebook implementation in ArcGIS Pro 2.8.0, I created a new feature layer and want to add metadata; the metadata object is read-only, though.
How do I address this problem?
All variables have been set accurately:

outFeatureClass references a valid feature class (i.e. myFeatureClass)
countyBoundary references a valid feature class (i.e. myCountyBoundary)
outFeatureClassFullPath had the complete path to the feature class (i.e., c:\temp\myGDB.gdb\myFeatureClass)

The call to CopyFeatures_management successfully creates the feature class referenced by outFeatureClass as demonstrated by the call
'''
tgt_item_md = md.Metadata(outFeatureClass)
'''
import arcpy
from arcpy import metadata as md

arcpy.CopyFeatures_management(countyBoundary, outFeatureClassFullPath)
#Add metadata
new_md = md.Metadata()
new_md.title = 'My Title'
new_md.tags = 'My Tag1, My Tag2'
new_md.summary = 'My summary'
new_md.description = 'My description'
new_md.credits = 'My credits'
tgt_item_md = md.Metadata(outFeatureClass)
if not tgt_item_md.isReadOnly:
    tgt_item_md.copy(new_md)
    tgt_item_md.save()
else:
    print("Failed to add Metadata to " + outFeatureClass)

The result is always: Failed to add Metadata to the feature class.


texts: Making SQL join query for NULL in search cursor in ArcPy?

texts:
I am trying to learn how to do some error catching/avoiding in ArcPy.  Here is a simplified version of what I would like to do.  I would like to take some points, and then check to see if there is a line in a line feature class with a matching ID to the point, and the same thing to match each point to a polygon.  If each point doesn't have a matching line and polygon, I want to skip that point and move on to the next.  Can someone tell me what I'm doing wrong?  Here is an example script showing how I am trying to do this:

# Import arcpy module
import sys, os, time, arcpy, smtplib
from arcpy import env
from arcpy.sa import *

env.overwriteOutput = True

env.workspace = r"D:\WorkingDirectory\Geoprocessing\Work.gdb"
PointFC = r"D:\WorkingDirectory\Geoprocessing\Work.gdb\points"
LineFC = r"D:\WorkingDirectory\Geoprocessing\Work.gdb\lines"
PolygonFC = r"D:\WorkingDirectory\Geoprocessing\Work.gdb\points"

Line_idField = arcpy.AddFieldDelimiters(LineFC, line_id)
Poly_idField = arcpy.AddFieldDelimiters(PolygonFC, poly_id)

try:
    rows = arcpy.SearchCursor(PointFC)
    for row in rows:
        Pt_id = row.pt_id
        PTpoly_id = row.poly_id
        PTline_id = row.line_id

        #This section if for skipping null values to help prevent crashes
        arcpy.MakeFeatureLayer_management (PolygonFC,  "Polylyr")
        arcpy.AddJoin_management(row, "poly_id","Polylyr", "poly_id","KEEP_COMMON")
        #This is where I get an error - it doesn't like the IS
        PolyIDnull = Poly_idField IS NULL
        arcpy.AddMessage('\nThis site does not have a polygon : ' + str(row.pt_id))
        if PolyIDnull:
            continue

        arcpy.RemoveJoin_manaement(row, "Polylyr")

        arcpy.MakeFeatureLayer_management (LineFC,  "Linelyr")
        arcpy.AddJoin_management(row, "line_id","Lineylyr", "line_id","KEEP_COMMON")
        LineIDnull = Line_idField IS NULL
        arcpy.AddMessage('\nThis site does not have a line : ' + str(row.pt_id))
        if LineIDnull:
            continue
        arcpy.RemoveJoin_manaement(row, "Linelyr")

except Exception, e:
import traceback
map(arcpy.AddError, traceback.format_exc().split("\n"))
arcpy.AddError(str(e))`


texts: Trimble 7X GNSS point Z returns negative elevations

texts:
I have a newly acquired Trimble 7X (with a centimeter kit).  I collected a point in an area with a clear view of the sky. The accuracy indicator showed 1" accuracy and post processing and a differential correction places the point spot on when displayed over imagery.  The settings on the Trimble are set to return Z values in heights above the ellipsoid in meters and output x and y in lat and long.

I have LiDAR data for this area referencing NAVD88. I would like to convert the Z values returned from the Trimble to NAVD88.  I tried this conversion using NOAA's vDATUM.

Here is the post process GNSS data in lat, long, & Z
(48.793770753, -122.612074293, -16.756m)

The LiDAR elevation at this locations is 4.673m NAVD88.

Below is a screen shot of the vDatum output.  The returned height is not even close to the LiDAR NAVD88 height in this location. My question is why is the GNSS Z value converted to NAVD88 so different from the vDatum value?  Am I doing something wrong in vDatum or am I expecting too much from the Trimble?


texts:
Sorry if my question feel stupid but I'm a new QGIS 2.4 user. My shape is composed of 7 polygons which are supposed to be adjacent (ie sharing sharing same nodes). Would like to join the polygons by their nodes and make sure there is no space in between two polygons. is there a way to do that?

Thanks.


texts:
I have a WMS layer and i want to get fid of polygon in it by click it. How can i do it?
I know that i can create WFS layer from this WMS and get feature attributes from layer or store. But its possible do this use just WMS? I dont want use WFS Layer couse it may contain ~30k features. But if i will use WFS Layer with visible=false selection control will be works?
I find SLDSelect control its can help me?
Maybe its possible to do like in gxp library?


texts: How do I reverse spplot colour key so the values are decreasing?

texts:
How do I reverse the colour key in a spplot so the that values are decreasing from the top to the bottom? Demo from https://edzer.github.io/sp/ below.

require(sp)
demo(meuse, ask = FALSE, echo = FALSE)
spplot(meuse, "zinc", do.log = TRUE,
    key.space=list(x = 0.1, y = 0.95, corner = c(0, 1)),
    scales=list(draw = TRUE))


texts: I added a label in the print composer of QGIS (3.6).  I rotated the label, but then the font is no longer aligned to a common line.  I think it is stair-stepping where the letters go, perhaps snapping them to discrete pixels or locations.  I messed around with the Guides and Grid settings, but cannot seem to fix the problem.  I have attached a screen shot, notice the label "survey grid".  It aligns well when I zoom in.


texts: OpenLayers Rule is not working as expected

texts:
I'm trying to add a rule to paint some roads of my GeoJSON file with a different style using OpenLayers.

The problem is that even after adding a property "featureType" to the feature in response to the event "featureadded", the rule doesn't seem to be working.

Any ideas why this happens?

This if a feature example in GeoJSON format after changing the "featureType":

{
    "type": "Feature",
    "properties": {
        "id": "42",
        "state": "texto",
        "featureType": "segmento_logradouro"
    },
    "geometry": {
        "type": "LineString",
        "coordinates": [
            [
                574951.47711036,
                9673700.3858644
            ],
            [
                574954.47530982,
                9673558.2712099
            ],
            [
                574814.15957496,
                9673551.0755312
            ]
        ]
    },
    "crs": {
        "type": "name",
        "properties": {
            "name": "EPSG:29193"
        }
    }
}

And this is the code I used to add the rules:

var ruleSegmentoLogradouroStyle = OpenLayers.Util.applyDefaults({
    fillColor : "#0000FF",
    fillOpacity : 0.9,
    strokeColor : "#FF0000",
    strokeOpacity : 0.9,
    strokeWidth : 5
}, OpenLayers.Feature.Vector.style["default"]);

var ruleSegmentoLogradouro = new OpenLayers.Rule({
    filter : new OpenLayers.Filter.Comparison({
        type : OpenLayers.Filter.Comparison.EQUAL_TO,
        property : "featureType",
        value : "segmento_logradouro",
    }),
    symbolizer : ruleSegmentoLogradouroStyle
});

var ruleDefault = new OpenLayers.Rule({
    symbolizer : OpenLayers.Util.applyDefaults({},
            OpenLayers.Feature.Vector.style["default"]),
    elseFilter : true
});

defaultStyle = OpenLayers.Util.applyDefaults({
    fillColor : "green",
    fillOpacity : 0.4
}, OpenLayers.Feature.Vector.style["default"]);

MapControl.editLayer.styleMap = new OpenLayers.StyleMap(defaultStyle);
MapControl.editLayer.styleMap.styles['default'].addRules([ ruleSegmentoLogradouro,
ruleDefault ]);


texts:
I have two steps here. First select layer by attributes. Second convert polygon to raster. It failed when it goes to second step and kept giving me error 999999: Error executing function. I pasted my codes below. Is there anyone who has the same problem. I am now using Arcgis 10.1 and has SP1 installed.

exec('''
import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("spatial")
import string
env.workspace = "G:\\countoverlap"

# Local variables:
FRAP_Fire111 = "G:\\countoverlap\\FRAP_Fire111.shp"

fldName = 'Year'
fcName = 'FRAP_Fire111.shp'
#set creates a unique value iterator from the value field
yearList = list(set([row.getValue(fldName) for row in arcpy.SearchCursor(fcName)]))
print(yearList)

 ''')

env.workspace = "G:\\countoverlap"
env.scratchWorkspace = "G:\\countoverlap\\overlap.gdb"
arcpy.MakeFeatureLayer_management(FRAP_Fire111, "lyr")

for year in yearList[1:len(yearList)]:

    print(str(year))

    # Process: Select Layer By Attribute
    arcpy.SelectLayerByAttribute_management("lyr","NEW_SELECTION", '"Year" = year' )
    print('Layer_Seclected')

    # Process: Polygon to Raster
    # Set local variables
    inFeatures = "lyr"
    valField = "Year"
    outRaster = "str(env.workspace)+'\\'+str(year)+'.tif'"
    assignmentType = "MAXIMUM_AREA"
    priorityField = "NONE"
    cellSize = 1080
    # Execute PolygonToRaster
    arcpy.PolygonToRaster_conversion(inFeatures, valField, outRaster, assignmentType, priorityField, cellSize)
    print('Raster_Created')


texts: ogr2ogr - Shapefile database SQL query; missing .cpg and .xml

texts:
I am using gdal and ogr2ogr to perform a SQL query on a database of shapefiles with the code below. The geoprocessing works and the appropriate shapefile is extracted from the query, but when I print boundary it shows 0 and the .cpg and .xml files are missing after the geoprocessing is completed. I also a get a warning: Value of field SHAPE_Area of feauture 0 not successfully written. Possible due to too larger number with respect to field width

import gdal
import subprocess

directory = r"E:\IrrigatedLands\FC_test\Water_Districts\Water_Districts.shp"
output_shp = r"E:\IrrigatedLands\FC_qgis\boundary.shp"

boundary = subprocess.call(["ogr2ogr", "-where", "NAME = 'Cache La Poudre River'", output_shp, directory])

print boundary


texts: Creating raster layers, weighting them, and combining them in QGIS

texts:
I'm looking to create raster layers from vector layers, rescale them, weight them then finally combine them for a final result.

I'm really not sure about the process to do this. On Arcmap it's really straight forward but on QGIS I'm having significant difficulty. This is the process I use on arcmap:
Vecter layer --> Eclidean Distance --> Rescale new layer using min-max score transformation on raster calculate --> combine all raster layers with Weighted Sum Tool.

I can't for the life of me figure out the same process in QGIS.


texts: Retail Huff Model & Network Analyst - Street/Travel Time Network error/issue/question

texts:
I am writing in regards to the Huff Model that Drew Flater posted on ESRI's website since I seem to be running into a problem whenever I try to run it with a network. I've tried emailing/contacting the author, but no luck.

The model is not generating outputs for all of my supply points, including my planned store location whenever I run the model with a network I've created. I've followed ESRI's tutorial on creating multi-modal/pedestrian networks.

I've got 172 supply points (cafes) + 1 (my potential site), but when I run the Huff Model toolbox the analysis only gives me outputs for 69 supply points? Why am I not getting outputs (probability & sales) for all 173 supply points? Screenshots are below.

The model runs fine whenever I do not include a network, but I need the Huff Model to account for travel time/multi-modal transportation networks...the example below is meant to account for solely walking since I'm modeling walking distances to and from local coffee shops.

Please note I have 141 Demand Points (polygons, neighborhoods) and 173 supply points in total:

NETWORK:

SUPPLY POINTS OUTPUTS AFTER RUNNING HUFF MODEL TOOL:

MAP & SUPPLY POINTS FULL LIST:


texts: Is there a way to have a single point feature with multiple records?

texts:
Is there a way to have a single point feature and have multiple attributes records  associated with it? Multipoint approach won't work as it gives me one attribute set (database row). I am trying to avoid having to stack points. Maybe my approach is flawed since I'm failing to find any relevant hits on a web search.

Use case: Access all details of personnel at company via a 'company building polygon centroid' for instance, by simply clicking on the point feature.


texts: QGIS 3.4.6 - symbology - Using Both size and colour to represent two pieces of information on single object

texts:
I would like to apply two different styles on the same vector layer (lines), from two different fields of the table.

To be more specific, I would like to make a graduate size of the lines from the X field, and a graduate color from the Y field.

There is a related question here : Using BOTH size and colour to represent 2 pieces of information on single point in QGIS

But I would like to know if the solution has been "updated" with QGIS 3.


texts:
I applied the Vector>Geoverarbeitungswerkzeuge>Auflösen (Dissolve) on a zoning plan (see picture zoning plan). I used the attribute village name as the dissolving criteria. The resulting geometry for the village (yellow colored) is sprinkled with small polygons or also larger polygons are not dissolved (see picture sprinkled geometry). In the attribute table the information is in one row.

Are there any approaches to obviate this? Changing any settings? I am using QGIS 1.7.2-Wroclaw.


texts:
I work with several MXD files. When opening those MXD files, i want to see the whole layout view, as this button function:
.

Is there a way, using arcpy, zooming to the whole page in the layout view, so in the next time i will open those MXD files, the default will be the layout view, zooming to the whole page?

I started to write this code, but i don't know how to continue with it:

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

env.workspace = r"G:\desktop\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
     print mxdname
     mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\\" + mxdname)
     df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
     mxd.activeView='PAGE_LAYOUT'
     print 'PAGE_LAYOUT'
     mxd.save()
del  mxd

This code don't zooming to the whole page, but change from data view to layout view, and that isn't my complete goal.


texts: Navigation/Routing in 3D View On Android

texts:
I need to build a navigation/routing application in 3D view like sygic, what are the available options if any?

Are there any existing open source or proprietary 3D mapping/libraries?

I have checked MapQuest, CloudMade and Openlayers2, but they don't seem to have this option.


texts: How to call a leaflet Javascript Function in dropdown menu bar?

texts:
I am having the following problem:

My target is a hybrid application with leaflet. I create a drop down menu with

html/CSS/JQuery.
Here is the code:

<image src="../img/Icon/info.png"/>

Interpretation

Kennblatt

            <image src="http://maps.ioer.de/cgi-bin/wms?MAP=S02RG_1000&&SERVICE=WMS&VERSION=1.3.0&SLD_VERSION=1.1.0&REQUEST=GetLegendGraphic&FORMAT=image/jpeg&LAYER=S02RG_2012_a"/>

<image src="../img/Icon/ebenen.png"/>

<image src="../img/Icon/kalender.png"/>

<image src="../img/Icon/lupe.png"/>

<image src="../img/Icon/helligkeit.png"/>

Now I like to call the Leaflet Geosearch Plug in from https://github.com/smeijer/L.GeoSearch at the point li id="search"

I know how to call the function within the JavaScript part, but not within the list element, which is outside the Leaflet function.

The code for the Leaflet Map is:

new L.Control.GeoSearch({
        position: 'topcenter',
        provider: new L.GeoSearch.Provider.OpenStreetMap(),
    }).addTo(map);

Which I like to implement in the drop down menu.


texts: Horizon line of a vector point from terrain raster - global mapper

texts: Is there a possibility to calculate the horizon line of a point with global mapper, not the viewshed, but just a line (or polygon) joining all the furthermost points from A and are at the same time seen from it.


texts: Adding multiple new values to MBTILES layer using `tile-join`

texts:
Using the tile-join function of tippecanoe, I have a MBTILES file with one point, and a CSV with one value to add to that point. This is what the CSV looks like:

I am able to edit the MBTILES layer by using the key to add the TIMESTAMP using the CSV file, although it removes the key value (OBJECT_ID) from the MBTILES file when it adds the TIMESTAMP in. I believe it's simply replacing one value with the other. Here's what I'm using in the command line:
tile-join -o updated_objects.mbtiles -x OBJECT_ID -c changes.csv objects.mbtiles

I've tried variations to this, such as removing the -x, although I haven't been able to keep the existing OBJECT_ID and the TIMESTAMP in the updated file.

Is there a way to add one or multiple values into an MBTILES file while still keeping the existing data?


texts: Calculating furthest distance from a pipeline line segment to a road line segment in QGIS

texts:
I am currently using the v.distance GRASS tool to calculate distances, but it naturally gives the nearest  (shortest) distances.  Oddly enough, I would like to find the furthest distance from a pipeline line segment to a road line segment.  Both shapefiles are Line (MultiLineString).
Is there any way for v.distance to be augmented to do such a calculation?


texts: Geoserver - change GeoWebCache Directory

texts:
can you help me please change caching directory of my geoserver instance? I've tried add to web.xml file in Geoserver\WEB-INF this :

  <!-- Setting GeoWebCache folder -->
  <context-param>
    <param-name>GEOWEBCACHE_CACHE_DIR</param-name>
    <param-value>\data\gwc</param-value>
  </context-param>
  <context-param>

But does not work ...
When I open Tile Layers page in Geoserver Web UI, there are still marked as N/A and when I look at GeoWebCache embedded page I saw my local storage is:

  C:\Users\Denis\AppData\Local\Temp\geowebcache

but in this folder Geoserver does not has permission to read or write, bacause I need my cache folder in geoserver installation folder.

I am using Tomcat 8.5 and Geoserver 2.10 on windows 10. Any ideas?


texts: Counting number of pixels/cells with particular value in raster using ArcPy?

texts:
I have a raster image containing values ranging from 1 to 10 in Integer format.

I want to count the number of pixel with value = 1.

Is there any method in Arcpy which returns the pixel count if given raster and value as input?


texts:
I want to divide a unique entity of SHP to multiple entities
The picture bellow shows a selected entity with area = 993ha
My objective is to divide this entity to multiple entities where each area is equal to 2ha

I am using QGIS.


texts:
[SOLVED] Good day!  I'm looking at either gdal_contour or gdal_polygonize.py to create a vector layer file from a raster file.  I can create the vector file using either method, but it's upside down when I plot it using QGIS.  Another list indicated this has to do with the world file.

I use gdal_translate and gdalwarp when transforming the orginal raster in lambert conformal to spherical mercator (for use in OpenLayers) and understand that the world file could be created using one of those utilities (probably when using gdalwarp).  However, I'm really unsure after that about the next step using either of the raster to vector utilities where the world file comes into play.  Do I have to also use ogr2ogr as well?  I'd appreciate a detailed list of steps as the GDAL documentation is a bit fuzzy here.

Just to note my current steps are:

Use gdal_translate to first geo reference the original raster in lambert conformal
Use gdalwarp to reproject the resulting raster image to spherical mercator
Use gdal_translate to create a PNG at a reduced size.  Becomes an Image Layer in OpenLayers.
Use either raster to vector utility on the raster produced in step 3.

I see where -co for the world file can be included in gdal_translate or gdalwarp.  And like I said, I can create the vector file once the original raster has been warped to spherical mercator.  Really just unclear here how the world file comes into play for the resulting vector file in step 4 so the image will be rightside up. Thanks for any help!

[SOLVED] As Matthew asked and then I commented on, I apparently was confused about the world file creation and in which step to do that in.

But, addding -co WORLDFILE=YES in step 3 above created a .wld file for the resulting PNG image.  Then used gdal_contour creating a shapefile and also gml, brought both into QGIS to view, and they're now rightside up!  Thanks all who viewed and Matthew who commented!


texts: How to create a dynamic distance calculator for selected line feature classes?

texts:
My situation:

I am relatively new to GIS so forgive me if this is a stupid question.  I am working on a project at a park that I am interning at where the goal is to create a single feature with roads and trails together with some way to show the total distance of selected roads and trails.

Currently I have merged the roads and trails features, so I have a single feature with both which contains separate fields with their respective distances.  In its current form I could simply add distances manually from selected trail and road intersections, but I was interested in automating a process to show this.  I was thinking some kind of model (made in ArcGIS model builder) that when run would return an output table of the summed distance of currently selected features.

Question:

1) How can I tell ArcGIS model builder to perform calculations on just selected features?

2) Is there a better way to accomplish the above mentioned task?

Thank you for any input.


texts: How can I zoom to featureset geometry item

texts:
I send GetFeatureInfo request to server with ajax and get the featureset data. I want to zoom a selected feature.

Object {type: "Feature", id: "mydata.39025", geometry: Object, geometry_name: "geom", properties: Object…}

Selected feature can be a point, multilinestring, polygon. So I can not cast it as an OpenLayers.Geometry.Point, or OpenLayers.Geometry.Polyline.

How can I zooom it?


texts: How is a "Map Package" created programmatically?

texts:
Given a loaded map document (mxd) how do I find the geodatabase objects related to the contained layers? I want to copy them with the map document so the map can be moved onto each user's local machine from a network share. Is there an API to do this? Can someone point to source code that does this?

To clarify (after all the great feedback): How is a "Map Package" created programmatically using the ArcEngine API? Is this even possible?


texts: qgis + spatialite: convert a list of points into lines and distances

texts:
I have a list of points stored in a spatialite database (track points from a gps log). I would like to transform these into a set of line segments and calculate the length of each segment. The original data is in CRS 4326 (WGS 84).
I think I have something that works, but I'd like to run it by folks here to make sure I'm not missing anything.
I start by importing the GPX file into QGIS by adding a vector layer. I then import this layer into a Spatiallite database. This results in a table that includes columns track_fid, track_seg_id, and track_seg_point_id for identifying each point (and a 3D POINT geometry column). During the import I transform the CRS from 4326 (WGS 84) to 26986 (NAD83 / Massachusetts Mainland), primarily because that way I get meters from st_distance.
I create a new table (track_lines) like this:
create table track_lines (
  id integer not null  primary key autoincrement,
  track_fid int,
  track_seg_id int,
  track_seg_point_id int,
  length int);
select addgeometrycolumn('track_lines', 'geom', 26986, 'linestring', 'XYZ');

Finally, I populate the track_lines table using the following SQL expression. In order to locate the endpoints of each line segment, I'm trying to "shift" the track points by joining on leftside.track_seg_point_id = rightside.track_seg_point_id-1.
insert into track_lines (
    track_fid, track_seg_id, track_seg_point_id, geom, length
    )
select
    a.track_fid as track_fid,
    a.track_seg_id as track_seg_id,
    a.track_seg_point_id as track_seg_point_id,
    makeline(a.geom, b.geom),
    st_distance(a.geom, b.geom)
from track_points as a join track_points as b
on (
    a.track_fid = b.track_fid
    and a.track_seg_id = b.track_seg_id
    and a.track_seg_point_id = b.track_seg_point_id-1
    );

This seems to work: the lines that result from the above query appear to match the track points. I'm a little suspicious of the join operation, because invariably when I show something like this to someone with more SQL experience they immediately reduce it in size by 70% while tripling the performance (fortunately in this case the number of points is low enough that performance isn't actually an issue).
Does this all seem sane?


texts: How can we extract the stream order branch?

texts:
Actually I am new to GIS, so I don't know much about technical terminology.

I have extracted a river network but now have a problem in distinguishing the  branch order of stream (from order 1 to 2, 2 to 3; side branches, not the main ones [say 1 to 1, 2 to 2]) according to the Strahler scheme. It is displaying only one column of stream order. Can anybody tell me how can I extract a separate column for side branch streams as I mentioned above?

Here is a Strahler ordering scheme. In the strahler order method, all links with no tributaries are assigned an order of one and are referred to as first order. When two first-order links intersect, the downslope link is assigned an order of two. When two second-order links intersect, the downslope link is assigned an order of three, and so on. When two links of a different order (say 1 to 2) intersect, the downslope link is assigned the order of the maximum of these two. My query is, can we make a separate field for the same order branch (1 to 1, 2 to 2) and different order branch (1 to 2, 2 to 3)?

I have added two tables. In first table we can see that all stream orders (Grid_code) are given in one column. I want to know how many streams of 1st order are joining with other 1st order, 2nd order, 3rd order, etc. Similarly, how many streams of 2nd order are joining with 2nd order, 3rd order, etc.

For example, in this table there are 12 '1st' order streams. I have to know that out of 12, how many are joining with another 1st order (1 to 1)? How many are joining with a 2nd order (1 to 2)? How many are joining with a 3rd order (1 to 3)? And so on.


texts:
In my company I am working with ArcGIS for Dektop on my PC with data stored at the company server (mapped path). I have trouble working with File-GDBs:

I can not import data to the GDB (only an empty FeatureClass is imported; Error:

  Error exporting data. Invalid function arguments Invalid function arguments Cannot acquire a lock. [GDB_Items] Cannot acquire a lock. The table was not found. [FeatureClassName] The table was not found. [FeatureClassName] The table was not found. [FeatureClassName] The table was not found. [FeatureClassName]

I can not add fields (rare) nor calculate fields (often); Error:

  Could not load data from the source. If you can correct the problem, press the refresh button to reload data. Possible problems can include bad network connection, invalid field, etc. Cannot acquire a lock. Cannot acquire a lock. Das Gleichheitszeichen zwischen Attribut und Attributwert fehlt.

The last sentence is German and means: The equals sign between attribute and attribute value is missing. However, I do everything right. If The GDB is on my PC (C:) everything works fine. I read somewhere, that this might be through a mix of mapped pathnames and UNC pathnames within a single mxd. However, I only have one kind of pathname (https://geonet.esri.com/thread/89084). I acctually have only one GDB in my MXD. Any suggestions?

I am currently using ArcGIS 10.4.1, but I had the problem before aswell (using 10.3).

Edit:
Currently, I am the only user of the GDB and MXD. Once I approve the GDB, several people need reeding access to it. The errors I mentioned, occur most of the time after that procedure:

I start Windows
I start my MXD, which uses one GDB on the server
I add a field (often but not always that works)
I calculate the field
I get the error message and cannot see the table anymore (it's emty now); the calculation worked, but I need to close and reload the table to see that.

I use Feature Datasets, but the errors occur in GDBs w/o Feature Datasets aswell. I am aware that there is a significant performance benefit to having FGDB data on the local filesystem. However, as other people need to have access to that data I don't see any other option to the server-option, which I am using for several years now. It always worked fine until now.


texts:
I am having trouble creating a buffer around a polygon of a county in QGIS. When I run the "fixed distance buffer" geoprocessing tool, no features are returned in the 'buffer' shapefile that is created. Has anyone had trouble with this as well? Or does anyone have suggestions of what I may be doing wrong?

Here is the polygon I am attempting to buffer:


texts: Openlayers 3 - get Geometry of a VectorTile to highlight features

texts:
I am trying to highlight features when mouse interaction is triggered over a VectorTile layer.

I have followed some examples but there is no one that makes what I would like to do. The most likely is this VectorTile properties example but it only shows how to get properties. At this point, I tried different things.

First, I have seen the OL3 hover vector layer example and I tried to apply to a VectorLayer without any success.

Then, I tried to get the geometry of the feature adding it to an empty VectorLayer which every time will contain only those features under mouse position. Then I could apply a custom style to highlight these features. I added the following function:

    ...

    map.on('pointermove', function(e) {

      for (var i = 0; i < highlightedFeatures.length; i++) {
        highlightedFeatures[i].setStyle(null);
      }

      highlightedFeatures = [];

      map.forEachFeatureAtPixel(e.pixel, function(feature, layer) {
          if (feature.setStyle) {

            // Here I change the syle for those features that belongs
            // to other type of layers (no VectorTiles)

            feature.setStyle(selectStyle);
            highlightedFeatures.push(feature);

          } else {
            console.log(feature.getProperties());
            console.log(feature.getType());
            // I get the polygons but they are very strange!
            console.log(feature.getGeometry());
          }
      });
    });

    ...

I realised about console.log(feature.getGeometry()); method is showing an object that seems having a strange coordinates arrays:

    hk {
    a : Array[4]
        0 : 16741
        1 : 28882
        2 : 28553
        3 : 4217
        length : 4
        __proto__ :  Array[0]
    b : Array[64]
        0 : 20751
        1 : 28942
        2 : 22853
        3 : 30064
        4 : 20135
        ... // a lot of "coordinates"
        length : 64
        __proto__ : Array[0]
    c : Array[1]
        closure_uid_123785828 : 23488
    f : Object
        codi : "CF7692N9826901"
        f_id : 5849
        layer : "coverage"
        __proto__ : Object
    g : "Polygon"
        __proto__ : Object
    }

Does anyone know why the coordinates are wrong? If I add the feature I get an error. Maybe someone knows how to get Geometries of VectorTile layers in a different way.

I am using OL3 v3.17.1 and I am using srs 3857.

I guess is not important but the highlighting style referenced variables before are:

    var highlightedFeatures = [],
    selectStyle = new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: '#0d47a1',
            width: 12
        }),
        image: new ol.style.Circle({
            radius: 5.5,
            stroke: new ol.style.Stroke({
                color: '#0d47a1',
                width: 12
            })
        })
    });


texts:
I would like to use some styles downloaded from the internet in .stylx format within QGIS.
Is there an easy way to do this?
This is the link to the style I would like to use to create maps:

watercolor maps


texts: Choosing video card in laptop to use with ArcScene

texts: I need a new laptop within the next couple weeks to use at grad school for GIS. I am a geology & geography major and use numerous aerial and satellite photos at once and often use ArcScene. I know that I should probably get an HP zbook with a workstation quadro series video card, but that is a tad out of my budget. Would a gaming laptop (such as alienware or ASUS Rogue) with an NVIDIA geforce card work almost as well?


texts: Applying decimal place change in ArcGIS Pro?

texts:
I am trying to visualize the decimals of the data in column TPvsTIM. This column's data is the result of the following:

TPvsTIM = Temps_TP / ToBreak (both being long integers with numeric values)

Here is a picture of my table


texts:
I want to generate a Cloud Optimized Geotiff (COG) using Python (with Anaconda environment). I wanted to use GDAL 3.2.0 since it has a COG Built in Driver since 3.1. Initially I cannot update GDAL to 3.2.0, it just keep installing 2.4.4 I am not sure why. So I downloaded a 3.2.0 package here: https://anaconda.org/conda-forge/gdal/files.
I installed it to my environment:
conda install gdal-3.2.0*.tar.bz2

Now it is now versioned 3.2.0 according to my
conda list
...
gdal            3.2.0            py37h6bca1a1_0    <unknown>
...

However when I try to use it, COG is still not recognized.
gdal_translate input.tif output.tif -of COG -co COMPRESS=LZW

It gives the error:
Output driver `COG' not recognised.

This is based on: https://gdal.org/drivers/raster/cog.html
My system is running on Debian 10.


texts: Feature Class To Geodatabase => ERROR 000210: Cannot create output

texts:
error messages:

Executing: FeatureClassToGeodatabase E:\jags20120321.mdb\ch_cs "Database Connections\Connection to localhost.sde" "Database Connections\Connection to localhost.sde"
Start Time: Sat Jun 16 13:32:28 2012
Running script FeatureClassToGeodatabase...
Failed to convert E:\jags20120321.mdb\ch_cs. ERROR 000210: Cannot create output Database Connections\Connection to localhost.sde\jags.DBO.ch_cs
Failed to execute (CopyFeatures).

Completed script FeatureClassToGeodatabase...
Executed (FeatureClassToGeodatabase) successfully.
End Time: Sat Jun 16 13:32:37 2012 (Elapsed Time: 9.00 seconds)

Those are the error messages.  I don't know what's the reason.  I use arcgis 9.3, arcgis sde 9.3, and sql server 2008.


texts: Relationship between DXF's entities and objects

texts:
I'm trying to parse a raw .dxf file. Inside the Entities section I can see the polylines I'm looking for (they contain bus routes) and I inside the Objects section I can see each route's data (bus front sign, for instance). My problem is: I can't see how the relationship of entities and objects is defined, so I'm unable to link a route to its data.

Can anybody here explain how entities and objects are related?


texts:
Currently working on a project that among other things requires displaying NOAA ENC for a web based intranet portal. We downloaded S-57 ENC data from NOAA site and successfully imported them into Postgres DB using "ogr2ogr". We now need to serve them up with the right colors and symbols that must follow S-52 standard.

I have looked at OpenCPN and ArcGIS-ArcMap and other desktop tools which have beautiful renditions of S-57 data with colors, daylight/ night and symbols. I would like to do the same for a web app that is served by GeoServer. The styling is hopefully accomplished using Styled Layer Descriptor (SLD) files.

I've looked at Atlas Styler, OpenGeo Styler etc. Nice tools but i'd hate to hand-craft them as there are a ton of layers and features, not to mention validating all of them to get it right. It would be nice if there is a way to download or import SLD files from good sources that correspond to S-57 layers. What is the least painful way to render ENC data on Web?


texts: Earth engine error: Layer 1: Layer error: Image.visualize: Parameter 'image' is required

texts:
For land use classification of my shapefile, When I use Landsat 4 collection, got an error:

  Layer 1: Layer error: Image.visualize: Parameter 'image' is required.

but in the case of Landsat 8, there is no such problem.

Codes here:

var image = ee.Image(ee.ImageCollection('LANDSAT/LT04/C01/T1_TOA')
    .filterBounds(table)
    .filterDate('1984-04-01', '1984-12-31')
    .sort('CLOUD_COVER')
    .first());

var visParams = {
      bands: ['B3', 'B2', 'B1'],
      min: 0,
      max: 3000,
      gamma: 1.4,
};
Map.addLayer(image, visParams);

How do I solve this?


texts:
I am trying to get the length of a line between two selected points. The points reside in table subs and are utilizing srid 4236.

I can verify the selection exists...

SELECT id, name, geom FROM subs where  name = 'DALLAM COUNTY'or name = 'DALHART';

Which returns the expected results...

However when I try to construct line with the same selection I get empty results?

SELECT ST_Length(ST_Transform(ST_MakeLine(geom), 4236)) FROM
    subs where  name = 'DALLAM COUNTY'or name = 'DALHART';

The sql returns empty results.

What am I missing?


texts: How to make osm2po 5.1 write an SQL file?

texts:
With osm2po 5.1, I can get the osm2po router running

D:\osm2po-5.1.0>java -Xmx1024m -jar osm2po-core-5.1.0-signed.jar
prefix=vie ..\vienna_austria.osm.pbf

but the resuting vie folder only contains the osm2po graph .gph and no .sql file for import into PostGIS. I don't get any warnings and couldn't find any new instructions in the documentation related to writing .sql files.

In 5.0, the .sql file is created as expected.


texts: The value is not a Feature Layer. Failed to execute (AddRuleToTopology)?

texts:
I am trying to create topology with ArcGIS python .
when I do the AddRuleToTopology its show the error as The value is not a Feature Layer. Failed to execute (AddRuleToTopology).
But I am tried Manually with AddRuleToTopology_management tool its working
but when copied snippet to python getting same error.
my code:

test10=arcpy.CreateTopology_management(r'D:\Velugoti\test.gdb\test','test10',0.01)
arcpy.AddFeatureClassToTopology_management(test10,r'D:\Velugoti\test.gdb\test\test',1,1)
arcpy.AddFeatureClassToTopology_management(test10,r'D:\Velugoti\test.gdb\test\test2',1,1)
arcpy.AddRuleToTopology_management(test10,"Must Not Self-Intersect (Line)",r'D:\Velugoti\test.gdb\test\test')

Its shows
 ERROR 000840: The value is not a Feature Layer. Failed to execute (AddRuleToTopology).


texts: Land Cover vector data for Saskatchewan / Manitoba Canada

texts:
Looking for Land Cover (aka Land Use) vector data for Saskatchewan / Manitoba Canada.
Using Saskatchewan's site I can see the data but not download it as seen here:
Another source is Data Basin but they also seem to only have it in viewer format.


texts: Put output of SAGA algorithm of Processing model in R script

texts:
I'm trying to combine the Variogram Cloud model of SAGA with an R script inside a Processing model.

First, the R script is:

##Vector processing=group
##showplots
##Layer=table
##X=Field Layer
##Y=Field Layer
plot(Layer[[X]], Layer[[Y]], xlab=X, ylab=Y)

Second, the Processing model:

where:

Shape is a generic shapefile
attribute is the attribute for the variogram calculation
distance is just a number for the default distance of the SAGA algorithm

The output of the SAGA algorithm is a simple csv table, here a screenshot:

What I want to do now is to put inside the R script (in the model) the Distance and Semivariance columns as x and y axis.

The only way I found is to add manually Distance and Semivariance in the model:

So the question is: how can I add default x and y axes in the Rscript window?

I think link a dropdown menu for both X and Y where I can choose the column of and algorithm output.


texts: Increase the number of bookmarks in a QGIS project

texts:
In QGIS 3.0 and 2.18 there is apparently a limited number of bookmarks, that can be stored. It is a limit of 250 bookmarks per project.
I find bookmarks great for an in-project team coomunication and I need to change that limit somehow, as we reached it.
Made a lot of how-to research and didn't found anything similar. The only path I found possible is to change "something" in the Advanced Settings Editor.

Does someone know how to change the limit of possible bookmarks in a QGIS project?


texts: Creating custom CRS based on local project parameters in QGIS?

texts:
I'm having trouble creating custom proj4 CRS files to work with existing .shp files on local project coordinates. (Local in that the coordinates are relative only to themselves and not the 'real world')

Right now we are wanting to export property lines from QGIS into AutoCAD.  The property lines were provided to us with a .prj file using NAD83 State Plane Coordinates.  In the past, we have utilized ArcGIS (or notepad) to create a projected local coordinate system based on our project parameters and assign its .prj to the data frame of our .shp file, and tell ArcGIS to use that .prj for the output. Done.

The idea is that we would like to be able to create custom CRS in proj4 and be able to utilize that to export 'data' from QGIS and have it relate to our local project coordinates. Vice versa, we would like be able to export CAD 'data' (obviously on the local coordinate system) and assign the custom CRS to have it relate to the rest of the world (or at least my state)

For example, I generally set up our project coordinate systems based around State Plane Coordinates on ground distances(for quantity purposes).  I generally "lop off" the first couple digits off from the SPCS in an effort to make the coordinates a bit easier to work with(knowing that most people we work with find it easier with 4 or 5 digits rather than 6 or 7).  The problem herein lays that our point of "origin"(aka the point where Northing = 0 and Easting = 0) does not fall down in Northeastern New Mexico where my state has its SPCS origin.  I need to be able to tell QGIS a false northing/easting, a scale factor, a rotation, and a lat/long origin to accurately project our project work onto the 'real world'


texts:
I have worked with .shp, geodatabases and WMS formats in QGIS but I am new to working with WFS format. I have QGIS 2.18 installed on my computer.

I am trying to open and view U.S. NOAA NHC WFSserver,  http://idpgis.ncep.noaa.gov/arcgis/services/NWS_Forecasts_Guidance_Warnings/NHC_Atl_trop_cyclones/MapServer/WFSServer, in QGIS.   Actually I get the connection to work,  I can add the layers I want to Layers Panel, I can view the attribute table and if I do a 'zoom to layer' the screen moves to approximate location of where I expect to see the data layer.  I just cannot see the actual features on the map.

I am not sure what the issue is? Has anyone worked with this data series before.


texts: Using multiple selection queries to calculate fields in ArcMap ModelBuilder

texts:
I am trying to calculate two fields based on a combination of select by attribute and select by location in model builder for road center lines.  In my current model I have it set up so that I initially select my boundaries using select by attribute and then perform a sub selection to only select the road center lines that are within the previously selected boundary.  From there I have it calculate two fields based on that selection.  However, I am noticing that once it runs the first query and calculation, the second query and calculation still contains the selected fields from the first query and thus overwriting the initial calculations.  The first image is the entire model while the second and third image refers the segments that involve my selection query and calculate fields.  Every time I run the model the the first query and selection results in over 100,000 fields selected which is correct.  The following query results in over 100,000 fields, but that is incorrect as the boundary that is used in the first part of the query is smaller in comparison the boundary used in the initial selection.  The results should be between 10,000-9,000 fields selected.  I am also running into cases that the sub selection from select by location does not produce any results so it will run without selecting any entries.

Any advice on how to fix this?


texts:
I've installed GeoServer and uploaded a shapefile in it and i want to export this added layer to shapefile using GeoTools library.

I've writtren some code which generates a shapefile(shp,dbf,shx,prj) that contains features without geomertic shape. in the other word it just exports attributes.

pls help me with this problem. where should i revise my code.

here is my code

try {
        String getCapabilities = "http://localhost:8060/geoserver/ows?service=wfs&version=1.0.0&request=GetCapabilities";
        Map connectionParameters = new HashMap();
        connectionParameters.put("WFSDataStoreFactory:GET_CAPABILITIES_URL", getCapabilities);
        DataStore data = DataStoreFinder.getDataStore(connectionParameters);

        String[] typeNames = data.getTypeNames();

        for (int i = 0; i < typeNames.length; i++) {

            if (!typeNames[i].contains("navahi"))
                continue;
            System.out.println(typeNames[i]);
            SimpleFeatureSource wfsFeatureSource = data.getFeatureSource(typeNames[i]);

            if (wfsFeatureSource instanceof SimpleFeatureStore) {
                SimpleFeatureType TYPE = wfsFeatureSource.getSchema();

                File newFile = new File(System.getProperty("user.dir") + typeNames[i].replace(":", "_") + ".shp");

                System.err.println(newFile.getAbsolutePath());

                ShapefileDataStoreFactory dataStoreFactory = new ShapefileDataStoreFactory();

                Map<String, Serializable> params = new HashMap<String, Serializable>();
                params.put("url", newFile.toURI().toURL());
                params.put("create spatial index", Boolean.TRUE);

                ShapefileDataStore shpDataStore = (ShapefileDataStore) dataStoreFactory.createNewDataStore(params);
                shpDataStore.createSchema(TYPE);
                String shpTypeName = shpDataStore.getTypeNames()[0];
                SimpleFeatureSource shpFeatureSource = shpDataStore.getFeatureSource(shpTypeName);

                SimpleFeatureStore shpFeatureStore = (SimpleFeatureStore) shpFeatureSource;

                SimpleFeatureCollection wfsFeatureCollection = wfsFeatureSource.getFeatures();

                Transaction transaction = new DefaultTransaction("create");
                shpFeatureStore.setTransaction(transaction);

                try {
                    shpFeatureStore.addFeatures(wfsFeatureCollection);
                    transaction.commit();
                } catch (Exception problem) {
                    problem.printStackTrace();
                    transaction.rollback();

                } finally {
                    transaction.close();
                }

            }

        }
    } catch (Exception e) {
        // TODO: handle exception }
    }


texts:
I have a polygon layer with a number as label. Now I want that my labeling is always in the left upper corner of the single polygon, like in the screenshot.The red rectangle is the position, where the label should be placed

If I zoom in the should be visible and placed like in the second screenshot.

If such a labeling is possible, how can I create that?
I tried the labeling tab of the gui, but it doesn´t work pretty well
With my settings the labels of other polygons where placed in their visible neighbours


texts:
I would like to delete the space in this number like this 'F00 089' with the field-calculator. Is theres a function in the field-calculator?

!
see my screenshot


texts: Saving file to .shp with have many attribute

texts:
Is there a way to save a file to .shp with have many attribute (I use QuickOSM to create a building Layer and I tried to save but it get error).
Is there a way to change .KML to .shp?
I get error message:

Could not make temporary scratch layer permanent. Error: Creation of
field addr:street_1 failed (OGR error: Can't create field addr_str_4
in Shape DBF file, reason unknown.)


texts: How can I combine these two commands? (Get a clipped raster ASC-II file from PostGIS)?

texts:
My objective is to clip a polygon from PostGIS raster table, and then save this raster (draw polygon shape) as an ASC-II in a folder. What I did now are 2 steps.

1st Step: Insert the clipped raster into another raster table using Pgresql in PHP:

    $preparedQuery = "INSERT INTO rastertable(rast) (SELECT ST_Clip(rast,1,geometry) FROM dtm, polygons WHERE polygons.polygonid = ?);";
    $sth = $dbh2->prepare($preparedQuery);
    $sth->execute(array(
    $geoid
    ));

2nd Step: Retrieve the raster data from this table using PostGIS dumper:

gdal_translate -of AAIGrid "PG:host=localhost port= dbname='' user='' password='' schema='' table='rastertable' where='rid=1'" postgis.asc ;

How can I combine these two processes? I mean, can I get the raster file directly without inserting to and retrieving from the raster table?


texts: Orfeo Toolbox (OTB) KMeans classification issues QGIS 2.4

texts:
I am having trouble processing an image through the Unsupervised KMeans image classification in the OTB Toolbox of QGIS 2.4. I have a .tif image in the Input Image box and keep all the defaults.

When I run it, this error immediately shows up:

  "Algorithm Unsupervised KMeans image classification starting...
  unsupported operand type(s) for +: 'NoneType' and 'str' See log for more details"

I have tried increasing the available RAM (Mb) option as well as saving the Output Image in a file rather than a temp. I am relatively new to QGIS so I apologize if I have overlooked anything simple. Can anyone shed light on this problem for me?


texts: GDAL calculation stop working python.exe any time

texts:
i try to use python GDAL raster calculation to take some raster result where i need.

i windows user and i work in OSGeo4W Shell.

first i try to work using batch file with success.

code work with batch file :

import os
os.system("C:/xxx/xxx/path/test.bat ")

batch file :

gdal_calc.py -A C:/xxx/xxx/path/a_layer.tif -B C:/xxx/xxx/path/a_layer.tif --outfile=C:/xxx/xxx/path/a_layer.tif --calc="(2*A)+(2*B)"

all work without error.

but if i try to wotk using python then any time before complete my script then python.exe crashed.any time

try using python script :

call('python C:/OSGeo4W64/bin/gdal_calc.py -A {0} -B {1} -C {2} --outfile={3} --calc="(2*A)+(2*B)"'.format(a_layer,b_layer,out_layer))

i want to use python script because i have many options for working usnig python and gdal.

but can someone to tell me why batch file work and not python script ?

edit

if print 'PATHEXT' environment variable in OSGeo4W Shell i dont have '.PY' extesion.but if is this hte problem can i define 'PATHEXT' inside some python script ?


texts: Defining atlas feature margin using PyQGIS3

texts:
I want to create a simple atlas composer map using PyQGIS3. My question is about atlas feature scale. I want to use this settings from image.

I know this select option but I don't know how to add some number percent for margin in first option ( layout.referenceMap().setAtlasScalingMode(QgsLayoutItemMap.Auto) )

layout.referenceMap().setAtlasScalingMode(QgsLayoutItemMap.Auto)
layout.referenceMap().setAtlasScalingMode(QgsLayoutItemMap.Predefined)
layout.referenceMap().setAtlasScalingMode(QgsLayoutItemMap.Fixed)

Any idea how to do this in PyQGIS3?


texts: Issue with using MultiLineString in Leaflet Animated Marker Plugin

texts:
This response on image below is coming from Geoserver as MultiLineString and I restructured it to LineString like on image because Animation plugin doesn`t support MultiPolyline json structure and because response has coordinates as lng, lat and I need lat, lng as Leaflet defines. This structure on image works nice with animation plugin.

Since Leaflet Animated Marker Plugin I use doesn't support MultiLineStrings I had to structure it as LineString. It accepts this LineString structure but coordinates are mixed after I apply coordsToLatLngs static method.
That makes my marker not appearing on first position and jumping back and forth. Until this line of code coordinates are in right order, as defined in SQL View on Geoserver and I could check it and prove in DB and in QGIS

This is view on Geoserver:

SELECT pgr.seq AS seq,
   pgr.edge AS edge,
   e.id AS id,
   SUM(e.distance) AS distance,
   st_linemerge(st_union(e.geom)) AS geom
  FROM pgr_dijkstra('SELECT id, source, target,onewayNum as cost, r_cost as
   reverse_cost FROM parking.n_streets_noded',%source%,%target%, true) AS
   pgr, parking.n_streets_noded  AS e
 WHERE pgr.edge = e.id
 GROUP BY  pgr.seq, e.geom,e.name,e.type, pgr.seq, pgr.edge,e.id
 ORDER by pgr.seq

This is function with ajax call:

function getFinalRoute(){
var urlRoute = `${geoserver}/wfs?service=WFS&version=1.0.0&request=GetFeature&
typeName=xxx:shortestpath&viewparams=source:${source};target:${targetN || targetE}&outputformat=application/json
&srsName=EPSG:4326`;
var routeLayer = L.geoJSON(null);
var flippedCoords;
$.ajax({
    url: urlRoute,
    async: false,
    success: function(response){
        console.log(response)
                var routeArr = response.features;
                var coordsArr = Object.keys(routeArr).map(key => {
                    return routeArr[key]
                })

                var xxy = coordsArr.map(function(feature){
                var obj = feature.geometry.coordinates[0];
                    return Object.keys(obj).map(function(key){
                        return obj[key];
                    })
                })
   //code above gives me coordinates in perfect order but after applying this method below their order mixes up
                var flipCoor = L.GeoJSON.coordsToLatLngs(xxy,1);

                var flippedCoor = flipCoor.sort().map(function(obj){
                    return Object.values(obj).sort().map(values => {
                        return Object.values(values)
                    })
                })

              var myArray = [].sort().concat(...flippedCoor);

                console.log(myArray);

                var multiLineString = L.polyline(myArray);
                var secondAnimated = L.animatedMarker(multiLineString.getLatLngs(), {
                                        distance: 100,
                                        interval: 500,
                                        iconSize:[16,16],
                                        iconAnchor: [7, 16],
                                        //autostart: false,
                                        icon: pulsingIcon
                                    });

                map.addLayer(secondAnimated);
                routeLayer = L.geoJSON(response)

            }

    })
    map.addLayer(routeLayer);
};


texts: Using Google Maps Coordinate System in QGIS?

texts:
I have some Points in Shapefile with BC Albers Coordinate Syetem and I need to Project them to Google Maps Lat and Long. I tried to do this by using

WGS 84 / Pseudo Mercator EPSG: 3857

and

WGS84 EPSG: 4326

But none of them returning correct result! I believe there is Coordinate for Google Maps in ArcGIS as "Web Mercator" to do this but I couldn't find it in QGis.

What Coordinate should I use for Getting correct points X and Y on google map?

I also used EPSG:900913 Google Mercator but when I tried to get X and Y coords through Geometry Tool/ Export add Geometry Columns I am getting numbers like :

Result using EPSG:3857:


texts:
I converted a shapefile to a SpatiaLite database so I could work with attributes with a length longer than 254 characters, and also not see some columns and change their order.

In QGIS 2.18.16 I have used DB Manager to connect to this SQLite database, selecting the columns from the table I want to bring into my project, and in the order I want to work with the columns.

I am using the following code in the SQL window to do this

    `SELECT  id, cpt, sub_cpt, name, area_hect, category, designat, grid_ref,
     species, plant_year, activity, actvtynote, years, restock, notes, GEOMETRY
     From "clientname subcompartments SQLITE "

It loads into the project fine, but the layer is not editable? I have the SQLite file on my local hard drive

How to load it as an editable layer?

I am sure it must be something simple that I am not doing, or just basic misunderstanding as I am not familiar with working with SQLite.

I have had a look through GIS Stack Exchange but cannot find anything.


texts: QGIS API (PyQGIS) doesn't properly read layers

texts:
I have a problem with reading qgs file by PyQGIS lib. I'm creating python standalone stript which externally use PyQGIS from installed QGIS desktop application. The problem is that PyQGIS sometimes read qgs file but WITHOUT READING LAYERS INSIDE.
I'm using Ubuntu 18.04 LTS, django, python 3.6.9, QGIS 3.4.14-Madeira

PyQGIS helps me with reading and updating .qgs files. I'm always using same start and ending code.

from qgis.core import QgsProject, QgsApplication

app = QgsApplication([], False)
app.initQgis()
app.exit()
project = QgsProject.instance()
project.read("/path/to/file.qgs")

(PyQGIS operations)

project.write()
app.exit()

And... the problem is that this code works fine... but never know how long, sometimes few days or few weeks but finally comes the moment when this code stops working properly. It still reads .qgs file but without layers inside, see groups but not layers. It should be impossible that code works fine for long time but some day it stops.

PyQGIS unfortunately works that way. The only thing comes to my mind is that the problem may be inside qgis desktop because PyQGIS is not common library. I belongs and works only with qgis desktop. This program maybe remember reading states or use multithreading with the way which i don't understand.

Of course, the problem may be completely different but i really don't have any more ideas how to change this code.

Does anyone know how to fix it?

Some logs from app which makes me wonder, maybe this will help.

WARNING: QApplication was not created in the main() thread.
Loading a file that was saved with an older version of qgis (saved in 3.4.13-Madeira, loaded in 3.4.14-Madeira). Problems may occur.
QObject::~QObject: Timers cannot be stopped from another thread
QObject::~QObject: Timers cannot be stopped from another thread
QObject::connect: Cannot queue arguments of type 'QDomDocument'
(Make sure 'QDomDocument' is registered using qRegisterMetaType().)


texts: How to harmonize two raster of different spatial resolution?

texts:
To harmonize two raster of different spatial resolution, what is most preferred and scientific to do:

To make coarse resolution finer

or

To make fine resolution coarser?

The datasets includes:

Landsat derived NDVI data at 100 m resolution and
Climate datasets on 'annual mean temperature' and 'Annual precipitation' at about 4 km resolution.

Goal: To develop some relationship between NDVI and Climatic variables


texts:
Let's say I have 3 layers:

point layer (columns: id, group),
table layer named materials (columns: point_id, mat_id),
table layer named mat_list (columns group, mat_id).

There is relation between point layer and materials table using id and point_id and value relation materials - mat_list.

In mat_list I grouped materials for example:

What I want to do is after I pick group in point layer it will limit choices to pick in materials layer.

EDIT:

More information maybe usefull

table/layer with material list

layer with points.

Relation between point and material.

In this screen I want to limit materials to group_2 so drop down list will show material_4, material_5, material_6. But group is picked in point layer not material. Is there any easy way to accomplish it.


texts: Restrict only a layer instead of giving the whole service?

texts:
Scenario (example):

I have a database called MYDB linked to GeoServer. In GeoServer I have a workspace called A where A has a store connected to my MYDB. In section layers I published four layers called waters, buildings, roads and people.

How can I provide an URL for each published layer separately.

For example, for layer waters I will have the following WMS URL:
http://localhost:8080/geoserver/geonode/wms?service=WMS&version=1.1.0&request=GetMap&layers=A:waters

When I connect this URL into QGIS it will connect to the whole workspace.

How can I ensure that it only will display one specific layer or group of layers?

Note that if I call in OpenLayers it will filter only the specific layer but in QGIS it will connect with the entire workspace.


texts: How to get QGIS layer information saved in JSON format?

texts:
I have called a URL to load imagery. Once the imagery is loaded in QGIS, I want the information of this layer (metadata), as seen in the below picture, to be stored in JSON format.
NOTE: It can be opened by right clicking the imagery layer --> Properties --> Information.

I have tried the below code but it returns QgsLayerMetadata object with null values:

a = iface.activeLayer().metadata()

How to save this information in JSON file?


texts: New feature class with matching records ArcGIS

texts:
I have two feature classes in a FGDB each with a set of time stamped points. Let's call them fc1 and fc2. Between these two feature classes one or the other does not consistently have a greater number of points, sometimes fc1 will have more and other times fc2 will.

I want to create two new feature classes (fc1_join and fc2_join) that only contain matching records based on the timestamp. In the end fc1_join and fc2_join should have an equal number of features, correct?

And what steps would be taken to do this?


texts: Selected feature disappears when selected again

texts:
I am having problem editing line layer with editor widget

If I click the selected feature for editing vertices  the
selected features disappears sometimes
If I select the graphics
and then select new selection from toolbar and select the already
selected feature  it disappears.

I can't understand why its happening, below is jsfiddle code of the Editor widget with a polyline layer.

https://fiddle.jshell.net/b5a20ua8/


texts: Using PyQGiS in custom application gives import qgis.core error

texts:
I'm using windows 7 and/or 8 64 bit with QGiS 2.01.  I'm wanting to us the pyQGiS tools in a custom python application outside of the QGiS environment (i.e by running a python executable located in C:/Python27.

To do this I have been following the instructions posted at http://www.qgis.org/en/docs/pyqgis_developer_cookbook/intro.html

For the first step (to remove the error "ImportError: No module named qgis.core") i set my pythonpath to "C:\Program Files\QGIS Dufour\apps\qgis\python"

As expected I now get an error:
File "C:\Program Files\QGIS Dufour\apps\qgis\python\qgis__init__.py", line 26, in
    import sip
ImportError: DLL load failed: The specified module could not be found.

The instructions say to add "the directories where the QGIS libraries reside".  I have tried adding a number of paths, with no success.  I also cannot locate the sip.dll file anywhere within the C:\Program Files\QGIS Dufour directory.

Does anyone know which paths I need to set?


texts: Different DEM range values in QGIS and ArcGIS

texts: I downloaded a raster image from USGS  for a part of Malaysia. I opened it by QGIS and ArcGIS Desktop as shown in the attached photo. The first problem is that some values are in a negative and the second is the elevation range in QGIS is different compared with ArcGIS Desktop.


texts: problem transforming CRS coordinates EPSG:4269 to EPSG:3857 in python

texts:
source_srs = get_srs(shp_abs_path)
source_epsg = source_srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
print source_epsg #prints 4269
target_epsg = 3857 # because these are headed to leafletjs
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(target_epsg)
osr_transform = osr.CoordinateTransformation(source_srs, target_srs)

then as I enumerate through the shapefile I transform the bboxes on the fly

sf = shapefile.Reader(base_path)
shapeRecords = sf.shapeRecords()  #will store the geometry separately
self.shapefile_records[fn] = shapeRecords
for i, shape in enumerate(sf.shapes()):
    bbox = shape.bbox
    ## OGR/SRS TRANSFORMATION
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(bbox[0],bbox[1])
    point.Transform(osr_transform)
    print point.GetX(), point.GetY()

    point.AddPoint(bbox[2],bbox[3])
    point.Transform(osr_transform)
    print point.GetX(), point.GetY()

which prints coordinates like this...

-13613708.5846 4623834.1438
-13613668.296 4623869.60132

these put me out in the ocean I think, but they're supposed to be in Napa county. Any suggestions are warmly welcomed. I can use almost any python lib.

NEW

I've identified an issue and I'm not sure if it's with osgeo or the 4269 to 4326 transformation I've defined. Have a look at this output:

>>> from osgeo import ogr, osr
>>> source_srs  = osr.SpatialReference()
>>> source_srs.ImportFromEPSG(4269)
0
>>> target_srs = osr.SpatialReference()
>>> target_srs.ImportFromEPSG(4326)
0
>>> osr_transformation = osr.CoordinateTransformation(source_srs, target_srs)
>>> point = ogr.Geometry(ogr.wkbPoint)
>>> point.AddPoint(6449750.845, 2001628.41312)
>>> point.Transform(osr_transformation)
0
>>> print point.Transform(osr_transformation)
0
>>> point.GetX(), point.GetY()
(6449750.845, 2001628.4131199997)

Nothing is modified. I've triple confirmed the data is in EPSG:4269 and included a shortened bitly link to it below in the comments.


texts:
I am trying to use the OpenLayers TimeLine library with a kml file.

I used the  available example for kml in github, but it shows points instead of polygon.

How can I display polygon the TimeLine library of OpenLayers?


texts:
I am trying to ingest a point shapefile to geomesa, I have added the schema and convertor part to conf file but when i am giving the ingest command usinf geomesa-fs it is giving me the below error.

**

  ERROR Warning: Missing dependency for command execution: Could not
  initialize class
  org.locationtech.geomesa.fs.FileSystemDataStoreFactory$

**

Anyone got any idea how it can be solved?


texts:
Actually, I want to vector bend a layer (not properly georeferenced), so it fits in my project.
When I add it in my QGis 3.14, it appears on a very strange location, i.e. its not even on planet earth.
Now, when I try to use vector bending, QGis crashes everytime and I think its happening because of the extreme distance, but I truly dont know about that.
Is there a possibilty to just add the new Layer right on top of my current projectlayers to use vector bending properly?
The first picture is the layer I want to add, the second my actual project. On the bottom you can see the coordinates


texts: MapBasic For Loop ignoring every other row

texts:
I have made a loop in MapBasic but is it only updating every other row in the table:

My code works from the bottom of the table upwards to generate start and end coordinates for polyline features:

    Dim nrows as integer

Dim USX, USY, UEX, UEY as float

For nrows = TableInfo(WorTab, TAB_INFO_NROWS) to 1 Step - 1
    Select * from WorTab where rowid = nrows
    Dim ObjTYPE as integer
    ObjType = ObjectInfo(Selection.obj, OBJ_INFO_TYPE)
    If ObjType = 4 then
    nPoints = ObjectInfo(Selection.obj, OBJ_INFO_NPNTS)
    USX = ObjectNodeX(Selection.Obj,1,1)
    USY = ObjectNodeY(Selection.Obj,1,1)
    UEX = ObjectNodeX(Selection.Obj,1,nPoints)
    UEY = ObjectNodeY(Selection.Obj,1,nPoints)
    Update Selection Set SX = USX, SY = USY, EX = UEX, EY = UEY
    nrows = nrows -1
    Close Table Selection
    Else
    End if

I'd be very grateful if somebody could point out what I'm doing wrong and why nrows -1 is the wrong choice.


texts: How do you display latitude \ longitude coordinates in ArcGIS?

texts: I would like to display latitude and longtitude coordinates (at the bottom of window) for my data which I am editing in different coordinate system. Is it possible?


texts:
I am trying to use the Fiona module in Python, but am having trouble getting it to work. After issuing conda install fiona, the installation using conda seems fine, but when I issue import fiona, I get the following error message: module 'fiona' has no attribute '_loading'
Has anyone else encountered this problem? The version of Fiona that I have tried to install is 1.18.18.
It seems related to the issue identified here, but I am not exactly sure what this means. I am running conda in an ArcGIS Pro environment, and so don't want to create any problems.


texts:
I am trying to add a clickable button on Leaflet map. Right now using this:
<ion-content padding scrollbar-y-auto>

    <button ion-button class="button-action" (click)="Action()"  block>Action!>/button>

</ion-content>

CSS:
page-map {
  @font-face {
    font-family: 'Minecraft';
    src: url('../assets/Minecraft.ttf')
    }
    .button-action{
      height: 40px;
      width: 200px;
      text-decoration: none;
      text-align: center;
      color: white;
      background: url('https://i.hizliresim.com/QpVGbg.png') no-repeat;
      font-family: Minecraft;
      line-height: 26px;
      padding-top: 6px;
  }
  .button-action:hover{
      background: url('https://i.hizliresim.com/a1zraR.png') no-repeat;
  }
}

This is what I got:

I can't click the button it's like a picture also hover option is not working too.
class="leaflet-bottom leaflet-right"> killing the button, or somehow button goes behind map so I can't click


texts: Reclassifying vector polygons based on text string to integer using Python Parser of ArcGIS Field Calculator?

texts:
So I am trying to reclassify a text string by searching for bits of strings and converting those strings to an integer value in a new field within the field calculator. e.g. I want to find all values in the field "text" that contain the word "cover" in them, for example, "Other deciduous trees with 11- 25% impervious cover". I then want to take that and assign an integer value to a new field "class". I tried to write the script (shown below), but no dice so far, anyone have any insight?

def reclass(C_TEXT):
if "cover" in C_TEXT:
    return 10


texts: Using Raster calculator in Modelbuilder ArcGis 10

texts:
I am using Modelbuilder to calculate several rasters with different parameters for each grid cell. One of the final steps is to calculate the log of one of the rasters using Raster calculator. Unfortunately, I get the error at that last final step:

  ERROR 000539: Error running expression: rcexec()  Traceback (most
  recent call last):   File "", line 1, in    File
  "", line 5, in rcexec TypeError: cannot concatenate 'str' and
  'int' objects

  Failed to execute (Raster Calculator).

My model looks that way:

I use the raster calculator to compute the ln

I can imagine that one of the reasons is that the arrow to the "raster calculator" is not solid. But I dont know why. The rasters derived from the previous step are fine. When using raster calculator outside modelbuilder to do that last step, all works fine. I am a bit at a loss here.


texts: Python scripts that run inside ArcMap vs ones that run outside?

texts:
I'm just starting to get into Python scripting for work.

I'm currently creating a script to automate a process.

Basically, it asks the user the client name, gets a projection if available, creates directory on C: drive for client, creates file geodatabase specific to client, creates required dataset, and creates feature classes specific to clients data. Eventually, it will also add the required fields to each feature class and probably some other things.

I started this not really knowing the proper etiquette of Python scripting for ArcMap. But what I've created so far will only run outside ArcMap I believe.

Is this acceptable?

Instead of getting user input through arcpy.getparamaterastext(), which I just found out about, I'm using raw_input().

Is this okay?

It does work, I'm just not sure if this is a proper way to do scripting.

Here's the code I have so far.

import sys
import arcpy
import os

#Records name of the client
client = raw_input("Enter the name of the client: (letters and underscores only) \n")

#Records filepath of client to be created
clientpath = "C:/" + client

#Inquires if projection file exists
projection = raw_input("Is there a .prj or .shp available with correct projection? Y or N \n")

#Records the projection location if available
if projection.upper() == "Y":
    spatialr = raw_input("Drag the .prj or .shp here to record the filepath \n")
    nspatialr = spatialr.replace('"', "")
elif projection.upper() == "N":
    alert = raw_input("You must add the spatial reference manually, hit enter to continue. \n")
elif projection.upper() != "N" or "Y":
    exit = raw_input("That is not a valid response. Try again. \n")
    sys.exit()

#Checks if client folder exists; if not, creates one
if not os.path.exists(clientpath):
    os.makedirs(clientpath)

#Variable for file geodatabase location
FGBpath = clientpath + "/" + client + ".gdb"

#Checks if client file geodatabase exists; if not, creates one
if not arcpy.Exists(FGBpath):
    arcpy.CreateFileGDB_management(clientpath, client)

#Variable for dataset location
FDatasetpath = clientpath + "/" + client + ".gdb" + "/Network"

#Checks if dataset exists; if not, creates one
if not arcpy.Exists(FDatasetpath):
    if projection.upper() == "Y":
        arcpy.CreateFeatureDataset_management(FGBpath, "Network", nspatialr)
    elif projection.upper() == "N":
        arcpy.CreateFeatureDataset_management(FGBpath, "Network")

#Variable for cable feature class location
FCcablepath = clientpath + "/" + client + ".gdb" + "/Network" + "/cable"

#Checks if cable feature class exists; if not, creates one
if not arcpy.Exists(FCcablepath):
    if projection.upper() == "Y":
        arcpy.CreateFeatureclass_management (FDatasetpath, "cable", "POLYLINE", "", "", "", nspatialr)
    elif projection.upper() == "N":
        arcpy.CreateFeatureclass_management (FDatasetpath, "cable", "POLYLINE")

#Variable for splice point feature class location
FCsplicepath = clientpath + "/" + client + ".gdb" + "/Network" + "/splice_point"

#Checks if splice point feature class exists; if not, creates one
if not arcpy.Exists(FCsplicepath):
    if projection == 'Y' or projection == 'y':
        arcpy.CreateFeatureclass_management (FDatasetpath, "splice_point", "POINT", "", "", "", nspatialr)
    elif projection == 'N' or projection == 'n':
        arcpy.CreateFeatureclass_management (FDatasetpath, "splice_point", "POINT")

exit = raw_input("\n\n File geodatabase, dataset, and the cable \n and splice point feature classes successfully created. \n\n Hit enter to exit.")

I still have some work to do, like adding the needed fields.


texts: ArcGIS Pro - reset to default (reset profile)

texts:
Is anyone aware how one can set ArcGIS Pro back to default settings?
In ArcMap etc it was done by following this tutorial


texts: How to submit data to modify feature using applyEdits of ArcGIS API for JavaScript?

texts:
I'm new to ArcGIS API. In a form, how can you submit data to modify a feature using applyEdits? No map is involved with this.

I have a form that populates the editable data with javascript. When the user clicks save changes, I want the changes to be processed using applyEdits.

So far, I have:

var Editable = new FeatureLayer("......", {
        mode: FeatureLayer.MODE_SNAPSHOT
    });

    Editable.applyEdits(null,[ {
      "attributes": {
        "bldg_prop_id": 6086,
        "bldg_water_utility_name": "New Data"
      }
    }],null);

This doesn't work but I'm not sure where to go with this.


texts: Multiple on-the-fly filtering based on markers' features on leaflet

texts:
I would like to filter the markers on my (geojson) layer according to two features: year and type of event, ideally using both filters (year and eventType) as check-boxes without using layerGroup at all. I would like to be able to click on a checkbox and update dynamically the markers by filtering them according to feature.properties.eventType too, in addition to the year displayed, something like

  <input type="checkbox" id="all" name="all" value="all">
  All
  <input type="checkbox" id="funfair" name="funfair" value="newsletter">
  Funfair
  <input type="checkbox" id="trade" name="trade" value="trade">
  Trade
  <input type="checkbox" id="exhibition" name="exhibition" value="exhibition">
  Paid exhibition

etc... which would trigger a

filter: function(feature, layer) {
  return feature.properties.eventType == value-of-checkbox;
}

I have five types of eventType. So, for instance, if 1890 is checked, the map is showing all 1890 events, no matter what kind of event it is (because by default it would show all events). If the user checks the 'funfair' and the 'trade' checkbox, though, I want the map to show only the 1890 events of the funfair and trade kind, and not other kind of events.

This is my complete script, which so far is working fine with filtering by year.

var oldmap = L.tileLayer('address', {
  maxZoom: 12
}),

var modernmap = L.tileLayer('https://api.tiles.mapbox.com/v4/{id}/{z}/{x}/{y}.png?access_token={accessToken}', {
  maxZoom: 12,
  id: 'mapbox.streets',
  accessToken: 'mytoken'
});

var map = L.map('map', {
  center: [40.1857864,-4.5536861],
  zoom: 7,
  layers: [oldmap, modernmap]
});

var baseMaps = {
  "Old map": oldmap,
  "Modern map": modernmap
};

var jsontest = {"type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates":  [ -16.5471268,28.4136726 ]
      },
      "properties": {
        "Name":"Point1",
        "year":'1892',
        "eventType":"funfair",
        "notes":""
      }
    },
    {
      "type": "Feature",
      "geometry": {
        "type": "Point",
        "coordinates":  [ 2.1753895,41.3767925 ]
      },
      "properties": {
        "Name":"Point2",
        "year":'1890',
        "eventType":"trade demonstration",
        "notes":""
      }
    }
  ]};

var mapLayerGroupsYear = [];

function onEachFeature(feature, layer) {
  var year = mapLayerGroupsYear[feature.properties.year];
  if (year === undefined) {
    year = new L.layerGroup();

    //add the layer to the map
    year.addTo(map);

    //store layer
    mapLayerGroupsYear[feature.properties.year] = year;
  }

  //add the feature to the layer
  year.addLayer(layer);
}

var myLayer = L.geoJSON(jsontest, {
  onEachFeature: onEachFeature
})

L.control.layers(baseMaps, mapLayerGroupsYear).addTo(map);

I am a beginner with leaflet so please forgive the dumb question, in case it is.


texts: ST_Distance in Virtual Layer query not working on QGIS on Mac but on Windows

texts:
I tried to use solution posted here, but ST_Distance in the Virtual Layer query doesn't work properly in Qgis (v3.2) on Mac, the result is showing blank while the same script is working properly in QGIS (v3.4) on Windows.

Here is my test script for virtual layer:

WITH distances AS(
  SELECT
    s.id site_id,
    a.Mem_ID mem_id,
    ST_Distance(
            s.geometry,
            a.geometry
            ) dist
  FROM site_a27 s
  CROSS JOIN A27_re a

)
SELECT * FROM distances

table for site_a27:

id,name,Response_X,Response_Y
1A27,bade,121.301412,24.968988

table for A27_re:

Mem_ID,Response_X,Response_Y
19960126000178,121.312175,24.987227
19960129000943,121.302705,24.978186
19960131001647,121.303628,24.957379
19960131002543,121.300855,24.984292
19960203000247,121.296974,24.968484

output of virtual layer on Mac:

output of virtual layer on Windows:

Information about my QGIS version on Mac:

How do I solve this in QGIS on Mac?


texts: Making a raster coverage that shows the average speed of points

texts:
I had a layer of points of trajectory and I wanted to create a map of the average speed of the movement. I have a column of speed in the table of points.
I created a big base raster coverage. Then, I used ST_PixelAsPolygons() to create a polygon of each pixel. Then, I calculated the average speed of points in each pixel. Finally, I used ST_AsRaster() to create a raster coverage of the speed map.

I used the following code (where passages is the point table and grid_pixel is the polygon of each pixel and sog(speed over ground) is the speed):

CREATE TABLE speed_map AS(

with speed_average AS (
SELECT g.id , AVG(p.sog) AS sog_avg
FROM grid_pixel g
INNER JOIN passages p
ON ST_Within(p.geom_point, g.geom)
GROUP BY g.id
),

speed_avg_geom AS (
SELECT g.id AS id, g.geom AS geom, sa.sog_avg
FROM speed_average sa
INNER JOIN grid_pixel g
ON g.id = sa.id)

SELECT id AS id , ST_AsRaster(geom, 1, 1, '8BUI', sog_avg) AS rast
FROM speed_avg_geom
)

Now, each small polygon of grid_pixel table is transfered to a raster tile (one row of raster table) and I have a very big table that has one pixel in each row.
The map is still correct and I can see the correct raster in QGIS. However, I think it is better to have tiles of 100*100 pixels.

I wanted to retile the raster and I found ST_reTile() function but I do not know how to use it to retile the raster to 100*100 rows.

Or maybe even my method is totally wrong from the beginning where I created a polygon for each pixel.
Can anyone please let me know what is the best way to create a raster of average speed of points.


texts: Converting a raster pixel by pixel to vector cells in Python

texts:
Problem

Today I needed to convert raster to a shapefile with rectangular polygons representing each pixel of the raster in Python, including the pixels with value 0.

Related questions:

Converting raster pixels to vector cells? does not consider Python solution.
Convert raster to vector based on each pixel offers a solution in SAGA and OpenJUMP
How to polygonize raster to shapely polygon does not polygonize zero pixels and if I understand it correctly, it groups the shapes by values - I explicitly want to have each recangular cell/pixel as a separate feature.

Solution attempt

I could not find an existing solution and so I wrote my own hacky solution using rasterio, shapely, and geopandas.

I am assuming that I know the raster file is in epsg:4326.

First import the libraries

import itertools
import rasterio
from shapely.geometry import box
import geopandas as gpd

Create rectangles using shapely

with rasterio.open(sample_raster) as dataset:
    data = dataset.read(1)

    t = dataset.transform

    move_x = t[0]
    # t[4] is negative, as raster start upper left 0,0 and goes down
    # later for steps calculation (ymin=...) we use plus instead of minus
    move_y = t[4]

    height = dataset.height
    width = dataset.width

    polygons = []
    indices = list(itertools.product(range(width), range(height)))
    for x,y in indices:
        x_min, y_max = t * (x,y)
        x_max = x_min + move_x
        y_min = y_max + move_y
        polygons.append(box(x_min, y_min, x_max, y_max))

Extract the data

data_list = []
for x,y in indices:
    data_list.append(data[y,x])

Combine in a GeoDataFrame using geopandas

gdf = gpd.GeoDataFrame(data=data_list, crs={'init':'epsg:4236'}, geometry=polygons, columns=['value'])

Question

This is a quick solution and I suspect not ideal. I have two questions

1) Is this solution correct? Are there some obvious or necessary improvements?

2) Is there a quicker / standard way to do this? Perhaps I have missed some library, or a function.


texts:
I have a problem where I need to, for each cell, calculate the distance to the nearest flow channel. This is likely a complex problem due to the fact that I need the nearest flow channel while considering topography (flow direction). A simple Euclidean distance does not suffice due to this. I cannot even think of terminology to use to google a little deeper. Anyone with any experience they would like to bestow on to me?

ArcGIS 10.3, all license levels

So, I have looked at both path distance and flow length and, correct me if I am wrong, neither work for my current scenario. The path distance tool does not look like it incorporates flow direction and the flow length tool only calculates distance to the last cell in the drainage watershed. I am looking for the least cost path sum of distance for each cell to the nearest source point considering flow direction. Am I missing something in either method?

I have now been looking at what @dru87 pointed to and am making progress. The method of subtracting the dem from 0 and computing flow accumulation upon this raster is almost what I need. There is just one further step that I need to address; the fact that I have defined a classification of the FA to represent stream channels and I need the distance from this classified layer. I feel this necessity might limit this method as editing of the FD layer would likely be required. Thoughts?


texts: How many random points are needed to statistically describe a DEM?

texts:
I'm running some statistical comparisons between several DEMs (SRTM, ASTER, ALOS...) through GRASS-GIS and Python. Dealing with millions of cells is not very practical, so I thought I could compare elevation values from random points. How many points should I use? 5000? 10000?

Of course, the more, the better. But then I went searching for any literature on the subject and couldn't find anything with a good answer.

Can someone point me a good article on this?


texts: Avoiding thin parts of polygons using ArcGIS Desktop?

texts:
I've used Simplify Polygon with 'Maximum Allowable Offset' 250m, and then Aggregate Polygons with 'Aggregation Distance' also 250m (opposite order also tested). Still I have thin parts of polygons as below:

I presume that I need some minimum angle tool.
But there are also thin inner parts like this:

Yes, the width is only 35m, while I want objects have at least 250m. In this case propably 'Split Narrow Points' would be helpful, but I don't have this tool.
Any ideas please?


texts:
When I try to use v.proj to reproject my map from OSGB to GoogleMaps coord system I get this error....

Reprojecting primitives:
pj_transform() failed: failed to load datum shift file
ERROR: Error in pj_do_transform

Can someone confirm I not doing something stupid to cause this error please ...

My file structure is...

e:\work\grass\GB\mapset1
e:\work\grass\UK_EPSG900913\mapset1

I.e. 'mapset1' folders contain the 'vector' and 'dbf' directories etc.

So I have two locations 'GB' and 'UK_EPSG900913'
The projections are set to OSGB and epsg900913 respectively.

I have a map in the 'GB' location that I wish to reproject to epsg900913
I understand that to do this I choose the location as 'UK_EPSG900913' and 'mapset1' when I launch GRASS
and then I type this command...

v.proj input=scot location=GB mapset=mapset1 output=scot1

And then a map should appear in my 'UK_EPSG900913/mapset1/vector' folder.

Is my understanding of how to re-project correct? If so can anyone suggest what might be going wrong?

I'm on Windows 7 64-bit, GRASS 6.4.2


texts: Adding * .prj file to output via Python script tool that creates fishnet?

texts:
I'm a new python user.

I wanted to make a tool box like fishnet in ArcGIS 10.2.

Everything, works fine. However, at output there is no .prj file, although spatial reference was defined at the tool.

How can I modify the script to get a .prj file in output according to defined spatial reference?

I used following function to define spatial reference and finally used fishnet to make the grid.

spatial_ref = arcpy.GetParameter(2)

arcpy.CreateFishnet_management(outFeatureClass, originCoordinate, yAxisCoordinate, cellSizeWidth, cellSizeHeight, numRows, numColumns, oppositeCoorner, labels, templateExtent, geometryType)


texts: KML from PostGIS Geometry into Openlayers

texts:
im getting the KML value from my postgis database.. st_askml()

And the output is like this:

"<Polygon><outerBoundaryIs><LinearRing><coordinates>-99.999999990686774,48.499999959021807 -99.999999990686774,47.999999988824129 -99.500000020489097,47.999999988824129 -99.500000020489097,48.499999959021807 -99.999999990686774,48.499999959021807</coordinates></LinearRing></outerBoundaryIs></Polygon>"

Now i need to display them on a Openlayers Map i got.

Layer:
var KML_Layer = new OpenLayers.Layer.Vector('KML');

Formatting code:

function GetFeaturesFromKMLString (strKML) {
var format = new OpenLayers.Format.KML({
        'internalProjection': options.projections.source, //EPSG:3785/900913
        'externalProjection': options.projections.display //EPSG:4326
    });

return format.read(strKML);
};

Adding the feature to the layer:

KML_Layer.addFeatures(value.st_askml);

Debuggin info:

value.st_askml = The output format i pasted earlier in topic.

After running it trough the GetFeaturesFromKMLString() its set to [] (No value here)

How can i get this to display correctly on the map?


texts: Counting points along polyline at given distance of every 40 Ft using ArcMap?

texts:
I've got a polyline representing a pipeline and have a pig run that has metal loss.

I am going to create a points shapefile in ArcMap 9.3.1 of the metal loss locations along the pipeline but want to do it like a sliding rule every 40 ft., giving me a count of the metal loss points per 40' section.

What would be the best way to perform this task in ArcMap?


texts: Error 000728 field does not exist within table - field already exists

texts:
I am running arcpy.AddField_management to add a field called "cnid" and update it with x centroid values.
For add.field I get the error "cnid" already exists. I tried print.field name and I cannot see "cnid" field. Now I am running calculate shape centroid and I get the error message "cnid" does not exist within table. I cannot understand, sometimes it works (if I close arcmap and open it again) but if I run both codes, add field and update it, it doesnt work.

    # AddField for new fields
arcpy.AddField_management("Connections","cnid","DOUBLE", "", "", "", "", "NULLABLE")

    #test fields
    for field in arcpy.ListFields ("Connections"):
    ...     print field.name

    xExpression = "float(!SHAPE.CENTROID!.split()[0])"
    inFeatures ="Connections"
    fieldName ="cnid"
    arcpy.CalculateField_management ("Connections","cnid",xExpression, "PYTHON")

I read something about in_memory and refresh bug from Why am I getting an ERROR 000728: Field *fieldname* does not exist within table even though field exists?, but how can I refresh it without having to close and open it again?

I am using ArcMap's Python window.


texts: Conversion of jpg image to greyscale and how to get attribute table?

texts: I am studying on concrete technology, I have jpg image of concrete which is taken by digital camera .I am working with arcgis 10.1 version. I need to convert it into greyscale and get attribute table for further analyasis and classification.How can I get that ?


texts: ArcGIS local server using javascript api still does http gets to services.arcgisonline.com (internet)

texts:
We have setup a local server hosting the ArcGIS server software (10.1).  We installed the ArcGIS JavaScript API on the same local server.  Our goal is to supply maps when the client has no internet connection, only a local area connection.
When we run the client (i.e. using the JavaScript API from the local server), we can see that all of the java script files are pulled from the local server.  However, once our java script code calls new esri.Map, we see that the javascript code starts to pull map data from http://services.arcgisonline.com.
I looked at the C:\ArcGIS\JavaScriptAPI33\arcgis_js_v33_api\library\3.3\jsapicompact\init.js code (this is where it resides on our local ArcGIS server) and I see many occurrences of http://services.arcgisonline.com.  I am thinking I need to change all of these occurrences to be the local server but this was not a step in the installation instructions.
Does anyone have any advice on how we can prevent the javascript from going to the internet for any files (map data, javascript data, css, etc).  I understand that the javascript is the controlling code that makes this happen but I thought that esri would have handled that in their javascript logic.
Edit
Thanks for the input.  I did follow only the steps given in the readme (I have changed nothing else).  As I said, all the java script files come from the local server.  Once I call "new esri.Map" in our java script code, I see the following "http gets" to services.arcgisonline.com:
 1. services.arcgisonline.com/ArcGIS/rest/info?f=json
 2. services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer?f=json&callback=dojo.io.script.jsonp_dojoIoScript1._jsonpCallback
 3. services.arcgisonline.com/ArcGIS/rest/services/Reference/World_Boundaries_and_Places/MapServer?f=json&callback=dojo.io.script.jsonp_dojoIoScript2._jsonpCallback

Also, you asked me "do we have any basemap or anything else in our map that comes from ArcGIS online?"  My answer is no.  But clearly (unless there is an ArcGIS bug) something is not working as we want it to because it ends up going to services.arcgisonline.com.
Thanks again for the response, we really appreciate the help.
As your response is implying, it seems I do not understand the basemap values since I am using hybrid for the local server.  Sounds like I need a different value there.
Once again, thank you very much for the help.
Here is the code snippet:
var l_Options =
{
zoom: 10,
sliderStyle: "small",
basemap: "hybrid",
infoWindow: g_Popup,
nav: true
};

map = new esri.Map("map_canvas", l_Options);
dojo.connect(map, "onLoad",
function()
{
    var l_ArcGISTiledMapServiceLayer = null;

    switch (g_strMappingEngine)
    {
    case 'esriLocalServer':
        /* ????
        l_ArcGISTiledMapServiceLayer =
        new esri.layers.ArcGISTiledMapServiceLayer(
        l_strMappingDataAddress
        );
        ???? */
        l_ArcGISTiledMapServiceLayer =
        new esri.layers.ArcGISDynamicMapServiceLayer(
        l_strMappingDataAddress
        );
        break;
    case 'esriInternet':
    default:
        l_ArcGISTiledMapServiceLayer =
        new esri.layers.ArcGISTiledMapServiceLayer(
        l_strMappingDataAddress
        );
        break;
    }

    map.addLayer(l_ArcGISTiledMapServiceLayer);
    ...