texts: Separating a shapefile with specific features into a new shapefile using PyQGIS

texts:
I have a shapefile. I would like to know if it is possible to classify the different layers in that shapefile (like highways, local roads, rail lines, etc.) with MTFCC and extract a separate shapefile from the existing one.

So for example, if I had a shapefile of our existing county with all vectors present could I make another, separate shapefile, of that same county containing only highways and local roads?

This is all using PyQGIS v2.14.0


texts: Exporting shapefiles with symbology using ArcGIS Pro for use in ArcGIS Runtime

texts:
I have a set of shapefiles with elevation values in its attribute. In ArcGIS Pro, I have placed these shapefiles into a Scene as a 3D layer. Then I have added extrusion, and the shapefiles are rendered correctly with height.
I tried to load the shapefiles into a Scene in ArcGIS Runtime, and did the extrusion in code. However, the symbology (or the color) is lost. I am aware that I can hard-code each shapefile's symbology in code, but that sounds wrong to me.
Is there a way for me to export these from ArcGIS Pro and use them in ArcGIS Runtime without hard-coding anything?
Update
This is what I am currently doing in code:
private async void LoadShapefiles(Scene sc) {
    foreach (var file in GetPaths()) {
        ShapefileFeatureTable myShapefile = await ShapefileFeatureTable.OpenAsync(file);

        FeatureLayer layer = new FeatureLayer(myShapefile);
        var symbol = new SimpleFillSymbol(SimpleFillSymbolStyle.Solid, Colors.Gray, null); // Hard-coded
        var renderer = new SimpleRenderer(symbol);

        renderer.SceneProperties.ExtrusionMode = Esri.ArcGISRuntime.Symbology.ExtrusionMode.Maximum;
        renderer.SceneProperties.ExtrusionExpression = "[HEIGHT]";

        layer.Renderer = renderer;
        layer.RenderingMode = FeatureRenderingMode.Dynamic;

        sc.OperationalLayers.Add(layer);
    }
}

private string[] GetPaths() {
    return Directory.GetFiles(@"D:\Test", "*.shp");
}

This causes all features to be rendered with the same color.


texts:
I'm basically trying (and failing) to replicate this in OpenLayers 3.

I've set up a fiddle here:
http://jsfiddle.net/vu8ecy8v/5/

I'm under the impression that OpenLayers will reproject a GeoJSON source 'on the fly' if the projection of the GeoJSON is specified - which I'm doing in layers: of ol.map:

    ...
    new ol.layer.Vector({
        title: 'Montana',
        source: montanaSource,
        projection: 'EPSG:4326',
    }),
    ...

This makes no difference. Am I missing something obvious or do I have to reproject the GeoJSON in some other way.

If I change the background / view to be EPSG:4326 (same as the source) the GeoJSON appears without a hitch. These alternatives are in the fiddle (commented out.)


texts: How do I figure out the grid references (Eastings and Northings) necessary to create a circular 1 mile radius around one location?

texts:
How do I figure out the grid references necessary to create a circular 1 mile radius around one location?

For instance:

Grid Reference: SO 47904 87484

X (Eastings):   347904

Y (Northings):  287484

I know the approach has something to do with referencing the number of degrees from the original point for instance 10. Then subtracting the output by the number of metres in the radius size (a mile, so 1609.344).

But I've forgot the actual approach. What is it, thanks?


texts: Visualizing frequency of pixel change between several rasters using ArcMap

texts:
I have several rasters at my disposal. They come from a digitalization work according to a pre-established typology.
I would like to visualize the frequency of change of each pixel over time. Let me explain.

Let's say I have 3 rasters corresponding to 3 dates (same area; same resolution, same typology):

If a pixel changes (of typology class) 3 times over time it obtains the value 3
If a pixel changes (of typology class) 2 times over time it obtains the value 2
If a pixel changes (typology class) 1 time over time it obtains the value 1
If a pixel does not change the typology class over time, it obtains the value 0

Do you have an idea / solution to carry out this manipulation using ArcMap?


texts: OpenLayers 5 multiple markers all end up at last position

texts:
I have a LineString in OL5 and I want to iterate over each segment and display a label at the segment end. I have this mostly working, but all the Overlay labels magically end up at the last position inserted.
I have this code:
function createMeasureTooltip() {
    let measureTooltipElement = document.createElement('div');
    measureTooltipElement.className = 'tooltip tooltip-measure';
    let measureTooltip = new Overlay({
        element: measureTooltipElement,
        offset: [0, -15],
        positioning: 'bottom-center'
    });

    return measureTooltip;
}

draw.on('drawend',
    function(evt) {
        let out = '';
        let measureTooltipElement;
        let geom = evt.feature.getGeometry();
        console.log('In drawend');
        console.log(evt);
        if (geom instanceof LineString) {
            geom.forEachSegment(function (a, b) {
                out = formatSegment(a, b);
                let measureTooltip = createMeasureTooltip();
                measureTooltipElement = measureTooltip.getElement();

                measureTooltipElement.innerHTML = out;
                measureTooltipElement.className = 'tooltip tooltip-static';
                measureTooltip.setOffset([0, -7]);
                measureTooltip.setPosition(b);
                console.log(b);
                console.log(out);
                console.log(measureTooltip);
                map.addOverlay(measureTooltip);
                measureTooltip = null;
            });
            measureTooltipElement.innerHTML = out + ', (' +
                formatLength(geom) + ')';
        }
        // unset sketch
        sketch = null;
        $('.tooltip-measure').remove();
        unByKey(listener);
    }, this);

the drawend code is from my sketch Interaction draw handler.
I get all the labels correctly generated, but they do not stay where they are initially positioned at, and are all located at the end of the LineString.
Any ideas on how to fix this? I've run out of ideas and nothing seems to help.


texts: CartoDB Torque script causing JQuery UI conflict

texts:
I am using the torque library hosted here: https://github.com/CartoDB/torque

I am using it in a viz.json-based visualization (so the torque layer is defined in the cartodb web console.) When I try to include a torque layer when I have imported JQuery, the visualization freezes with the following stacktrace:

 Uncaught TypeError: undefined is not a function          cartodb.mod.torque.js:6
 cdb.geo.ui.TimeSlider.cdb.geo.ui.InfoBox.extend.updateSlider               cartodb.uncompressed.js:676
 Backbone.Events.trigger           cartodb.mod.torque.js:5
 L.undefined.L.TorqueLayer.L.CanvasLayer.extend.setKey       cartodb.mod.torque.js:5
 (anonymous function)      cartodb.mod.torque.js:3
 r.time      cartodb.mod.torque.js:3 r._tick

The torque visualization only works if I remove my import of JQuery, which breaks other parts of the website, such as the tabs() function.

Below are the imports I am using, the first works for Torque but breaks tabs(), the second works for tabs() but breaks Torque.

For Torque:

 <script src="dist/torque.uncompressed.js"></script>
 <script src="http://code.jquery.com/ui/1.11.1/jquery-ui.js"></script>
    <script src="http://libs.cartocdn.com/cartodb.js/v3/cartodb.uncompressed.js"></script>

For no Torque:

    <script src="dist/torque.uncompressed.js"></script>
    <script src="http://code.jquery.com/jquery-1.10.2.js"></script>
   <script src="http://code.jquery.com/ui/1.11.1/jquery-ui.js"></script>
    <script src="http://libs.cartocdn.com/cartodb.js/v3/cartodb.uncompressed.js"></script>

The tabs error occurs at:
$('#windowtabs').tabs();
Uncaught TypeError: undefined is not a function

You should be able to simulate this problem by trying to load the different sets of libraries as I posted -- on the one hand the torque visualization will freeze (if jquery is loaded), on the other hand, the torque visualization will work, but .tabs() will create an error.


texts: Scraping the vectorlayer/shp/SVG/GeoJSON from this webmap?

texts:
I try to download the data from this webmap for further analysis. I tried to look at the network tab, but didn't find any such data. I looked into the JavaScript files, and didn't find any either. Yet, I found some array data calling such as "LAYERS" on the JavaScript file, but I can't find any URL to the files.

Link: https://www.atrbpn.go.id/Peta-Bidang-Tanah

Specifically I want to get the "Zona Nilai Tanah" layer, but I think it's not a separated data from other layer.


texts: Moving end point of multiple lines in same shapefile using ArcMap?

texts:
Using ArcGIS 10.2

I require a method of editing multiple lines (selection) that share a common end point and to move the selection to an alternative common end point.

In the picture below I would like to move the light blue selection to a green point at the bottom of the picture.

Hopefully there is a tool for this that I am unable to find. That can quickly answer this question.

The lines are all in the same .shp file.


texts: Getting Response "XML request isn't valid" when inserting a feature using wfs-t

texts:
I am using Tinyows to perform wfs transactions.

Here is my snippet of tinyows configuration file.

<layer retrievable="1"
   writable="1"
   ns_prefix="tows"
   ns_uri="http://www.vizexperts.com/"
   name="target_entity"
   title="Target Entities"/>

and this is my xml template that I'm posting to insert a target.

    <wfs:Transaction version="1.0.0" service="WFS"         xmlns:wfs="http://www.opengis.net/wfs" xmlns:gml="http://www.opengis.net/gml"
xmlns:ogc="http://www.opengis.net/ogc">
  <wfs:Insert>
    <tows:target_entity xmlns:tows="http://www.vizexperts.com/">
      <tows:name>smallOne</tows:name>
      <tows:description>adding target</tows:description>
      <tows:category>whatever</tows:category>
      <tows:grid_no>15</tows:grid_no>
      <tows:latitude>67.777</tows:latitude>
      <tows:longitude>66.55</tows:longitude>
      <tows:elevation>100.0</tows:elevation>
      <tows:geom>
      <gml:Point xmlns:gml="http://www.opengis.net/xml" srsName="EPSG:4326"><gml:coordinates decimal="." cs="," ts=" ">66.55,67.777</gml:coordinates
          ></gml:Point></tows:geom>
      <tows:img_source>Some Image</tows:img_source>
    </tows:target_entity>
  </wfs:Insert>
</wfs:Transaction>

So after making the post request I'm getting the response that the xml is not valid.

Here is the exact response from server.

<?xml version='1.0' encoding='UTF-8'?>
<ows:ExceptionReport
 xmlns='http://www.opengis.net/ows'
 xmlns:ows='http://www.opengis.net/ows'
 xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance'
 xsi:schemaLocation='http://www.opengis.net/ows http://schemas.opengis.net/ows/1.0.0/owsExceptionReport.xsd'
 version='1.1.0' language='en'>
    <ows:Exception exceptionCode='InvalidParameterValue' locator='request'>
        <ows:ExceptionText>XML request isn't valid</ows:ExceptionText>

    </ows:Exception>
</ows:ExceptionReport>

Please can anyone help me to figure out what is wrong with the xml template ?

Here is the result of DescribeFeaureType.


texts: I'm building a windows 8 js application with an Iframe that uses openlayers with google as the basemap. Everything is working except for the map navigation. When I try to pan or zoom with my finger, the whole page moves, not the map. How do I fix this?


texts: What controls the coordinate value in the singleclick MapBrowserEvent argument?

texts:
I am trying to understand what is affecting the "singleclick" callback argument: MapBrowserEvent in my OpenLayers map (v6.4.3).
While I always receive the pixel coordinates back, I don't always get a coordinate back on the MapBrowserEvent object.
Here's one where the coordinate was returned:

And one where I do not get a coordinate back:

Once I receive an event with empty an coordinate, all subsequent events have an empty coordinate.  But no errors in the console.
I am wiring up events like this:
this.map.on("singleclick", function (evt) {
  console.log(`MapComponent.setupMap - Map - click - evt=`, evt);
  component.onMapClick(evt);
});

this.map.on("dblclick", function (evt) {
  console.log(`MapComponent.setupMap - Map - dblclick - evt=`, evt);
  component.onMapDoubleClick(evt);
});

I want to understand this behavior, before implementing a workaround (which I assume is just going to be converting the pixel coords to lat/lon).


texts: Given a polygon (no multipolygons, or holes) what is the best way to explode it into lines using T-SQL. The end results is I need a record for each line that makes up the polygon.


texts: Clip raster data by multiple shapefile in ArcPy using Clip data management tool

texts:
I have one raster data of a state and also I have 14 shapefile in same area. Now I want to clip 14 shapefile in arcpy using for loop. The issue is when I run the script its showing TypeError:'NoneType' object is not iterable. Here is the code
import arcpy

rasterlayer = (r"L:\\Arcpy\\Data\\Mosaic")
output = (r"L:\\Arcpy\\Data\\clip")
vector_layer = arcpy.ListFeatureClasses(r"L:\\Arcpy\\Data\\vector_data\\selected_district")

for vectordata in vector_layer:
  arcpy.CalculateValue_management("\"%Name%\",.split(\"_\")[0]", "","Variant")
  arcpy.Clip_management(rasterlayer, "", output, vectordata, "0", "ClippingGeometry", "NO_MAINTAIN_EXTENT")


texts: Creating new field considering duplicate and numeric values

texts:
I'd like to create a new field considering duplicate and numeric values.

For example, I have an attribute table with duplicate values in the field1 and numeric values in the field2. I'd like to create a field3 in which there are the suffix _1, _2...on the field1 attributes, when they are duplicates and when "field2">0.

In case "field2"=0 or not duplicated value in the "field1" -> "field3"="field1".

The suffix should not be set randomly, but following the spatial direction (or left or right); and for each group of duplicated values.

Before:

After:

Any suggestion to use a string in the field calculator to do this in QGIS?


texts: Error when using spatial query in arcgis silverlight api

texts: I'm developing a webgis using ArcGIS Silverlight API 3.1 . I used esri interactive SDK sample  for spatial query in the webgis .The problem is when i draw a freehand polygon and the polygon intersect itself (fig1), the error "Query failed: Error code '400':unable to complete operation." raised (fig2) . How can i solve the problem ?


texts:
I followed this methode https://geomappando.com/2016/11/08/come-realizzare-una-heatmap-con-leafletjs/  to get a heatmap from my GeoJSON file. But, I got it without having the value field included.

Here's my Geojson:

var lyrmesures = {
"type": "FeatureCollection",
"name": "mesures_test",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "id": 1, "_uid_": 1, "id_mesure": 17, "cm": 11627, "gisement": "EL HALASSA", "niveau": "Couche 3", "element1": 5.42, "element2": null, "element3": 6.558, "element4": null, "element5": 5.427, "element6": 5, "puissance": 1.9, "longitude": null, "latitude": null }, "geometry": { "type": "Point", "coordinates": [ -6.888298798254256, 32.639886473063939 ] } } ]}

Here's my code :

            geoJson2heat = function(geojson) {
            return geojson.features.map(function(feature) {
            return {lat:parseFloat(feature.geometry.coordinates[1]),
                    lng:parseFloat(feature.geometry.coordinates[0]),
                    value:parseFloat(feature.properties.element6)};
            });
            }

            var mesures = geoJson2heat(lyrmesures, 1);

            var myconf= {
                radius: 20,
                maxOpacity: 0.6,
                scaleRadius: false,
            };

            var heatmapLayer = new HeatmapOverlay(myconf);
            heatmapLayer.setData({max:1, data:mesures});
            mymap.addLayer(heatmapLayer);

Here bellow my result :


texts: How to Handle Negative Coordinates When Converting Long/Lat to Pixel Values

texts:
I'm writing a short script that draws white lines on a PNG file to represent roads. I'm using the coordinates from a .geoJSON file to draw the individual line segments. Unfortunately, because some of the coordinate values are negative, my image is stretched and thus becomes inaccurate compared to the actual data. How should I deal with these values?
Here is the conversion equation I'm using:
x = image_width * ((long-min_long)/(max_long-min_long))
"long" is the longitude coordinate in the geoJSON file
min_long and max_long are the mins and maxes of the input coordinates
Command:
python draw_roads.py C:/Users/'user'/Downloads/roads.geojson --world_bounds -74.016536 40.703922 -74.00874 40.701813 --image_size 2000 2000
Inputs:
FilePath: Location of .geoJSON file
--world_bounds: Real-world coordinates of the bounds of geoJSON file
--image_size: Size of output image


texts:
I'd like to identify mountain-peaks from an SRTM TIFF (created/composed from the original HGT ones).

How I'd wish it would work: Check for each pixel whether it is the highest value within a given radius (and I'd like to play with this radius). I know... it will generate some noise, but I will deal with it by having a larger "moving-frame" (i.e. bigger radius) up to 1-2 km.

This entry: how-to-calculate-peaks-in-a-dem asks basically the same and especially the (second) comment below the question from user/whuber phrases the problem perfectly - but I'd need a description (somewhat step-by-step) how to execute it.

The first (and chosen) answer would have been a nice start, but executing r.param.scale only throws an error ("Datenquelle ungültig" - "data-source invalid"), even though the same layer/file perfectly works with other GRASS tools.

Edit "Solution"

Turns out r.param.scale works perfectly, but you need to work in a projected coordinate system. So simply reprojecting to such a projected coordinate system solved my troubles.

within the r.param.scale one has to choose  for the 'Morphometric parameter' and this will categorize all pixels from 1 to 6, being 1 = planes, 2 = pits, 3 = channels, 4 = passes, 5 = ridges, 6 = peaks.

As expected the results vary vastly depending on the chosen parameters.

/Edit

I am aware of this entry too: how-to-extract-mountain-peak-points-from-a-given-digital-elevation-model But there is not much practical advice given.

It seems this is a somewhat popular question brought up sporadically by amateurs like myself...


texts:
I have a road layer, with two locations points. I have created a shortest path analysis between these points.

Now, I would like to delete some segment from the road network. I would like to delete (some) segments that the shortest path (displayed in bright red) lies on. Once this lines are deleted, I want to run the shortest path again.I want to do this in order to have two "shortest" paths that I can compare later on.
I have tried to create a multistring, then edit, select segment and delete. However, the segments remain in the layer.
How can delete segments, so my second shortest path analysis run through a different route?


texts: Calculate distance in QGIS, wrong geographic projection?

texts:
I'm trying to do the following: I have 2 shapefiles, one containing all the public schools in NY, another one containing the hotspots in NY (both in points form).
Now what I'm trying to do is calculate the distance from the schools to their closest hotspot. So that, when I merge it to my "SchoolPoints" layer and import it in TileMill I can visualise for each school the distance from the closest hotspot (for example: coloring the schools close to hotspots green, the ones that are far from a hotspot red).

To do this I used the distance matrix tool, with following settings:

(Note: pid is just a unique column I added, containing just integers)

After that I linked the data to the "SchoolPoints" layer (with as mutual field InputID and pid).

That results in a column, with a distance to the closest hotspot for each school, as it should be.

Well, not exactly. The problem I'm facing is: when importing my "SchoolPoints" layer in TileMill and after writing a small piece of code to plot the color differently according to the distance. For example:

[AfstandA7_ < 141.030340279] {marker-fill:#02D722;}
[AfstandA7_ > 1250] {marker-fill:#008A14;}
[AfstandA7_ > 2500] {marker-fill:#D70537;}
[AfstandA7_ > 3750] {marker-fill:#FF1F55;}
[AfstandA7_ >   5096.78247742] {marker-fill:#8A0123;} )

It shows wrong results. Points close to a hotspot that should be green show up red, and points far from hotspots show up green. Some points appear to show a correct color, so the problem is not in my CSS code.

A friend suggested that it might be related to the projection system I'm using. I already used two, both of them giving wrong results:

EPSG:31300    Belge 1972 / Belge Lambert 72
EPSG:2155     American Samoa 1962 / American Samoa Lambert (deprecated)

Now my question is:

What projection system should I use for NY? And also: might that solve my wrong distance values issue?

Here's a link with the shapefiles I used if you'd like to check it out yourself.


texts:
i want to start use Geoserver REST api, but have no idea how to use it.
In Geoserver docs i found some examples but then use cURL. Google says what cURL its console comand utilite. So i dont understand how to use it. For example if i want to make user autification using ExtJs or another javascript framework how can i use cURL. Or i gonna do somthing else?


texts: Which open source software to create desktop interface for Geo-processing tools?

texts: I want to create an desktop interface to create geo-processing tool and that to using Open source GIS software. I will be using Postgresql and postgis too to save and retrieve data for operations.  I  want to save resulted geo-processing data with naming conventions in Postgresql. For that i want to know Open sources Software to create interface, geo-processing tool and how tO save and retrieve same. Is it possible?


texts: Adjacency matrix from polygons in GeoJSON

texts:
I have a set of polygon features in a GeoJSON file. They are neigbourhoods in Chicago. I would like to produce an adjacency matrix A, in CSV format or any other form, where the i,j element of A is 1 if Pi is adjacent to Pj (i.e. they are touching).

Is there any tool that can easily read the file and do this?


texts: Refreshing QgsMapCanvas after removing selected layer using PyQGIS

texts:
I am working on a PyQGIS (QGIS 3.10) application to overlay data on OpenStreetMap raster tiles.  I am adding a function that removes the most recently added layer from the project and refreshes the map.  It looks like this:
def rmv_layer(self):

    if not self.layer_names:

        print("No overlay to remove.")

    else:

        name = self.layer_names[-1].strip('[]')
        pstr = "Removing overlay: %s" % name
        print(pstr)
        layer = self.project.instance().mapLayersByName(name)[0]
        self.project.instance().removeMapLayer(layer)
        self.QMC.refreshAllLayers()
        self.layer_names.pop(-1)
        self.layers.pop(-1)

    return

where self.layer_names is a list that is updated with the layer name after one is added to the map canvas in a separate add_overlay class while the name of the removed layer and the layer itself is erased here using pop.  This code works fine the first time around.  However, when adding a new layer to the map canvas using the following code:
def add_layer(self):

   layer = QgsRasterLayer(self.new_raster,self.layer_title,"gdal")

   if ( layer.isValid() ):

        self.layer_names.extend(layer_title)
        self.project.instance().addMapLayer(layer)
        self.map_canvas.setExtent(layer.extent())
        self.layers.insert(0,layer)
        self.map_canvas.setLayers(self.layers)

I get this error:
TypeError: index 1 has type 'QgsRasterLayer' but 'QgsMapLayer' is expected

I'm not sure what the root issue is since QgsMapLayer is the base class for QgsRasterLayer, QgsVectorLayer, etc.  The self.layers list is set up to where it always retains the base OpenStreetMap tiles in the 'registry'.  This error only occurs when I remove and refresh.  If I keep adding layers it's happy, and I have to refresh the map canvas otherwise the layer will still be displayed (e.g., commenting out self.QMC.refreshAllLayers() inherently does nothing).  This error also occurs without the self.layers.pop(-1) line in the first block.
Is there a better way to clear the most recent layer than this?


texts: Shapefile of continents not displaying correctly in GEE

texts:
I have a shapefile containing polygons that displays correctly in QGIS:

I would like to upload this as a shapefile to GEE and use as a mask to clip the JRC water extents to remove  the ocean water pixels near all the coasts.

However, when I upload the shapefile and plot the resulting table, it looks like this:

It apparently covers the entire globe, with Antarctica being inverted (not covered). The original shapefile had polygons for each continent, so I tried to dissolve this into a single geometry using QGIS's Dissolve tool, but I got the same result. The polygons have holes and are quite complex in some places.

I can share the shapefile if anyone wants to take a crack at it, but does anyone have ideas? I found a few notes about the clockwise/counterclockwise orientation of the vertices, but nothing definitive.

UPDATE:
Following a commenter's suggestion, I used Shape Escape to create a fusion table, and the result was as expected:

However, the original question remains: why doesn't GEE ingest/display the shapefile properly?


texts: About 71% of the Earth's surface is covered with water and the remaining 29% is landmass*. If you could theoretically cut the planet in half, with the objective of producing a half with the most landmass possible (and another with the most water possible), where would the cut be? More specifically, where would the centers (in the surface) of both halves be located?


texts:
Take a look at this map entitled Sunrise around the world:

They have animated where the sun is shining around the world as a moving layer in their torque animation. How would I go about doing this?


texts: strange zoom/pan behavior in QGIS on layers extracted from USGS Geopdf topo

texts:
Starting with a USGS geopdf topo 7.5 quad, I extracted the layers to a spatialite database using ogr2ogr:

ogr2ogr -a_srs "WGS84" -f sqlite fisher.sqlite UT_Fisher_Towers_20170216_TM_geo.pdf

I then load these layers into a new map in qgis (2.18). Zooming to layer extent does not actually zoom to extent, but some seemingly random portion of the map. Same for zoom full. Zooming in with the scroll wheel sometimes zooms in, sometimes out. Panning with the pan (hand) tool can cause the zoom level to change. The scale shows 591,657,527 for multiple different zooms on the screen.

I've tried this on two different topos and have the same results. The behavior is the same whether I load one layer or all the layers in the map.

Not sure if it matters, but the ogr2ogr conversion/extract generates multiple errors, each of which reads:

ERROR 1: Max token size reached


texts: Get latitude and longitude of a geography column


texts: QGIS Expression Calculator Field Shows Null

texts: When I go into the expression builder under the "Field and Values" drop down on the left side, instead of a list of field for the layer all it says is NULL. Is there something that would cause this instead of a list of fields?


texts: PyQGIS error when running zonal statistics?

texts:
I want to extract the sum of my raster cell values with a specific area through zonal statistics and have prepared the following code (see below). When I run the code I get the following error message (see bottom).

Any suggestions?

from qgis.analysis import QgsZonalStatistics

#specify polygon shapefile vector
polygonLayer = QgsVectorLayer('C:/Users/Anna Nilsson/Desktop/AIRPOLIM-ES/GIS data/Kenya/cpp_KE_buf1.gpkg', 'zonepolygons', "ogr")

QgsProject.instance().addMapLayer(polygonLayer)

# specify raster filename
rasterFilePath = 'C:/Users/Anna Nilsson/Desktop/AIRPOLIM-ES/GIS data/Kenya/KEN_worldpop/KEN_popmap15adj_v2b.tif'

# usage - QgsZonalStatistics (QgsVectorLayer *polygonLayer, const QString &rasterFile, const QString &attributePrefix="", int rasterBand=1)
zoneStat = QgsZonalStatistics (polygonLayer, rasterFilePath, 'pre-', 1,QgsZonalStatistics.Sum)
zoneStat.calculateStatistics(None)

Error message:

Traceback (most recent call last):
  File "C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\code.py", line 90, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 12, in <module>
TypeError: QgsZonalStatistics(): arguments did not match any overloaded call:
  overload 1: argument 2 has unexpected type 'str'
  overload 2: argument 2 has unexpected type 'str'
  overload 3: argument 1 has unexpected type 'QgsVectorLayer'


texts: How can I efficiently use a list comprehension on an ogr.layer inside a loop in Python?

texts:
Problem statement

I'm using python 2.7 and ogr 1.11.0 on Windows 7.

I have two shapefiles: one for streets and one for junctions. The streets shapefile has fields listing the IDs of the start and end junctions, as well as additional information about the streets.

I need to loop over each junction, and collect information about all of the streets that are bounded by it. I started off with a spatial filter, like this:

for junction in junctions:
    streets.SetSpatialFilter(junction.GetGeometryRef())
    connected_streets = [street for street in streets]

But this was slow (about 1 second per junction, and I've got 191,000 junctions). So I tried to take advantage of the SRT_JUNC and END_JUNC fields in the roads shapefile, to see if this would work faster.

for junction in junctions:
    junctionID = junction.GetField("ID")
    connected_streets = [street for street in streets if street.GetField("SRT_JUNC") == junctionID or street.GetField("END_JUNC") == junctionID]

The list comprehension works okay in the first iteration of the loop, but for all of the following iterations it returns nothing (or a list with length zero?).

The only way I could get it to work is to Open the driver and GetLayer inside the loop. But I've got a large shapefile and the performance of this method is also too slow.

I show some simplified examples below.

Can I make multiple list comprehensions work for an ogr layer, without reloading the entire layer each time? If not, why not? (I would also accept suggestions for an entirely different approach to the main problem, but I still am looking for an explanation regarding list comprehension.)

Example 1: List comprehension without ogr layer

The following basic script

my_list = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
for x in range(3):
    print len([a for a in my_list])

produces the expected result:

10
10
10

Example 2: List comprehension with ogr layer

The following script is parallel to Example 1, except it uses an ogr layer instead of a simple list.

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(file, 0)
streets = dataSource.GetLayer()
for x in range(3):
    print len([street for street in streets])

It produces an unexpected result:

191017
0
0

Example 3: move GetLayer inside the loop

I moved the GetLayer call into the loop:

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(file, 0)
for x in range(3):
    streets = dataSource.GetLayer()
    print len([street for street in streets])

But still it produces the unexpected result:

191017
0
0

Example 4: move driver.Open inside the loop

I also moved the driver Open call into the loop:

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")
for x in range(3):
    dataSource = driver.Open(file, 0)
    streets = dataSource.GetLayer()
    print len([street for street in streets])

And it produces the expected result:

191017
191017
191017


texts: Change entire label style after certain extent

texts:
I need to design a dynamic label in ArcMap that changes its entire style after a certain zoom extent. The research I've done suggests to me that there is currently no way to use VBscript to manipulate certain label stylings. I haven't discovered a way to, for instance, programatically add a callout box and a label offset after I zoom out past 1:150. Below are two label styles I'd like to fuse together into one dynamic style.

Zoomed in past 1:150, and then zoomed out past 1:150:

Will this only be doable in Python?


texts: QgsAttributeDialog in PyQGIS Standalone Application

texts:
I am trying to implement a desktop application with PyQGIS 3, I have not been able to program a layer info using the attribute dialog that the QGIS API.
This is part of the code that I use:
def info(self):
    self.vlayer = QgsVectorLayer('provincias.shp', "Provincias", "ogr")
    QgsProject.instance().addMapLayer(self.vlayer)
    layer = self.vlayer
    features = layer.getFeatures()
    feature_list = list(f for f in features)
    feat = feature_list[0]
    attribute_dialog = QgsAttributeDialog(layer, feat, True)
    attribute_dialog.show()

shows the following error:
Municipio Failed to create widget with type 'TextEdit',
Provincia Failed to create widget with type 'TextEdit'

How to solve this problem?


texts: map does not zoom to a feature - openlayer's zoomToExtend

texts:
I'm using Geoserver 2.1.3 and openlayers 2.12. I created a vector layer and I implemented a "Search" functionality using CQL. After the CQL search, I succesfully render new features on the layer. And then I use the following code, so the map will zoom to the new features.

map.zoomToExtent(myLayer.getDataExtent());

If I use this code with the other vector layers, works just fine. But for this particular layer, does not zoom to the new features. Keeps zooming on a certain, irrelevant area.

Any tips/hints?

Thanks


texts: Get raster coordinates after warping with C# GDAL

texts:
I'm currently trying to read a raster file with GDAL and C#. The reading part is fine, but I have troubles with reading the correct bounding box of the raster to map it correctly to the screen. If the coordinates are in lat, long then it's fine. However, if they're in pixel values (for example with EPSG:3857 Web Mercator), I can't read them out correctly. I know Gdalinfo include the original bounding box (see below).

Driver: GTiff/GeoTIFF
Files: F:\ArcGis\Gdal\NE1_50M_SR_W\DensityMaps\Pop1Warpped.tif
Size is 141988, 60942
Coordinate System is:
PROJCS["World_Mollweide",
GEOGCS["GCS_WGS_1984",
    DATUM["D_WGS_1984",
        SPHEROID["WGS_1984",6378137.0,298.257223563]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.017453292519943295]],
PROJECTION["Mollweide"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
UNIT["Meter",1.0]]
Origin = (-17619594.547443531000000,8750529.461868489200000)
Pixel Size = (250.000000000000000,-250.000000000000000)
Metadata:
    AREA_OR_POINT=Area
Image Structure Metadata:
   COMPRESSION=LZW
   INTERLEAVE=BAND
Corner Coordinates:
   Upper Left  (-17619594.547, 8750529.462) (  4d35'30.71"W, 83d37'50.88"N)
   Lower Left  (-17619594.547,-6484970.538) (107d 3'58.73"E, 55d59' 2.97"S)
   Upper Right (17877405.453, 8750529.462) ( 15d11'38.88"E, 83d37'50.88"N)
   Lower Right (17877405.453,-6484970.538) (103d21'55.35"W, 55d59' 2.97"S)
   Center      (  128905.453, 1132779.462) (  1d17'47.23"E,  9d10'35.76"N)

As you can see the upper left includes the coordinates in lat, long as well. Is there anyway I can read it out programmatically?


texts: Accessing existing instance of QGIS in stand-alone Python script

texts:
For my own data analysis pipeline, I am running a split screen setup with IPython notebooks on one side, generating variables from source data, pushing them to a CSV file which QGIS then reads/joins to an appropriate geography, instantly visible on the other side of the screen.

I would like to push the control of the running instance of QGIS to the ipython notebook, in effect allowing a super quick calculation/visualization cycle.

Right now I can do this in the QGIS python console with something like this:

canvas = qgis.utils.iface.mapCanvas()
cLayer = canvas.currentLayer()
canvas.refresh()

but I would like to:

connect to an existing instance of qgis from an external Python
environment
at a minimum refresh the map, after each re-write to the csv from ipython notebook, but preferably
re-calculate breaks (using quantiles, natural breaks) based on the newly joined values

How can I achieve this?


texts: How to get the second nearest neighbor between two point patterns in R?

texts: Is there a way to get the distances for the second nearest neighbor between two point patterns in R? The spatstat package has a function called nncross but it only applies to the nearest neighbors between two patterns and I need the distances to the second nearest neighbors.


texts: How to display a large series of images to the map with a for loop in Earth Engine

texts:
I need to view imagery from a collection as separate layers. Since I have to do this for several sites and the number of images in the collection therefore varies, I tried to apply a for loop to solve this.

My working code without for loop looks like this:

// Load a landsat image and select three bands for Mount Spurr

// HSV-based Pan-Sharpening of Landsat 7 TOA images.
var panSharpenL8 = function(image) {
  var rgb = image.select('B3', 'B2', 'B1');
  var pan = image.select('B8');
  // Convert to HSV, swap in the pan band, and convert back to RGB.
  var huesat = rgb.rgbToHsv().select('hue', 'saturation');
  var upres = ee.Image.cat(huesat, pan).hsvToRgb();
  return image.addBands(upres);
};

var dataset = ee.ImageCollection('LANDSAT/LE07/C01/T1_TOA')
                  .filterDate('2004-04-01', '2006-06-30')
                  .sort('system:time_start')
                  .map(panSharpenL8)

// short form
var Image_filtered = dataset.filterBounds(geometry_Mount_Spurr)
.filterMetadata('CLOUD_COVER','less_than',10)
print(Image_filtered.size())

var listOfImages = Image_filtered.toList(Image_filtered.size());
print('List:',listOfImages);

var img1 = ee.Image(listOfImages.get(0));
var img2 = ee.Image(listOfImages.get(1));
var img3 = ee.Image(listOfImages.get(2));
var img4 = ee.Image(listOfImages.get(3));
var img5 = ee.Image(listOfImages.get(4));
var img6 = ee.Image(listOfImages.get(5));
var img7 = ee.Image(listOfImages.get(6));
var img8 = ee.Image(listOfImages.get(7));
var img9 = ee.Image(listOfImages.get(8));
var img10 = ee.Image(listOfImages.get(9));
var img11 = ee.Image(listOfImages.get(10));
var img12 = ee.Image(listOfImages.get(11));
var img13 = ee.Image(listOfImages.get(12));
var img14 = ee.Image(listOfImages.get(13));
var img15 = ee.Image(listOfImages.get(14));
var img16 = ee.Image(listOfImages.get(15));
var img17 = ee.Image(listOfImages.get(16));
var img18 = ee.Image(listOfImages.get(17));

// add pansharpened bands

Map.addLayer(img1, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img2, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img3, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img4, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img5, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img6, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img7, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img8, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img9, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img10, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img11, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img12, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img13, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img14, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img15, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img16, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img17, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
Map.addLayer(img18, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');

//Map.addLayer(Image_filtered, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');

As you can see, there are many Map.addLayer expressions. I wish to not copy and paste the lines to add each image to the Map, so I built an ee.List and looped like this:

var listOfImages = ee.List(Image_filtered);
print('List:',listOfImages);

var len = listOfImages.length;

for (var i=0; i < len; i++) {
  var img = ee.Image(listOfImages.get(i));
  Map.addLayer(img, {bands:['red','green','blue'], max: 0.3}, 'Pansharpened');
}

However, it does not add layers to the map and does not result in an error message.

How do I solve the problem?


texts: WFS layers disappear when zooming out, does not appear back when zooming in

texts:
I work with a wfs layer. Because they are lots of features in this layer and wfs visualization can not handle them all, I unchecked the "cache" box.

It works all fine except that some entities disappear when zooming out (scale-dependent rendering), and do not come back when zooming in again. Seems like the wfs protocol is looking for new entities only at zoom out, and that the ones that disappear due to rendering issues are not coming back.

Does anyone as an idea on whether this analysis is correct, and how to overcome the problem ?


texts: Checking if point is inside irregular polygon?

texts:
I'm using Javascript (https://github.com/substack/point-in-polygon) to test if a point is inside polygon, however I think it's adapted from https://www.ecse.rpi.edu/~wrf/Research/Short_Notes/pnpoly.html

Following are the list of points,

let points = [
            [ 42.34624, -71.06024 ],
            [ 42.33558, -71.06616 ],
            [ 42.32632, -71.05835 ],
            [ 42.32987, -71.05428 ],
            [ 42.34732, -71.05432 ],
            [ 42.34618, -71.05973 ],
            [ 42.34624, -71.06024 ]
        ];
let testPoint = [
            [42.33288, -71.05804]

        ];`

let isInside = pp(testPoint, points);
expect(isInside).to.equal(true); // always fails

I'm not sure why it fails though. Is it because the algorithm doesn't cope with this "type" of polygon? Any hints please? Please point me to any  alternatives to this approach as well.


texts: Adding duplicate labels within a polygon - QGIS

texts:
Is there a way to include multiple duplicates of the same label within a single polygon in QGIS?
So a single polygon(s) with a label from a single field, repeated at a set distance within the polygon - shown in the image as made by Arcmap Maplex engine.


texts: How can I "reset" SketchUp origin after importing DXF with "Preserve Drawing Origin" option enabled?

texts:
I can successfully import all of my 3D DXFs into a common location so that they are referenced correctly (in reference to each other). TO do this I simply check the Import DXF "Preserve Drawing Origin" option on.

However, when I do this DFX models are placed very far away from the Sketchup "default" origin. It becomes very difficult to zoom to the location of the 3D DXFs and navigate around them.  Zoom Extents zooms out all the way showing the default origin and the imported DXFs as a spec in the distance (see screen capture). To even zoom in to the imported model is nearly impossible.

How do I "reset" the origin and correct the extent so that sketchup zooms in onto the model when "Zoom Extent" is clicked?


texts: Is the Azimuth on the equator equal to the one not on the equator?

texts:
I am trying to fully understand the concept of Azimuth and I am encountering some inconsistencies (or maybe it is my error).

I show you some examples that do not match, hoping somebody can explain me how this really works.

I show the coordinates in EPSG:900913, in PostGIS and using my own JavaScript function.

MY FUNCTION

/* Difference between the two longitudes */
var dLon = lon2 - lon1;
/* Y value */
var y = Math.sin(dLon) * Math.cos(lat2);
/* X value */
var x = Math.cos(lat1) * Math.sin(lat2) - Math.sin(lat1) * Math.cos(lat2) * Math.cos(dLon);
/* Calculates the azimuth between the two points and converts it to degrees */
var angle = Math.atan2(y, x) / Math.PI * 180;

EXAMPLES

/* Same Y, not on the equator */
Point A: (-81328.998084106, 7474929.8690234)
Point B: (4125765.0381464, 7474929.8690234)
Result in PostGIS: 90 degrees
Result in my JS function: 74.232 degrees

/* Same Y, on the equator */
Point A: (-81328.998084106, 0)
Point B: (4125765.0381464, 0)
Result in PostGIS: 90 degrees
Result in my JS function: 90 degrees

I understand, that on the equator, the Azimuth is 90 (or 270) for a horizontal line.  think that if you draw a horizontal line a bit North (or South) of the equator, then the Azimuth is not 90 degrees anymore. But... PostGIS tells me that it is always 90 degrees when we have the same Y.

Additionally, this calculator also show that the Azimuths of horizontal lines are not 90 degrees when Y != 0 (not on the equator).

How is it correct?

Thanks


texts: Add diferents animated widgets with diferents layers at the same map

texts:
I'm making a map with 4 layers that we want to represent them as animated maps.
I think Carto doesn't let to have two animated widgets in the same map.
Does anybody know if it is possible to make two or more animated widgets?


texts: Creating ArcGIS report with multiple related tables?

texts:
Is there a way to create a report in ArcGIS 10.3 that includes data from multiple related tables in a file geodatabase?

I'll give an example using soil profile data. I have a point feature class representing soil cores. Each point has data specific to that site that the core was extracted from. This table has a common ID and a relationship class, with a one to many relation, with another table that has a separate record for each layer in the soil core as it gets further down the core. Is there a way to select one soil site point, make a report for that point with the ArcGIS report wizard, and include the soil horizon records for that point in the report?


texts: Subtract each raster in a list by one before it

texts:
I have a list of rasters named like so:

19940506.tif
19950506.tif
19960506.tif
19970506.tif

and I want to subtract 19950506 by 19940506 and 19960506 by 19950506 etc.

I can do this manually with something like this:

arcpy.env.workspace=r'F:\pathway_to_rasters'
May1List=arcpy.ListRasters('*.dat')
outpath=r'F:\outpath'
outRaster = Raster(May1List[1]) - Raster(May1List[0])
out='19950506_19940506'
outpath=os.path.join(outpath,out + '.tif')
outRaster.save(outpath)

but if would like to try and find a way to automate this in a loop if possible.


texts: Layer disappears when zooming out past certain point in ArcGIS Online

texts:
I have shared several layers from ArcGIS Pro to ArcGIS Online, and one of the layers disappears when I zoom out past a certain point in the map (1:320,000). All of the layers work as they should in ArcGIS Pro. The other layers (with the similar settings) resize when I zoom in and out, and all of the layers have the same visibility range. The layer I have problems with changes symbology from a thin line (1:160,000-1:80,000), a strange lump-shape (1:40,000-5,000), towards correct shape and extent at 1:2,500.
I have checked the Spatial Index and Check Geometry.
What else can I do?


texts: Using Custom Legend Symbol in ArcGIS Online?

texts:
I am currently creating an ArcGIS Online Journal story map. In one of the web maps, I have a point feature class, where the points are transparent and they provide labels for regions (with no clear boundaries, hence no polygons). I need the legend to show one of the labels instead. This way, the legend symbol would be an example of the label text.

Is there a way to do this?


texts:
I'm trying to write a python script that adds to a csv file the date, time and mxd path.

I tried to run the script from the Python IDLE without the mxd path and I succeeded to add the date and the time but when I trying to run the script from a tool box in ArcMap it does not work.

Below is my code so far:

import csv
import os
import time
import arcpy
date=time.strftime("%d/%m/%Y")
time1=time.strftime("%H:%M:%S")
mxd=arcpy.mapping.MapDocument("CURRENT")
path=str(mxd.filePath)
with open('log.csv', 'a') as fp:
    a = csv.writer(fp, delimiter=',')
    data = [[date, time1,path]]
    a.writerows(data)

It work now but i still have two problems:

My mxd path have Hebrew letter so i received - "'ascii' codec can't encode characters" i tried to use .encode('utf-8') but it didn't work.
In the csv file i getting a row blank between row to row, How can i fix it?


texts: How do I cut a donut hole in a raster in QGIS

texts:
I have a DEM(1) of a bigger area and a higher resolution DEM(2) in a small part of the same area. I would like to cut a hole in DEM1 and replace it with DEM2.

Is there a non python way of cutting a hole from a mask in a raster layer in QGIS?


texts:
We are trying to get single wms layer from Geoserver. We know how to pick out the wms but we want to be able to pick out wms for layergroup that we have created within our geoserver.

http://docs.geoserver.org/stable/en/user/services/virtual-services.html

Here we can pick out wms from a layer or a workspace but not from layergroups.

Is it possible to create wms from layergroup in Geoserver?

Sjonni, Iceland


texts: how to create thiessen polygon around specific areas of interested area?

texts:
I wish to create Thiessen polygon within all sub_basin of my interested area individually, meaning that to create Thiessen for every sub_basin one by one. I have to layer one is point layer and another is the shapefile of the study area which my sub_basins is located within it. I know it can be done using manually and also using ModelBuilder.

I prefer to use ModelBuilder.


texts: How to conveniently generate a Google Maps link from a PostGIS geometry?

texts:
When working with a spatial database, it is often convenient to make a sanity check for a location of a point by viewing it on Google Maps (or any other online mapping service).

Copying the lat/lon is tedious - it requires me to manually copy-paste data from the console to the browser, and since st_asText order (lon, lat) is the opposite of Google (lat, lon) I have to do it twice and add a comma in between. Impractical for fast checking a dozen or so points.

How do I generate a valid Google Maps or a Bing map link from a PostGIS geometry?


texts: Comparing between layers name and merging layers having same name using PyQGIS?

texts:
I'm trying to write python script in QGIS to merge layers having the same name.

To do that I specify the name of layer in the script, so all layers having this name, will be merged.  However, since  I have a lot of names this is not practical. I tried to convert the names type from unicode to list to be able to iterate in the list of names.

In this script, I specify layer name. I want to iterate over all names and compare them to finally merge layers having same name.

import qgis.core
import processing
def __init__(self, iface):
    self.iface = iface
lyr_list=[]
count=1
layers=iface.mapCanvas().layers()
for i in layers:
    if i.name() == 'layer_name':
        lyr_list.append(i)

processing.runandload("qgis:mergevectorlayers", lyr_list, "memory:merged")


texts: Defining the Laskowski Tri-Optimal in QGIS

texts:
I'm fairly new to the world of GIS software and I'm trying to project the Laskowski Tri-Optimal projection in QGIS. It has a proj.4 definition (+proj=lask), but no parameters are listed for it on spatialreference. I did find one website that used the string

+proj=lask +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m

but even this gave me the error "the proj4 definition is not valid" when entered into QGIS's CRS menu. I've successfully entered other projections, so I'm fairly certain my issue isn't a procedural one with the definition menu. Is it possible that QGIS or even GDAL don't support the Laskowski?


texts: How to show a correct scale for mercator-based maps in print composer?

texts:
To be more specific:

I have created a map in Pseudo-Mercator EPSG: 3857 using Google Satellite.
So my layers, and my project, are all set to Pseudo-Mercator.

How can I create a map in print composer which shows me a right scale?

I mean, when I create a map of my area, I know that the scale is wrong when compared to official maps of that area (my map shows a smaller area of the world).

I cannot indicate the scale (e.g. 1/10000) because it is wrong.

What can I do? How can i create maps in typical scales like 1:10000 or 1:100000 which are correct?

I am using QGIS 2.18.2.


texts:
I am having some trouble with GDAL (v2.1-dev).  I build it from source today using Poppler and OpenJPEG.  All the command line utilities work fine, but none of the Python plugins work at all, producing this error message:

Traceback (most recent call last):
  File "/usr/local/bin/gdal_merge.py", line 37, in <module>
import gdal
  File "/usr/local/lib/python2.7/dist-packages/gdal.py", line 2, in <module>
from osgeo.gdal import deprecation_warn
  File "/usr/local/lib/python2.7/dist-packages/osgeo/__init__.py", line 21, in <module>
_gdal = swig_import_helper()
  File "/usr/local/lib/python2.7/dist-packages/osgeo/__init__.py", line 17, in swig_import_helper
_mod = imp.load_module('_gdal', fp, pathname, description)
ImportError: libgdal.so.20: cannot open shared object file: No such file or directory

Now, I've seen plenty of threads that say update LD_LIBRARY_PATH or include /usr/local/bin/ in my ld.so.conf file. I've done both of those things and they don't help.  Furthermore, those posts all refer to libgdal.so.1 not being found. So, I added the location of libgdal.so.20 to ld.so.conf and still no dice.

Any recommendations?

Update (11/2/15)
I am running Ubuntu 14.04.03 32bit, cleanly installed just a few days ago.
I followed these instructions to download, make, and install GDAL and OpenJPEG in a folder I made called /home/ronn4031/src/.


texts: How to Optimize Openstreetmap Tile Rendering?

texts:
I am currently trying to render Water (3,198,333 records) and detailed world borders (processed_p) for the world. I am successful in rendering both layers, but I am not happy with performance.

My current setup is as follows:

Amazon EC2 instance (m2.2xlarge)
34.2 GB of memory
13 EC2 Compute Units (4 virtual cores with 3.25 EC2 Compute Units each)
850 GB of instance storage
64-bit platform
I/O Performance: High
Mapnik, 4 threads
Postgresql settings based on these benchmarks
GiST index in place for all of my layers

Some of my Ideas to Increase Performance:

Simplifying my features
Don't render water tiles - when loading with openlayers can set load failure to blue color
Split world water features into various sections and explicitly tell mapnik its bounds. For example, if I split water features from North America, when Mapnik is generating tiles for Australia it should not query North America to see if any of its features are in the tiles.
My second idea is a bit more crazy. Since Mapnik performs a spatial query for each tile in order to find features to process, I can perform this query beforehand for all of my layers. I would assume that this would save some processing power.

My Question:
What else can I do to optimize tile rendering?


texts: Creating average polyline from set of polylines

texts:
I am trying to plot the average of several polylines (up to 40). I have several shorelines of a short section of beach plotted from different dates and would like to plot their average position.

Is there a tool that will allow me to do this?


texts: Remove Holes Existing Within a Single Polygon in ArcGIS 10

texts:
I have a shapefile that consists of a single polygon.  I have gaps within the polygon that need to be eliminated.  The resultant dataset shoulde be a single, solid polygon.

I have tried implementing topology rules (rule: must not have gaps) on the dataset (it's in a geodatabase and within a feature dataset), but I don't think this is the appropriate solution/technique to accomplish this task. Any ideas?

I'm on ArGIS 10 with an ArcInfo License.


texts:
I'm trying to teach myself to make my first QGIS plugin. I used the Plugin Builder to create the base python files, and the plugin seems to work fine when I press the toolbar icon. I'm now trying to allow the user to use a keyboard shortcut to trigger the plugin as if they were pressing the "ok" button in the plugin dialog - assuming the plugin dialog is already open.  I've done something wrong and now my menu item is repeated in more than one place and I can't remove them:

I had trouble getting the method described here to work, and only later realized I wasn't unloading correctly, or at all.

My plugin loading in the initGui is certainly wrong, but so is my
current unload wrong:

def unload(self):
    """Removes the plugin menu item and icon from QGIS GUI."""
    for action in self.actions:
        self.iface.removePluginMenu(
            self.tr(u'&Copy Attributes'),
            action)
        self.iface.removePluginMenu(
            self.tr(u'&Please Copy the Attributes'),
            action)
        self.iface.removeToolBarIcon(action)
    # remove the toolbar
    del self.toolbar
    self.iface.unregisterMainWindowAction(self.keyAction)

I've used the Settings -> Configure Shortcuts manager to remove shortcuts, but I have multiple entries there as well. I originally tried to create a shortcut for the plugin by just using the Shortcut Manager, but it didn't work, so that's why I started trying to do it in the initGui.

I'm not sure of the correct method to unload a menu or toolbar item, but will ask this in a different question.

Update: I've managed to figure out how to remove the duplicate menu items. I simply restarted QGIS without the plugin loaded. Then reloaded it from the Plugin Manager. The duplicate menu entries have now disappeared.


texts: pgRouting pgr_tsp() not returning results when used with SQL parameter


texts: Calculating elevation profile for a polyline in ArcGIS 10?

texts:
My elevation data is a floating point tif. I have a road for which I need to calculate the elevation profile, but I can't see the tools to do this.

Can I do this with ArcGIS, with spatial analyst/3D analyst?


texts: Convert information in a text file to a map using ArcGIS Runtime SDK?

texts:
I am new to the ArcGIS Runtime SDK. What we are trying to do is to convert information in a text file (contains country name and some information) to a map (visualize the information as a heat map for example). Do we have to convert the text file to a shapefile?

If yes, does the Runtime SDK provide API's for it or we have to go back to Engine SDK?


texts: ArcPy Clip Raster ERROR 000964: Specified extent is invalid

texts:
I am trying to write an ArcGIS script, however, I can't seem to get Clip Raster to cooperate. This is what I have so far:

import arcpy, math

if __name__ == '__main__':
    inDEM = arcpy.GetParameterAsText(0)
    points = arcpy.GetParameterAsText(1)
    bufferSize = arcpy.GetParameterAsText(2)

    arcpy.AddMessage("Creating buffer...\n")
    arcpy.Buffer_analysis(points, "bufferOutput", bufferSize + " meters", dissolve_option = "ALL")

    arcpy.AddMessage("Projecting buffer...\n")
    arcpy.Project_management("bufferOutput", "bufferProject", arcpy.Describe(inDEM).SpatialReference)

    arcpy.AddMessage("Clipping DEM using buffer...\n")
    arcpy.Clip_management(inDEM, "bufferProject", "clippedRaster", clipping_geometry = "ClippingGeometry")

Every single time it throws the following error when it gets to the Clip Raster part.

Traceback (most recent call last):
  File "C:\Users\User\OneDrive\Research\Scripts\script.py", line 22, in <module>
    arcpy.Clip_management(inDEM, "bufferProject", "clippedRaster", clipping_geometry = "ClippingGeometry")
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 16854, in Clip
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\management.py", line 16851, in Clip
    retval = convertArcObjectToPythonObject(gp.Clip_management(*gp_fixargs((in_raster, rectangle, out_raster, in_template_dataset, nodata_value, clipping_geometry, maintain_clipping_extent), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError:  ERROR 000964: Specified extent is invalid
Failed to execute (Clip).
 Failed to execute (script).

Now the first thing I checked was the coordinate systems, but they are the same. Heck, even in my code you can see I set the coordinate system to match. And when I run this exact same setup manually in ArcGIS it works without a hitch.


texts:
I have a map where a user can add features or arrows to certain areas to "mark" the map.  However, upon zooming in/out of the map, the arrow coordinates are centered for the center of the arrow, not the arrow tip.  How could this be fixed?

I think we are using this to set up the layers:

_this.featureOverlay = new ol.layer.Vector({
    source: new ol.source.Vector({
        features: _this.features})
})

I should also clarify the arrows are font awesome arrows, so they are text.

They are styled like so:

style = new ol.style.Style({
    text: new ol.style.Text({
        text: _this.iconUnicodes[shape],
        font: 'normal ' + fontSize + ' FontAwesome',
        textBaseline: 'middle',
        fill: new ol.style.Fill({
            color: color,
        }),
        rotation: rotation
    })
 });

Added new code.


texts: Copy Layers with groups and layers with data from one MXD to another

texts:
Is there a way of copying one group of layers with data (001_Group) from one MXD(group1.mxd) to another(group2.MXD) with arcpy. As far as I understand all the layers in ArcMap are referenced as standalone layers in arcpy, so arcpy wouldn't know if for example 010_Subgroup is beneath 001_Group.

I have hundreds of large groups which have the same structure as 001_Group but layer_with_data layers have different data sources. I already have a script to update all the layers within 001_Group (in group1.mxd) with 002_Group datasources and save as group2.MXD.

However, I want to keep all the groups within the same MXD, so 001_Group, 002_Group and so on are listed from top to bottom. There was one question posted how to move all layers with looping through with arcpy.mapping.ListLayers, but this does mess things up and layer structure is inverted or not in correct order, as you can see from the example TOC below that some layers need to be left empty or lots of subgroups.

Here is the code which produces Group_002.mxd, Group_003.mxd and Group_004.mxd separately, which I'm trying to make into a single MXD with Groups being organised from top to bottom:

import arcpy
########################################################################################################
mxd = arcpy.mapping.MapDocument(r"C:\temp\Proj\Group_001.mxd")
subGroupPath_010 = r"C:\temp\Proj\subGr_010"
subGroupPath_020 = r"C:\temp\Proj\subGr_020"
subSubGroupPath_021 = r"C:\temp\Proj\subGr_020\ssubGr_021"
subSubGroupPath_022 = r"C:\temp\Proj\subGr_020\ssubGr_022"

grp_noArr = range(2,5)
grp_nameArr = ["Group_002", "Group_003", "Group_004"]
########################################################################################################
for grp_no, grp_name in enumerate(grp_nameArr, grp_noArr[0]):
    fileName111 = str(grp_no) + "_" + grp_name
    srcOutName50d = str(grp_no) + "_" + grp_name

    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == "001_Group":
            lyr.name = grp_name
        if lyr.supports("DATASOURCE"):
            for lay_i in range(1,5):
                if lyr.name == "10" + str(lay_i) + "_layer_with_data":
                    lyr.replaceDataSource(subGroupPath_010, "RASTER_WORKSPACE", fileName111 + "_subgrp010_" + str(lay_i))

                elif lyr.name == "21" + str(lay_i) + "_layer_with_data":
                    lyr.replaceDataSource(subGroupPath_020, "RASTER_WORKSPACE", fileName111 + "_subgrp020_" + str(lay_i))
                elif lyr.name == "22" + str(lay_i) + "_layer_with_data":
                    lyr.replaceDataSource(subGroupPath_020, "RASTER_WORKSPACE", fileName111 + "_subgrp020_" + str(lay_i))
                elif lyr.name == "23" + str(lay_i) + "_layer_with_data":
                    lyr.replaceDataSource(subGroupPath_020, "RASTER_WORKSPACE", fileName111 + "_subgrp020_" + str(lay_i))

    mxd.saveACopy(r"C:\\temp\Proj\\" + grp_name + ".mxd")


texts: map.pause() stopping LocationService with 10.2 SDK

texts:
I have an app in which I want the LocationService to remain running even if onPause() is called. I have a service running in the foreground with a handler & broadcast receiver that is bound to an activity. The service/activity work together to updates a trackline at an X time interval, both on the map & in a SQLLITE database on the device.

When I call map.pause() in the onPause() method, the locationservice stops listening for location change until I unpause the map again in onResume().

I initialize the LocationService using map.getlocationservice(), so it makes sense that it would stop the LocationService, but I don't want this to happen.

Is there a work around for this? Or perhaps something I'm not understanding correctly?

cross post: http://forums.arcgis.com/threads/99435-map.pause()-stopping-LocationService-with-10.2-SDK


texts: Perhaps a basic question, but can php code be used in editing the CartoDB infowindow HTML? We want to show or hide a category depending on the value.


texts:
In the past I was able to create a layer and specify interactivity and which fields would be returned on featuerOver etc. The code below worked in the past but now returns and error of "Interactivity is read-only in NamedMaps"

cartodb.createLayer(map, 'link to my viz')
.addTo(map)
.on('done', function(layer) {
  layer.setInteraction(true);
  layer.getSubLayer(0).setInteractivity("cartodb_id,_name")
  layer
    .on('featureOver', function(e, latlng, pos, data) {
      console.log(data);
    })
}).on('error', function(err) {
  console.log("some error occurred: " + err);
});


texts: What object argument to use in subclass definitionIn if subclassing ArcPy Spatial Reference class?

texts:
I want to do something like the following:

class SubClassOfSpatialReference(arcpy_spatial_reference_object):
    def __init__(self, input_data):

... etc.

Can I even do this and, if so, what is the parent object that I need to provide in the class definition?

If you type() a spatial reference object, you get arcpy.arcobjects.arcobjects.SpatialReference.  I tried that and got AttributeError: 'module' object has no attribute 'arcobjects'.  No surprise, really.  I imagine that there is also a problem with the fact that the arcpy SpatialReference class invocation takes an argument of a WKID, spatial reference name, or feature dataset from which to copy the spatial reference; there isn't any way to create a new, blank spatial reference object.

Here's another thought: can one subclass an instance, as in:

class SubClassOfSpatialReferenceInstance(arcpy.SpatialReference(54004):

starting in this case with the World_Mercator spatial reference?

What I am ultimately after is a convenient interface for a modifiable spatial reference.  Contrary to ESRI documentation of the arcpy SpatialReference class, a majority of the properties are read-only.  I figured that, if I could subclass the arcpy Spatial Reference class, then I could override the name(), standardParallel1(), etc. attributes and use WKT string manipulation to make the spatial reference modifiable.  I suppose I could make a generic class that deals only with the WKT and implements a CreateSR() method that uses the loadFromString() method to make a spatial reference after all the properties have been set, but it just wouldn't be as slick.


texts: Filter null bands from Image Collection in GEE

texts:
I have created an Image collection (say Var MIC) of monthly mean data from a larger image collection (say Var LIC).
The created image collection MIC have some images with no bands due to no scene presents in initial LIC image collection.
How to filter out images from MIC with no image/bands.


texts: Understanding value range for Sentinel-2 images

texts:
I have a question regarding Sentinel-2 images.
I do not understand how the spectral reflectance values are saved/converted.
After the radiometric specification page:

The equipment senses the light intensity with a resolution of 12-bits
The integer light intensity values are then converted to a 16-bit reflectance integer value (How it is converted is not really specified)

Looking at the Sentinel-2 Product Specification one can see that:
Reflectance(float) = DC / QUANTIFICATION_VALUE (QV = 1_000 or 10_000)
This was also mentioned in a similar question on this forum.
But that doesn't explain how the light intensity is converted to the 16-bit reflectance value/Digital Count.
I was looking at some L2C images and saw that the provided values exceeded 2^12, which surprised me.
(Not including the Reserved max value)
In the previously mentioned post, they explained that they changed it to a 15 (16?) bit encoding scheme and increasing QV from 1_000 to 10_000. I thought the satellite provides the data as 12-bits and just saves them as a 16-bit integer, as 2 Bytes are used either way. So how could they increase the values up to 15-bits?
Was the "max" value increased? Did the domain stay the same? Or was some kind of scaling involved?
If the data previously used 2^12 bits and converted the values by dividing by 1_000 everything is fine for me. But how were they able to go up to 15-bits? Does that mean that the images are more detailed?
I am still very new to this field, so sorry if I am missing something obvious :)
Thanks!


texts: ArcPy CopyFeatures_management() naming error

texts:
I'm trying to reproject and copy some shapefile featureclasses to a file geodatabases and I'm having problems with the names imported. The code Im using is

# -*- coding: iso-8859-1 -*-
# -*- coding: latin-1 -*-

fcItem = "gsd_ap.shp"
outWorkspace = "database.gdb"
outCS = "1330"
fcList = arcpy.ListFeatureClasses() for fcItem in fcList:
    outFC = os.path.join(outWorkspace, fcItem.strip(".shp"))
    arcpy.Project_management(fcItem, outFC, outCS)

Unfortunately the outFC name in the file geodatabase becomes gsd_a, the last letter disappears. The same happens if the last letter is a "s". There are more letter combinations that don't work as well.

Why does this happens? CopyFeatures_management() does exactly the same.

Doing the above inside ArcGIS Desktop using toolbox, this does not happen, only when I try a standalone script.


texts:
Given an odometry, a direction and a point where I'm, I would like to change my current position according to these variables.

Let's say dLat and dLon are the difference of latitude and longitude between where I'm and where I would like to go. Let's assume the radius of the Earth equals 6,371,000.0 meters.

With papers and a pen, I found some equations to solve my problem : (ok, my background in maths could be better)

First of all, I calculate the degrees of each angle :

degreesAngleLat = 90 * dlat / (abs(dlon) + abs(dlat))
degreesAngleLon = 90 * dlon / (abs(dlon) + abs(dlat))

Then, I'll use some trigonometry to find the distance of the latitude and longitude :

odoLat = cos(degreesAngleLon * pi / 180) * odometry // I know, lattitude and longitude are reversed
odoLon = cos(degreesAngleLat * pi / 180) * odometry

Once I compute the direction, I want to know how many meters I have made in latitude and in longitude. I draw a rectangle, 1 side is my current position, the opposite is where I should be (with the odometry and the angle). So, I split it in 2 right-angle triangles. I know the hypotenuse, the angle, I can calculate the adjacent side. If I'm not wrong, cos(angle) = adjacent/hypotenuse.

After that, I convert this odometry to an angle :

deltaAngleLon = asin(odoLon/2/earthRadiusM) / pi * 180 * 2
deltaAngleLat = asin(odoLat/2/earthRadiusM) / pi * 180 * 2

Why * 2 ? Between 2 points on the earth's surface and the centre of the earth, an isosceles rectangle can be drawn. I use only half of the odometry to get a rectangular triangle (hence the / 2), then I multiply the angle obtained by 2 to match the correct odometry.

I have the odometry (or distance) for the latitude and longitude, And I would like to convert them into an angle. I draw my isoscele rectangle between the center of the earth and the 2 points in the surface. I draw 2 right-angle triangles in this triangle. I assume the segment is the size of the radius of the earth (in mathematics, I know it's wrong, but in this case, I think it can only add more accuracy ?) then, I have the size of the hypotenuse, the size of the side of the opposite angle (odometry / 2), I can calculate the angle : sin(angle/2) = odometry/2/radius.

Finally, I "just" have to add these delta to my current position (and modify the signs if needed).
Using maths, it's working... I thought...
Yes, it's working in the position 0.0, 0.0, but not in Europe...

Let's say :

odometry = 100
dLat = 1
dLon = 1

With these steps, I can find :

deltaAngleLat = 0.000636
deltaAngleLon = 0.000636

I used a website to calculate the distance, between the null point and this coordinate, it gives me 100m, perfect. (and an angle of 45°, nice)
But, when I add them to my current position, in France in example, I find 83.9m and an angle of 42.5...

My current position : 50.426753, 2.828529

With deltaAngle added : 50.427389, 2.829165

What am I doing wrong ? I know the radius isn't the same everywhere, but given this huge error, I don't think it's linked ?

This is my little program in Go : https://play.golang.org/p/znC98p7xAM9


texts: How do I properly import python sklearn to QGIS PyQt plugin?

texts:
I'm somewhat new to python and QGIS plugin development. I've created a plugin structure for QGIS using the Plugin Builder on QGIS and I've added some controls to the gui and managed to update the appropriate module binding to the gui.

I'm also trying to use sklearn library nearest neighbor in my plugin application. Everything works well until sklearn first runs and I get an error message as follows:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/Users/c_aaalam/.qgis2/python/plugins\FemtoCoverage\femto_coverage.py", line 282, in run
    distances, indices = nbrs.kneighbors(subs)
  File "C:/Users/c_aaalam/.qgis2/python/plugins\FemtoCoverage\sklearn\neighbors\base.py", line 332, in kneighbors
    return_distance=return_distance)
  File "binary_tree.pxi", line 1295, in sklearn.neighbors.kd_tree.BinaryTree.query (sklearn\neighbors\kd_tree.c:10386)
  File "C:/Users/c_aaalam/.qgis2/python/plugins\FemtoCoverage\sklearn\utils\validation.py", line 117, in array2d
    if sp.issparse(X):
AttributeError: 'NoneType' object has no attribute 'issparse'

This is the script I'm running in my plugin:

nbrs = NearestNeighbors(n_neighbors = 1, algorithm='auto').fit(array1)

distances, indices = nbrs.kneighbors(array2)

It looks like the scipy sparse variable is "NoneType" and therefore the "issparse()" function is not recognized. I have placed the sklearn folder in the plugin directory (windows user) and added the necessary import sklearn to my plugin module.

Does anyone have experience working with sklearn in PyQt applications or specifically with developing QGIS plugins and can help me resolve this error?

Maybe I'm loading sklearn incorrectly?


texts: How can I export a large zone from OpenStreetMap?

texts:
I would like to export San Francisco in OSM so I can use it in Blender.
The problem is that the zone is too big for OpenStreetMap so I have an error each time...

Do I have another way to export the whole city?


texts: Add prefix to properties names using ogr2ogr sqlite

texts:
I am using ogr2ogr to create a geojson file with props from sqlite Db with this
ogr2ogr -f GeoJSON -sql "SELECT * FROM censusblockgroups2 shp INNER JOIN 'output-z' csv ON csv.name = shp.censusblockgroup" joined_output.geojson myjoindb.db
It will create a geojson with this structure
{
"type": "FeatureCollection",
"name": "SELECT",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
{ "type": "Feature", "properties": { "id": 1, "attr": "ipsum" }, "geometry": { "type": "Point", "coordinates": [ -78.653, 35.7874 ] } },
{ "type": "Feature", "properties": { "id": 2, "attr": "lorem" }, "geometry": { "type": "Point", "coordinates": [ -78.6298, 35.7873 ] } },
{ "type": "Feature", "properties": { "id": 3, "attr": "ipsum" }, "geometry": { "type": "Point", "coordinates": [ -78.6408, 35.7795 ] } }
]
}

But is it possible to create a prefix for properties name depending on SQL table name? like
"properties": { "id": 3, "output-z_attr": "ipsum" } or by alias name
"properties": { "id": 3, "csv_attr": "ipsum" } or even like
"properties": { "id": 3, "csv": {"name":"ipsum", other props from table} }


texts: Can I create tiled Qgis2threejs 3D model

texts:
Is able to create 3D model with Qgis2threejs and tiled it for better loading over net?

For example 3D model will be tiled for 100 m X 100 m


texts: How to copy the color styling to a color column in a layer's attribute table?

texts:
I have a polygon layer in QGIS that I have applied a random color style to each polygon. Under Layer Properties -> Style I selected a Categorized style and then generated a random color ramp. I added a color column to the layer table. Is there a way to automatically copy the color assigned in the style to the color column for each polygon in the form "#ff0000".

Ultimately, I want to export it as a GeoJSON layer and import it into a leaflet map. The color column will set the color in leaflet.


texts: Listing feature classes in multiple geodatabases in folder using ArcPy?

texts:
I have a folder with numerous geodatabases in it. Within the geodatabases are feature classes. I would like to get a list of the feature classes in these geodatabases using a Python script. I've only gotten so far as listing the geodatabases. Does anyone have a suggestion on listing the feature classes in the geodatabases?

The print item line in my script gives me the following results:

C:\output\data.gdb
C:\output\otherdata.gdb
C:\output\somethingelse.gdb

I want to generate a list of the feature classes in the geodatabases above. Below is my script so far.

    import arcpy, os, sys
    from arcpy import env

    arcpy.env.workspace = "D:\\output"
    inWorkspace = arcpy.env.workspace

    workspaces = arcpy.ListWorkspaces("*", "FileGDB")
    for item in workspaces:
         print item #This part gives me the print statements I shared above
    # fcList = arcpy.ListFeatureClasses() #I haven't figured this part out
    #I want to list the feature classes in the geodatabases


texts:
I have a DB in SQlite with several layers. I would like to perform a spatial analysis, by buffering several water lines with 30m from origin.

I am using DB Manager in QGIS 2.14 and the SQL statement is:

SELECT PK_UID, buffer('Geometry',30)
FROM cursos_agua;

Here's the resulting column. As you can see, empty...

What am I doing wrong?


texts: How to program select by attribute unique values using ArcGIS silverlight API?

texts: I'd like to create a select by attribute window in my webGIS using silverlight api . this is like select by attribute in arcgis desktop . i create  select layer combobox ( cbolayers)and select fields combobox(cbofields) to select fields of selected layer  . the application needs a listbox to list unique values of selected field like arcgis desktop unique values listbox . i surf in the internet but couldn't find any good source to help me about How can i create unique values listbox .could you please help me about this subject?


texts: Working with two rasters which have different resolutions?

texts:
I have two rasters who correspond to the same place, one is a thermal orthomosaic and the other is a DEM, because each one was done with different cameras (the thermal cameras have low resolution) the cell size is different.

The rasters are georeferenced.

How could I work with these two rasters?

My intention is compute the height and the temperature at the same time adding an equation with the raster calculator, this means the cells need to have the same size and be in the same position.

Which is the tool in ArcGIS for that?


texts:
I uninstalled QGIS and removed all files on the computer with QGIS in them. I then reinstalled version 2.8. I still get the same problem that google streetmaps does not display. When I try one of the other options in the openlayers plugin I get this.

These little question marks suggest a problem with the graphics of this layer. How do I get Openlayers working?

EDIT: UPDATE: It works on a home pc not at work which has a proxy. We have tried fiddling with the proxy but to no avail.


texts: Displaying RGB raster on QGIS with no enhancement with Python Console?

texts:
I'm working with a RGB raster on QGIS 2.6.1, on a Windows 7 machine.

I need to set the Contrast Enhancement option to No Enhancement with the Python Console (not through the Layer Properties window).
I have developed a simple script, inspired by this discussion: QGIS Scripting problem updating contrast enhancements for raster layers using bandStatistics

layer = iface.activeLayer()
renderer = layer.renderer()
ContrastEnhancement = QgsContrastEnhancement.NoEnhancement

for Band in range(3):
    if Band == 0:
        myBand = layer.renderer().redBand()
    elif Band == 1:
        myBand = layer.renderer().greenBand()
    elif Band == 2:
        myBand = layer.renderer().blueBand()
    myType = layer.renderer().dataType(myBand)
    myEnhancement = QgsContrastEnhancement(myType)
    myEnhancement.setContrastEnhancementAlgorithm(ContrastEnhancement,True)
    if Band == 0:
        layer.renderer().setRedContrastEnhancement(myEnhancement)
    elif Band == 1:
        layer.renderer().setGreenContrastEnhancement(myEnhancement)
    elif Band == 2:
        layer.renderer().setBlueContrastEnhancement(myEnhancement)

layer.triggerRepaint()

Unfortunately, it does not work and the layer becomes not visible.
It seems that it is not handling correctly the min/max values of the bands. Any suggestion?


texts:
I am creating a database model with several table storing data from a data producer.
They produce their data with ArcGIS/QGIS and so produce an .xml file compliant ISO 19139.

We now want to store this data in a PostgreSQL Database and we would like to have the metadata in this system.

As far as I understand it, I would like to create a "metadata" schema (or table ?)
and store in an XML field format the content of the .xml file along with the name of the table, and then use it in the next part of the workflow.

Is there any guide about handling metadata in PostgreSQL ?
I looked around and didn't find a clear answer for me.

I already saw this post : Editing metadata in PostGIS database?


texts: Can't delete features from TAB layer table in QGIS if no geometry

texts:
I have some layers containing DAT,MAP and TAB files, when I loaded these layers into QGIS some layers have check boxes before layer name and some layers won't like showing in below image.

In the above image you can see that boundary layer have the check box and compare_layer shows a icon instead of check box.

Now the problem I am facing is I am not able to delete the features from the compare_layer table and also from the layers similar to this kind.

When am trying to delete a feature from the this layer table and trying to save the layer, it throws an error like Could not commit changes to layer and compare_layer: OGR error deleting feature 15:

How can I solve this issue?


texts: 16bit grayscale raster displays all black in Geoserver

texts:
I have a 16-bit raster (MrSID) that I want to display in Geoserver 2.0.2 but it is displaying as all black.  If I create a new style similar to the default raster.sld but with graychannel element the image displays correctly.  My question:  Is there a way to modify the default raster.sld so that it handles BOTH RGB and grayscale images (8 and 16-bit)?
Thanks


texts: I have a geoserver installed on a server and I want to move the geoserver to another server. What is the easiest way to do that?


texts: GML Filename gives layer name of "PREDEFINED" when imported as vector layer

texts:
I am new (today) to QGIS and I have hit a problem almost immediately. Basically, I am trying to read in a number of GML files, each with unique names. However when I import them into QGIS they get renamed to "PREDEFINED". The problem with this is that I want to use the python console to loop some functions, however, all the layers are named the same, so I can't do it.

I have a lot of GML files (345) and hence, to manually rename them is not an option.

I don't know what technical information you might need, I am using windows 8 and QGIS Desktop.


texts: QGIS 3.10 - Python problem with double for loop

texts:
I am new to Python and I am trying to write a custom function which iterates through a table of points and adds a field with all of the lines that touch it.
At the moment I am getting a partially right answer but there is some repetition going on.

An extract of the code:

 Homelayer = iface.activeLayer()
 Hfeatures = Homelayer.getFeatures()
 layerSet = {layer.name():layer for layer in qgis.core.QgsProject.instance().mapLayers().values()}
 targetLayerName = 'PRECABLE'
 targetFieldName = 'PF_PERE'
 for feature in Hfeatures:
    tst = []
    for feat in layerSet[targetLayerName].getFeatures():
        geom = feature.geometry()
        if geom.touches(feat.geometry()):
            run = feat.attribute(targetFieldName)
            tst.append(run)
            res = ','.join(tst)
            print(res)

The result I get is:

  XX1710

  XX1710,XX1711

  XX1710,XX1711,XX1712

  XX1710,XX1711,XX1712,XX1713

  XX1710,XX1711,XX1712,XX1713,XX1714

  XX1710,XX1711,XX1712,XX1713,XX1714,XX1715

  XX1710,XX1711,XX1712,XX1713,XX1714,XX1715,XX1716

  XX1710,XX1711,XX1712,XX1713,XX1714,XX1715,XX1716,XX1717

  XX1717

  XX1717,XX1721

  XX1717,XX1721,XX1722

  XX1717,XX1721,XX1722,XX1723

I just need the lines in bold.


texts: How to display images as cross sections in 3d

texts:
I want to display a base map and several vertical cross sections in 3D. The base map image is easy to rectify as a geotiff. How can I also display the cross sections?

The cross sections are at present just rectangular images. I can work out the depth/pixel conversion and the surface xy locations for each column of pixles. I have no DTM so the base map is flat and the cross sections simply reside directly underneath.

My preferred options are (1) anything written in Python (2) open source GIS software. eg QGIS, uDig etc. (3) perhaps point cloud software like meshlab?

Regards


texts: Why can't I select any of the tools under "edit"?

texts: Why can't I click any of the tools under "Edit?" I want to merge selected shapefile layers together, but all of the tools under "edit" are greyed out.


texts: How to correctly uncompress and re-compress an msd file using 7zip from command line

texts:
I am trying to uncompress and recompress both the msd and sd file so as to tweak some things I want.
In fact I want to get inside the layers folder and for each layer.xml file to string replace the the name of the dbconnection.
My goal is to create, with a bat file, 5 sd files genarated from an original sd file - for 5 users that are connected to the same oracle xe database, with the same password.
It looks ambitious but It works.  I tried channging the db names by hand, by getting inside the sd and the the msd packages, using 7zip following a method described by @ericchiasson here. In the same article @jpmc26 states that you can recompress and uncompress the sd, but you have to follow specific compression properties.
 @jpmc26 says:

  By the way, you can recompress an sd file. You have to use the right compression options. Options that work for me are Format: 7z, Word Size: 16, Solid Block size: Non-solid. I also set Dictionary Size to 64KB, but I don't think that matters. All other options were left to their defaults. (I made some educated guesses based on the fact 7-zip reports the compression method as "LZMA:16" for an sd file.) –  jpmc26 Jun 18 at

I followed his instruction and used 7z from command line like this: 7za a -t7z -m0=lzma -md=64k -ms=off %MAP_NAME_MXD%.7z
However these properties work for compressing the sd file.
Inside the .sd there is one compressed file  which is of .msd type. I try compressing it with the same configutaration, after having uncompressing that so as to tweak the .xml files inside of it.

When I try to publish my generated sd it is published correctly but it fails to start. I tried replaceying with the original msd the generated msd inside of my newly generated sd, and both service publish and start worked. This made me thinking that the method of compressing the msd package is the problematic one.

Can anyone tell me what are the correct compression options for the msd files ?
An alternative would be to avoid uncompression and recompression by injected in some way the tweaked layer xml files inside of the msd. But I dont know if that is possible


texts: Edge detection from Google Earth basemap in QGIS

texts:
I would like to detect boundaries of agricultural fields for 70 different polygons spread out across a country.
Is it possible to run an edge detection algorithm using the google earth (or equivalent) basemap in QGIS to avoid having to download many separate satellite imagery files?


texts: Using CropCoverage Render Transformation in Geoserver SLD (WPS)?

texts:
I have yet to find an example on the web (or in code) of anyone using the "Crop Coverage" Render Transformation from a Geoserver SLD.  Does anyone have an example or at least be able to point me to the syntax of using an existing Geoserver hosted layer as a the "crop geometry?"

gs:CropCoverage is one of the WPS Processes used to crop raster layers.

This answer to Geoserver Raster layer cropped by polygon? provides an example of using gs:CropCoverage with a literal crop geometry, but I need to use an existing GeoServer layer as the crop geometry.

It is also possible to chain more WPS processes in a SLD but I cannot find an example of how to use gs:CropCoverage in a SLD.

At the moment I am using Geoserver 2.4.3 on Tomcat 6.


texts: Leaflet.markercluster, markerClusterGroup into another markerClusterGroup

texts:
I would like to add markerClusterGroup into another markerClusterGroup. I tried but without success, no marker are displayed.

Why? Because I would like to separate markerClusterGroup for multiple zones. (For each postal code of each country of each continent).

Or maybe there is another way to achieve this?

Edit:
Here is a example where I try to add multiple cluster inside a main cluster:
https://codepen.io/Inglebard/pen/jOWENgO

var marker_cluster_test= [];
var main_cluster_test= new L.MarkerClusterGroup();
var clusters= [];
for (i = 0; i < 5; i++) {
  clusters[i] = new L.MarkerClusterGroup();
  clusters[i].addTo(main_cluster_test);
  marker_cluster_test[i]=[];
  for (j = 0; j < 5; j++) {
    marker_cluster_test[i][j] = L.marker([49.0435987,2.0838575]) ;
    marker_cluster_test[i][j].addTo(clusters[i]);
  }
}
main_cluster_test.addTo(map);


texts: Is it possible to reproject KML using ESRI's Silverlight API?

texts:
My question is similar to this one, but for Silverlight.

I need to reproject an entire KML file to NAD83, add it to the map with all the layers, symbology, etc. in tact.

Is there anyway to do this with the Silverlight API?


texts: Scripted map rendering (multiple similar maps from one shapefile) using PyQGIS?

texts:
I would like to ask for advice on creating a Python script that will provide the opportunity to render a plenty of similar maps (with the same extent, projection, layout, color scheme and composition, but with differing titles and data used for rendering).

Say we have a single shapefile or Spatialite DB with a huge attribute table inherited from csv file like this:

I want to turn them into a set of choropleth maps like this:

So I guess that I need a script that will just:

Iterate through specified range of columns in the attribute table;
Call some (yet not written?) rendering function that will produce the choropleth map with a legend and title;
Dump the output to the vector formatted-file (PDF, EPS, etc.).

I'm going to use QGis for that purpose and would be glad if anyone could outline the main ideas and key QGis API functions that I can rely on at first.


texts:
I am not sure whether this question belongs here. I think it belongs more to Geo-technical, but anyway, I just give it a shot.

During earthwork development, earth terrain has to be cut ( if the terrain is higher than platform) or fill ( if the terrain is lower than platform) in order for the construction of platforms such as roads, housing area, pond and whatnot. Usually, one would have to make sure that the cut and fill within a development area is roughly the same, in order to reduce the import or export of the soil.

The process of calculating the height of the platforms with relative to the terrain so that the cut/fill is roughly the same is called auto-balancing.

How do engineers ( or earthwork software) do the auto balancing in practice?


texts: Use FME to delete records from PostGIS before loading new records

texts:
I have some shapefiles. Each shapefile contains the data on a topic for a county. I have data for 26 counties.

In the filename of the shapefile, the county ID is given.

I load these shapefiles into PostGIS using FME. The data for all counties is loaded into one table (so 26 shapefiles get loaded into 1 PostGIS table).

This data gets updated every now and then, and I can request the data for the county where the data has been changed. I get supplied with a new shapefile with all data for that county (not just the changed records).

I would like to delete the data in the PostGIS table for the county that I have the new data for before I load the new data.

How do I do this? Is this possible if I get the county ID from the filename of the shapefile and each record in the PostGIS table has a county ID attached to it?

My plan is to just add the new shapefile to a specified folder, and run my FME workspace.

Edit 31/05/2013

Just to clarify, I have multiple records for each county (i.e. all roads, buildings etc in a county). I basically want to bulk delete a counties data, and replace it with the new data.


texts: ogr2ogr import fgdb to postgis and prefix added table name

texts:
I have multiple fgdb with many feature classes.

I am importing them to my PostGIS and this is the code i used

  ogr2ogr -overwrite -progress -f "PostgreSQL" PG:"host=localhost port=5432 dbname=postgres user=postgres password=postgres^" MY.gdb

Yet the table names in my PostGIS comes same as my fgdb feature classes.
I want to add prefix when I am importing them to PostGIS
like "buliding" feature class in FGDB becomes Test_building in PostGIS

"-nln" option does not seem to have that function (or I just can't find?)


texts: Is there a coordinate system specifically designed for oblique map viewing?

texts:
I've long been impressed with maps that have an oblique perspective that also provides the viewer with a sense of the earth's curvature.  National Geographic maps for decades often have had this type of view, wherein scale varies as the earth recedes into the distance.  Similarly, here's a 2013 example from the New York Times:

(Source)

I get the sense that I'm looking down from a spaceship.  Google Earth gives a similar impression:

All of my previous mapping has used a traditional vertical view, so I thought that I would expand my cartographic skills by learning how to create that "oblique, bulging-earth effect" with Arc Desktop, using the New York Times map above as a reference.

I couldn't find any advice on how to proceed, so I experimented with a variety of coordinate systems.  the best that I came up with was the Vertical Perspective (World) (ESRI also refers to it as the World Vertical Perspective):

(For this example I modified the Longitude of Center to -150.0 to give the map the off-center, angled view towards the northeast.)

I think that this is not too bad.  I included latitude and longitude lines to confirm that the image indeed curved and bulged.

This was an interesting process which helped me to see how a variety of coordinate systems (and their settings) affected the map.

Now that I've found something that appears to work I'm wondering... are there any coordinate systems specifically designed to create the "oblique, bulging-earth effect" that I'm looking for?

It occurs to me that there may be some who object to my use of the phrase "Coordinate System", preferring that I use "Spatial Reference" or something else.  To fend off such diversions, I am merely following ESRI's nomenclature, as pictured below.


texts: Web AppBuilder - Trigger field calculation when application is closed

texts:
We have a public facing app intended for the public to go and add points to map. The feature service they are editing is public facing. I only want points visible that the user creates during their session.
For example, user A comes in and creates a few points. Then user B access the app but does not see any points till they create their own.
One possible solution I've come up with is to have a field populated with 'Visible' and 'Not Visible'. Using a def query in the map I can hide all points till they are populated with 'Visible'. Obviously Visible will have to be the default for users to see the points they create. But was wondering if it was possible to calculate all points as 'Not Visible' after the user closes the application. The data is coming from an enterprise geodatabase on SQL Server. Is this possible using Web AppBuilder or WAB dev edition?