texts: Iterate over features and years to export images for each feature-year

texts:
For each year, I need to take the max. of spring months from the representative rasters, and then clip the resulting "max raster" to each region in a feature collection. Ultimately, I need to automate the export of each of the resulting rasters (i.e. a raster for each year and each feature). So let's say...

feature_list is my Feature Collection
raster_collection is the Image Collection from the online catalog

var years = ee.List.sequence(1984, 2015);

I assume I have to iterate over the years, such that I can...

year = ee.Number(year).toInt();
var startdate = ee.String(year).cat("-04-01");
var enddate = ee.String(year).cat("-06-30");
var WaterSpring = raster_collection.filterDate(startdate,enddate).max();

Then I need to clip by each "feat" in feature_list...

var WaterSubset = WaterSpring.clip(feat);
var Water = WaterSubset.gt(1);  // make Water binary
Water = Water.updateMask(Water.neq(0));

I have found both these questions...

Iterating over years for features in feature collection using Google Earth Engine?

How can I export a set of images from Google Earth Engine?

I can't quite figure out how to "reverse engineer" them though because it seems like I might be getting a multi-dimensional list as the output, or something like that.

UPDATE: If this makes the wording/concept simpler, I basically need to iterate or map over a Feature Collection and an Image Collection, such that I have a separate clipped image for each feature, for all images in the Image Collection.


texts: How do I convert an elevation shapefile into a CSV file with its XYZ data?

texts: I am new to using QGIS and LeapfrogGeothermal Modeling.  I have a shapefile for Elevation and Altitude data and need to extract XYZ coordinates out of it? Can I do that? and how? and what kind of file do I need to change it to, to be able to use it?  I'm attempting to use this shapefile to create a solid topography surface but need the coordinates to upload into the program.


texts: API for getting only "essential" OSM source data for large zoom levels?

texts:
I'd like to implement a custom client-side rendering for OpenStreetMap. For that, I'd need the source data to render. I could use the OSM XML source - but not at higher scales, as the amount of data becomes too high. And most of that data is not needed, as it makes sense to show fewer features at higher zoom levels.

Size of already rendered data is roughly proportional to (bbox area)/zoom^2. Is there an API, and a server providing OSM source data with it, which would also return only about a fixed amount of data for a "fixed viewport size"?

Ideally, I'd like to select the features which I want to see at particular zoom levels, but a server which selects it for me would be also a great start.


texts: How to highlight leaflet marker cluster on mouse over?

texts:
Is there a way to highlight a leaflet marker cluster when mouse over as you can do for a simple marker?

// ADD YOUR BASE TILES
var baseLayer = L.tileLayer('https://dnv9my2eseobd.cloudfront.net/v3/cartodb.map-4xtxp73f/{z}/{x}/{y}.png', {
    maxZoom: 5
});
// DEFINE THE CLUSTER LAYER
var markers = L.markerClusterGroup({
    showCoverageOnHover: false,
    maxClusterRadius: 10,
    iconCreateFunction: function (cluster) {
        var markers = cluster.getAllChildMarkers();
        return L.divIcon({html: '
' + '
' + cluster.getChildCount() + '
', className: 'mycluster', iconSize: L.point(30, 30) });
    },

});

// CALL THE CARTODB SQL API HERE IN URL FORMAT
$.getJSON('http://soundgoods.cartodb.com/api/v2/sql?format=GeoJSON&q=SELECT * FROM soundgoods_mixtape_map_1', function(data) {
    var geojsonMarkerOptions = {
        radius: 8,
        fillColor: "#FFCC00",
        color: "#323232",
        weight: 2,
        opacity: 0.5,
        fillOpacity: 0.4
    };

    var geojson = L.geoJson(data, {

        pointToLayer: function (feature, latlng) {
            var popupOptions = {maxWidth: 500};
            var popupContent = '' +
                    feature.properties.soundcloud +
                    '
' + "Post & tracklist" + '
' + '';
            var circle = L.circleMarker(latlng, geojsonMarkerOptions).bindPopup(popupContent,popupOptions);

            // Highlight the marker on hover
            circle.on('mouseover', function(){
                circle.setStyle({ color: '#ff8d00'});
            });

            // Un-highlight the marker on hover out
            circle.on('mouseout', function(){
                circle.setStyle(geojsonMarkerOptions);
            });

            return circle;
        }

    });
    markers.addLayer(geojson);

    // CONSTRUCT THE MAP
    var map = L.map('map').setView([0,0],3);
    baseLayer.addTo(map);
    markers.addTo(map);
});


texts: Is there a maximum number of polygons that can be plotted using Folium/Leaflet.js from a JSON file?

texts:
I have 15,000 parcels that I am trying to plot on a Folium map that is being served via Flask. I've imported the geojson as a layer, but when I include all 15,000 polygons, the map doesn't load. If I downsize that to 4,000 polygons, it loads great.

Is there a limit to the number of polygons that can be rendered? Are there any ways around this limit?


texts: Is it possible to not show node in editing mod QGIS

texts:
I actually have a plugin that allows me to draw on a vector layer (I made it based on the FreeHand Editing Plugin). When a mouseMoveEvent is detected, it creates a new point on the vector layer and create a line between it and the last point created. My problem is that I want to show in editing mod only the lines and not the nodes. Does anybody know if this is possible?

EDIT:

This is not what I meant, sorry if you didn't understand well, I'll explain better.

I have a vector layer (and an option in my plugin to create one)
In edit mod, I can freely draw lines with the mouse on this vector (the lines will correspond to the drawing, check the screenshot below when it's not in editing mod).

The problem is that in editing mod, I don't want to display the nodes cross image (check screenshot below) on the screen like it actually is, I just want to show the lines like how it is in non-editing mod


texts: Map exported to PDF background color is black ol5

texts:
When I am exporting map to PDF, it is all looking black in backgroud.

Here is my code

map.once('rendercomplete', function(event) {
    $('.filter').css('display','none');

    var canvas = event.context.canvas;
    event.context.fillStyle = "white";
    canvas.lineWidth = 2;
    canvas.strokeStyle="#FF0000";
    //showText(event.context,'left-top',dist,tahsil,village);
    //drawLegend(event.context,'right-bottom');
    drawNorthArrow(event.context,'right-top');
      WriteScaletoCanvas(event);

    var data = canvas.toDataURL('image/jpeg',1.0);
    var pdf = new jsPDF('landscape', undefined, format);

    pdf.addImage(data, 'JPEG', 0, 0, dim[0], dim[1]);
    pdf.save(pdf_name);
    document.title=pdf_name;
    window.print();
    $('.filter').css('display','block');

    map.setSize(size);
  // map.getView().fit(extent, size);  //edit2: not needed
    mapView.setZoom(currZoom);  // edit2: set original zoom
    exportElement.disabled = false;
    document.body.style.cursor = 'auto';
});

Here is the image

i need to have white color in background.

by applying this code

 event.context.fillStyle = 'white';
    event.context.fillRect(0, 0, event.context.canvas.width, event.context.canvas.height);

result is


texts:
I have thousands of Raster data concerning the ice coverage at the pole at different times (days) over several years. Now I would like to perform a cluster analysis. I want to split the data into different periods (e.g. 4 seasons) with respect to the ice cover.

Can I perform such a cluster analysis in R? And if, how? (I am a beginner in R, so please keep it as simple as possible. Thanks.)

Or is it possible to use Grass GIS or QGIS (considering the large amount of data)?


texts: PyQGIS3 retrieving points from a shortestpath graph

texts:
I'm upgrading some code which worked with QGIS 2.18 to now work with QGIS 3.4.5 and had the following snippet working to retrieve a set of points from a graph created by the QgsGraphBuilder function:

p = []
while (idStart != idStop):
    l = tree.vertex(idStop).inArc()
    if len(l) == 0:
        break
    e = tree.arc(l[0])
    p.insert(0, tree.vertex(e.inVertex()).point())
    idStop = e.outVertex()

(With idStart & idStop being tiedpoints from the graph)

This now does not work, I understand that inArc has changed to inEdges, so my revised code is as follows:

p = []
while (idStart != idStop):
    l = tree.vertex(idStop).inEdges()
    if len(l) == 0:
        break
    e = tree.edge(l[0])
    p.insert(0, tree.vertex(e.inVertex()).point())
    idStop = e.outVertex()

. However this gives me an error:
 AttributeError: 'QgsGraphVertex' object has no attribute 'inEdges'

This seems to contradict the information from the API.


texts: Correct usage of a counter using arcpy.da.UpdateCursor

texts:
In the Python code inside a Calculate Value tool in my model, I have the following:

getNum("%Selected Features%")

def getNum(ftr):
  lst = [4, 11, 15] #values for SQL expression
  x = 20
  fields = ('TYPE_A', 'TYPE_B', 'NUM_ID')
  for l in lst:
    y = x + 1
    whereclause = """"NUM_ID" = %s """ % l
    with arcpy.da.UpdateCursor(ftr, fields, whereclause) as cursor:
      for row in cursor:
        if row[0] ! = row[1]:
          row[2] = y
        cursor.updateRow(row)
    y += 1

What I would like to do is, for every element in lst, create a cursor based on a selection of records where NUM_ID = l. The values in TYPE_A and TYPE_B should be compared for each record and updated if necessary.

The script is failing after the cursor is created, with the error Item not in collection. If I remove the whereclause, it executes. Once the records have been selected and the values compared and updated if necessary, I need to increment the counter y by 1. I think I'm doing it in the wrong place though, because if I run the script without the whereclause it does execute correctly.When I check the records with the changed values, they are now all set to y, which indicates that y was not incremented.


texts: Translate Cartesian coordinates to WGS84, preferably in pyproj?

texts:
I have some Cartesian coordinates:

x = 3862525.511
y = -109738.522

I need to translate them into WGS84 coordinates, preferably using pyproj. Is this possible, or have I misunderstood the nature of Cartesian coordinates? I don't have any height info.

cartesian = ?
wgs84 = pyproj.Proj(init='epsg:4326')
new_x, new_y = pyproj.transform(cartesian, wgs84, x, y)

In case it's useful, these Cartesian coordinates weren't created by me, but I'm pretty sure they were created by copying and pasting OSGB36 eastings and northings on this site.

If it can't be done in pyproj, I can use gdal.


texts: How to print/export QGIS map composition as PDF with Google baselayer in it without spatial shift?

texts:
I encountered a problem in QGIS 2.18.21 with printing or exporting my Map Composition into PDF.
The layer composition has:

one vector polygon layer (either 4326 or 3857)
one baselayer (either Google Road/Terrain or OSM)

I tried all sort of combination:

print and export
my polygon in 4326 and 3857

But as long as the baselayer used is Google (either Road or Terrain) the resulting PDF has shift between polygon layer and Google baselayer, as shown in the picture below:

Looking closely into the shift, it appears that the polygon layer and Google baselayer has different scale, as going eastward and southward the shift grow bigger.

With OSM as baselayer, the PDF looks fine (no shift), as shown in the picture below:

Now to the question: how do I print/export my map composition with Google baselayer in it?


texts: Conversion of (lat, lon, height) to another CRS with pyproj results in no transformation

texts:
I have latitude, longitude, and ellipsoidal_height measurements in the IGS14 CRS. I want to convert/transform/project this data into latitude, longitude, and ellipsoidal_height measurements in the WGS84 CRS.  How do I do this with pyproj in python?
If I understand correctly, this is a transformation from EPSG:9018 to EPSG:4979. However, I am getting the same output as I input with the following:
from pyproj import CRS
from pyproj import Transformer

igs14 = CRS.from_epsg(9018)
wgs84 = CRS.from_epsg(4979)

lat = 78.592836594
lon = 18.189419118
e_h = 556.290578

crsTransformer = Transformer.from_crs(igs14, wgs84)
crsTransformer.transform(lat, lon, e_h)

>> (78.592836594, 18.189419118, 556.290578)

I checked my igs14, wgs84, and crsTransformer objects to find the following:
<Geographic 3D CRS: EPSG:9018>
Name: IGS14
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: IGS14
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

<Geographic 3D CRS: EPSG:4979>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- h[up]: Ellipsoidal height (metre)
Area of Use:
- name: World (by country)
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

<Transformation Transformer: noop>
Description: Ballpark geographic offset from IGS14 to WGS 84
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)

igs14 and wgs84 are clearly different CRS (with different datums and different ellipsoids), so why is the transformer a noop?
Thank you so much!
Python: 3.8.3 via anaconda
pyProj: 2.6.1 via anaconda


texts:
In my application, I am consuming a WMS layer. However, I would like to be able to gracefully fail if the service doesn't respond to a request. Instead of a pink tile, I would like a transparent tile and a message displayed to the user saying that the layer was not loaded properly.

What I am currently doing:

1) Setting the visibility of the .olImageLoadError div to hidden:

.olImageLoadError {
    visibility: hidden !IMPORTANT;
}

2) Creating a function that uses jQuery to see if there are any divs that have the class "olImageLoadError." It displays a message if one or more tiles fail to render. I link this function to the loadend event of the layer I am concerned with.

myWMSLayer.events.register('loadend', this, onloadstart);
function onloadstart(evt) {
    if ($(".olImageLoadError").size() > 0) {
        alert("Some assets failed to load");
    }
}

In searching for possible solutions, I found this (and several similar) site. It discusses the OpenLayers.Util.onImageLoadError property, which is a function that is called when a tile fails to load. However, this is no longer available in the newest version of OpenLayers.

Is there something similar in the current version of OpenLayers (at time of writing, 2.13.1)?


texts: Removing pixel values below -1 using ENVI?

texts:
I have downloaded NDVI from MODIS website and pixel values go up to -3000 to 1. I want to remove all values below -1. How do I do this?

I am using ENVI 4.3


texts: Calculate distance matrix based on network

texts:
I have the following layers in QGIS:

origins points
destination points
a layer of lines that is the route

I want to calculate the matrix with the distance from origin to destination using Network Analysis.

I've tried to do a python script using this, but I'm not seeing any results.

I honestly do not know anything about gis or python. But I need to make a data analysis. This is what I have done, I can not even read the layers.
I'm also not sure if the cost represents the distance traveled

#Definition of inputs and outputs
#==================================
##[my scripts]=group
##origin_layer=vector
##origin_id_column_index=field origin_layer
##destination_layer=vector
##destination_id_column_index=field destination_layer
##route_layer=vector
##output=string

#Algorithm body
#==================================

from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

from processing.core.VectorWriter import VectorWriter
import processing
from scipy.spatial.distance import cdist
import numpy as np
from math import sqrt

origin_layer = processing.getObject(origin_layer)
origin_id_column_index = origin_layer.fieldNameIndex(origin_id_column_index)

destination_layer = processing.getObject(destination_layer)
destination_id_column_index = destination_layer.fieldNameIndex(destination_id_column_index)

route_layer = processing.getObject(route_layer)

# don't use information about road direction from layer attributes,
# all roads are treated as two-way
director = QgsLineVectorLayerDirector(route_layer, -1, '', '', '', 3)

for origin_feature in origin_layer:
    pStart = origin_feature.geometry().centroid().asPoint()
    for destination_feature in destination_layer:
        pStop = destination_feature.geometry().centroid().asPoint()
        tiedPoints = director.makeGraph(builder, [pStart, pStop])
        graph = builder.graph()
        tStart = tiedPoints[0]
        tStop = tiedPoints[1]
        idStart = graph.findVertex(tStart)
        idStop = graph.findVertex(tStop)
        (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)


texts: Filtering Label in SLD based on Elevation

texts:
How can I filter the label in SLD, I want to have a label for all the contour-lines where the elevation is between 200-400. I tried with this SLD, but it doesn't apply anything.

<Rule>
  <Name>rule 6</Name>
  <Title>Label (200-400)</Title>

  <ogc:Filter>
    <ogc:PropertyIsBetween>
        <ogc:Function name="IEEERemainder">
            <ogc:Function name="int2ddouble">
                <ogc:PropertyName>value</ogc:PropertyName>
            </ogc:Function>
            <ogc:Function name="parseDouble">
                <ogc:Literal>50.0</ogc:Literal>
            </ogc:Function>
        </ogc:Function>
        <ogc:LowerBoundary>
            <ogc:Literal>200</ogc:Literal>
        </ogc:LowerBoundary>
        <ogc:UpperBoundary>
            <ogc:Literal>400</ogc:Literal>
        </ogc:UpperBoundary>
    </ogc:PropertyIsBetween>
</ogc:Filter>

  <MaxScaleDenominator>120000</MaxScaleDenominator>
  <TextSymbolizer>
    <Label>
      <ogc:Function name="round">
        <ogc:PropertyName>value</ogc:PropertyName>
      </ogc:Function>
    </Label>
    <Font>
      <CssParameter name="font-family">Arial</CssParameter>
      <CssParameter name="font-weight">Normal</CssParameter>
      <CssParameter name="font-size">12</CssParameter>
    </Font>
    <LabelPlacement>
      <LinePlacement/>
    </LabelPlacement>
    <Halo>
      <Radius>
        <ogc:Literal>2</ogc:Literal>
      </Radius>
      <Fill>
        <CssParameter name="fill">#FFFFFF</CssParameter>
        <CssParameter name="fill-opacity">0.9</CssParameter>
      </Fill>
    </Halo>
    <Fill>
      <CssParameter name="fill">#ff0000</CssParameter>
    </Fill>
    <Priority>2000</Priority>
    <VendorOption name="followLine">true</VendorOption>
    <VendorOption name="repeat">300</VendorOption>
    <VendorOption name="maxDisplacement">50</VendorOption>
    <VendorOption name="maxAngleDelta">30</VendorOption>
    <VendorOption name="spaceAround">30</VendorOption>
  </TextSymbolizer>
</Rule>


texts: How to make a .pgw file for a .png depicting a 2D profile

texts:
My question is more of an inquiry about whether I am trying to do is possible. I would like to reference an image depicting a 2D Georadar profile. I know that world files are usually used for georeferencing but I am curious about whether I could do it for a vertical slice. Now, the pixel dimensions of my picture are 3715x356 and the coordinates are:
top left corner: 0, 8.53
bottom right corner: 628.5, -9.17
Remember these coordinates indicate horizontal distance along the profile and depth in relation to sea level (negative value means below sea level).

Is it possible to create a .pgw file for this image and how would it look like? What is the norm for making one? I tried to make one myself but the .pgw file although recognized, doesn't affect the coordinates of the referenced image. The referenced image I get has the pixel dimensions as length per axis.


texts: Does vertical datum get projected? I feel a bit confused with the concept of the vertical datum and doing a projection.


texts:
After import some data from OSM file into Postgre, the coordinates numbers are too big and having trouble to cast to Geography:

Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

In my table:

select ST_AsText( the_geom )  from routes limit 4

"LINESTRING(-4902130.86 -2265840.91,-4902083.04 -2265907.16)"
"LINESTRING(-4903296.22 -2266295.18,-4903363.69 -2266333.39)"
"LINESTRING(-5705594.67 -3518944.1,-5705523.22 -3518937.37)"
"LINESTRING(-5484056.84 -1884629.16,-5484093.47 -1884510.24)"

I think the numbers should be between (-90 / 90, -180 / 180) ...

Info:

1: select ST_SRID( the_geom )  returns 900913 for all lines.

2: select UpdateGeometrySRID('public', 'routes', 'the_geom', 4326) does not converts the numbers, but SRID now shows 4326.


texts: Can an ArcGIS Addin be aware of whether it was installed from a shared folder or installed though the the Addin file?

texts:
I am distributing an Addin for ArcGIS, and I love when you point the Addin manager to a folder, updates automatically get applied on the next restart. However, when you click on an Addin, it takes away this ability because it doesn't know where the source is.

In my Addin, I would like it to figure out if it is installed the shared way, or the double-click-on-the-file way. Then if it is installed in the double-click way, I can prompt the user if they would like to upgrade and show them the path to the upgrade.

It appears the Add-In manager knows where the Addin came from because it shows "My Add-ins" (which are from the double-click method) and the "Shared Add-Ins" (the add the folder method).

Anyone have any ideas which ArcObject this may be obtained or another round about way?


texts: Converting only one color in a raster to polygons

texts:
I would like to use the image of the map above to convert all the dark purple color (ie existing manufacturing land) to polygons.  Other shapes and colors I would like to ignore.  Furthermore, this is a map of Chicago, so I would like to be able to superimpose the jpg on the location of Chicago and then have the polygons created in the correct location.  I don't mind trying to do this either in ArcGIS or Python.  However when I try the basic tools like Raster to Polygon I don't have an option to select the color.  Also I don't know what my license is but many of the options in my toolbox (like Spatial Analyst) say that I don't have a license for it.

Update:

I have reclassified and converted the raster to polygons.  Unfortunately as I didn't know what the values meant in the original reclassify I only used the Red Band as the values in that seemed to correspond with the Red value in the RGB value.  After converting to polygons I got something like this:

As you can see there are big polygons in many of the right places (I superimposed the two pictures) but lots of small dots in other places especially where there are black (but not deep purple) lines like streets and boundaries.  Can I do any better?


texts: Can concurrent processes be run in a single model?

texts:
In model builder (for Arc GIS 10) is it possible to run two or more geoprocessing operations concurrently?

For example I have to run repair geometry on 2 files and then merge them. Currently it repairs topo of file A and then file B and then merges them. Each file takes about 30h to repair. Can I get it to concurrently repair A & B, wait for both to be complete and then merge them.

At the moment I have to open two instances of arc catalog and run the repair tool on a file in each and once complete do a merge which means that I have to come in to office to check when both are complete...on a quad core, each repair process only uses 1 core at a time.

The above example is simple but if multi-threading is possible where multiple operations are run concurrently (if they don't depend on each other) and wait for the depended files to be complete before the next process is run.

Thanks,


texts: Is it possible to remove small transparent outline effects between polygons?

texts: after using the vectorizing tool I get the following result (picture). Around the cells you can see a small "boundary". This thin line is transparent and annoying. Is it possible to remove these small transparent areas?


texts: Adding "wikimedia unlabelled layer" as basemap option in ArcMap

texts:
In QGIS, I have used the OpenLayers plugin to make a number of basemaps, including the Wikimedia labelled and unlabelled layers available. I am seeking to find a way to use these Wikimedia layers as basemaps in ArcMap.
Is there a way this layer can be accessed from ArcGIS Online or could somebody show me how to download tiles?
I am using ArcMap 10.8.


texts: Use Ushahidi api in another map using openlayers

texts:
I have a ushahidi deployment and i am using Openlayers + ushahidi api to show the reports in a separate map along with other features.

i have queried the data from ushahidi api and made a vector layer. the layer has data but it not showing up.

This part gets the data from ushahidi

var roadreportrequest = OpenLayers.Request.GET({
url: '../roadreport/api/?task=incidents',
async: false,
success: function(Request) {
    // alert("success");
    roadreportjson = JSON.parse(Request.responseText);
},
failure: function() {
    alert("failed to load data from roadrepoter");
    roadreportjson = {}
}

});

and this part converts it to vector layer

var roadreport_vector = new OpenLayers.Layer.Vector('RoadReporter Reports');

for (incident_index in roadreportjson.payload.incidents) {
// console.log(incident_index, roadreportjson.payload.incidents[incident_index]);
incident = roadreportjson.payload.incidents[incident_index];
roadreport_vector.addFeatures([
    new OpenLayers.Feature.Vector(
        new OpenLayers.Geometry.Point(incident.incident.locationlongitude, incident.incident.locationlatitude),
        incident)
])
}
map.addLayer(roadreport_vector);

There is data in the layer, i have checked the console and this code to pan the map also works.

map.panTo(new OpenLayers.LonLat(roadreport_vector.features[2].geometry.x, roadreport_vector.features[2].geometry.y))

My issue is that i am expecting the data to be shown up as on the map(the default orange circles), but it is not showing up.


texts: Saving vector layer into KML across 180 longitude line

texts:
I am trying to save a vector layer into KML, however in this case some of the shape cross the 180 longitude line and the shape just breaks and become long lines.. (they are supposed to be circles)

How do I save it so that the shape stay as the original image (2nd image)?
I would like to export the file into google map with the information in the tables.


texts:
I am trying to query all the nodes which are in Potsdam city and with postcode 14480.
I know that it can be coded as follows:
    area["ISO3166-1"="DE"][admin_level=2]->.germany;
(
  node["addr:city"="Potsdam"]["addr:postcode"="14480"](area.germany);
)->.r;

but I would like to use an if condotion with evaluators as shown below in the code.
However, it seems that there is something wrong in the code below
How can I fix it?
Code:
    [out:json][timeout:180];
area["ISO3166-1"="DE"]["admin_level"="2"]->.germany;
(
  node(
    if: (t["addr:city"] == "Potsdam") && (t["addr:postcode"] =="14480")
  );
);

  out geom;


texts: Getting odd values in ArcGIS attribute table?

texts:
I'm very new to/teaching myself ArcGIS.

The goal:  extract raster data using multi values to point

The process:  In short, I overlaid lat/long coordinate data on top of a raster, and used the Spatial Analyst tool to extract the values from those points.  The raster contains data related to percent wet nitrogen deposition, so the values should be between 0 - 100.  Prior to overlaying the coordinates on the raster data, I changed the raster datum to NAD83, to match the coordinate datum. I also changed the projection to match as well.

The issue:  Many of the cells in the resulting attribute table contain this value: "-3.40282E+38".  What does this mean?


texts: Open source tool to view DTED and DEM raster-data

texts:
I am new to GIS field. Is there any open source GIS tool to view DTED and DEM tiff-files?

Update 1

I have some DTED files and DEM files. I want to open and see those files using some open source GIS tools. Can I use MapWindow GIS or QGIS for this? Are there any other tools available?

Please advise.


texts: I'm not able to create a line in ArcGIS. I started an edit session and began with straight segment. All I want to do is draw a line out of three single points - lets say the coordinates for the points are (20,30 ; 50,30 ; 70,50). Every time I finish sketching, ArcGIS is creating a polygon instead of just creating a line without closing it.


texts:
I'm using the mapbox library and I have this function that processes a response of type json, but I want the color of every single feature to be dynamic based on a property of the feature itself:

function doSome(data) {

    map.addSource('hrd-src', {
        'type': 'geojson',
        'data': data
    });

    map.addLayer({
        'id': 'hrd-polygons',
        'type': 'fill',
        'source': 'hrd-src',
        'layout': {},
        'paint': {
            'fill-color': '#088',
            'fill-opacity': 0.8
        }
    });

    // Set filter to first hour of the day
    // 0 = mezzanotte del giorno corrente (0-23)
    filterBy(0);

    document.getElementById('slider').addEventListener('input', function(e) {
        var hour = parseInt(e.target.value, 10);
        filterBy(hour);
    });
}

Is it possible to set a dynamic value for the fill-color property?


texts: Multiple layer line selection based on multiple layer point, then calculate the length of lines

texts: I need to measure water network from a source to multiple meters. Source point is connected to water main line, then water main line is connected to service line, service line is connected to meters. i have 4 shapefiles.......... source point, meters, water main line and service line. I have to calculate distance from source to each meter individually. can anyone help me in VBA (ArcGIS 9.2 or 9.3).


texts: How to improve the display of color images?

texts:
I've found some really nice color ramps for portraying rasters, but in ArcGIS and QGIS, they don't seem to be easily usable for a finished product because their complexity causes obvious steps rather than a smooth gradient.  The raster renderer appears to be stuck at an 8-bit palette (only 256 levels), rather than a mathematical transform involving 24-bit standard RGB or deep 48-bit RGB.

There's probably an awkward workaround involving splitting up the raster into multiple layers and symbolizing them separately, but I'd like to know which desktop and web GIS systems can present stretched/colorized gradients like this natively in a higher precision format.


texts:
I'm new in QGIS and python, using QGIS 2.14 Essen.
I'm wrote following script in QGis:

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

import processing
from processing.core.Processing import Processing

uri = QgsDataSourceURI()

uri.setConnection("localhost", "5432", "dbName", "user", "pass")

uri.setDataSource("public", "shapes", "geom")
layer1 = QgsVectorLayer(uri.uri(), "tbName", "postgres")

uri.setDataSource("public", "point", "geom")
layer2 = QgsVectorLayer(uri.uri(), "tbName", "postgres")

processing.runalg("qgis:joinattributesbylocation",{
                       "TARGET": layer2,
                        "JOIN": layer1,
                        "PREDICATE":u'intersects',
                        "SUMMARY":0,
                        "KEEP":0,
                       "OUTPUT":"path"})

vl = QgsVectorLayer("path", "name", "ogr")

saveUri = "dbname='dbname' host=localhost port=5432 user='user' password='pass' key=gid type=POINT table=\"public\".\"test\" (geom) sql="
crs_id = 4326
crs = QgsCoordinateReferenceSystem(crs_id, QgsCoordinateReferenceSystem.EpsgCrsId)
error = QgsVectorLayerImport.importLayer(vl, saveUri, "postgres", crs, False, False)

This script works well when run in QGIS, but i would like to run it from outside. When I am trying to write python script in PTVS, add:

imports

qgs = QgsApplication([], True)
qgs.setPrefixPath("path", True)
qgs.initQgis()

code..

QgsApplication.exitQgis()

First of all my processing.algList() contains only modeler: and grass: algorithms, second i can't connect to db. What should I do? Can I run the script from QGIS remotely?


texts: How to generalize shape with PyQGIS without much Distortion in the shape?

texts:
I was trying to generalize the the inner boundaries of the polygon. Looking for a PyQGIS solution for that any other possible python solution will be fine.

I tried convex hull and simplify geometry that available in qgispy but the results create distortion in the sape. For Convex hull i got overlapping polygon which is not suitable for analysis and for general simplification the shape distorted by a lot make the polygon unuseable.

So i seek help if anyone knows any way to simplify the inner boundaries and make those boundaries with less ups and downs,more like a straight line. any link or suggestion will do a lot for me.


texts: Direction of the orientation of shortest lines by Geometry Generator - Archaeological hachures symbol

texts:
I am creating, by Geometry Generator, an archaeological graphical style for lines. The goal is to be able to create automatically "hachures" that connect the external perimeter to the internal line (base of slope layer called "line" in my case).

I am working on creating a point layer at the vertices of the line using the Extract vertice tool. The points are symbolised with a complex style:

lines created with Geometry generator using the expression

shortest_line($geometry, aggregate('Lines', 'collect', $geometry))

a triangular simple marker at the beginning of the line

Both symbols are oriented using the expression

angle_at_vertex(shortest_line($geometry,  aggregate( 'Lines','collect', $geometry)),0)

The system works quite well with isolated geometries.

My problem is that when two polygons are too close, the shortest line representing the hachures  is created not in the direction I would like to, but connecting the external line to the closest line "break of slope" that in many cases is the wrong one.

Any suggestion on how to fix this error and how to be able to control the behaviour of the lines?


texts: Is QGIS working on new Apple computers, Silicon M1 Processor

texts: Does QGIS work on the new M1 Apple processors? Is it working natively or with Rosetta emulation?


texts: ArcInfo/ArcGIS Server Spatial Analyst not licensed

texts:
It might be me but there is something I’m not able to figure out.
I created a model in ArcGIS 10 that uses (among the others) “Extract by Mask” and “Zonal Statistics as Table” as functions (therefore Spatial Analyst extension) which runs smoothly on the desktop side. On this machine runs also ArcGIS Server (10) that I use to create a service that I want consume over the web.

Based on this thread, How to check the ArcGIS server type, I think I am running an Advance License because when I check the keycodes  file found in C:\Program Files (x86)\ESRI\License10.0\sysgen I find:

arcgisserver,100,[……],none, [……]

svradvanced,100, [……],none, [……]

arcsdeserver,100, [……],none, [……]

svrenterprise,100, [……],none, [……]

The server side should include the Spatial Analyst extension, right? (svradvanced) So, why do I get

ERROR 000824: The tool is not licensed.
FAILED

When I try to consume the service?

I probably should mention that when I check ArcGIS Administrator I have two permanent licenses for ArcInfo Desktop and one permanent for Spatial Analyst


texts: What is missing from this GeoTIFF to allow me to convert to EPSG:3857?

texts:
I have a source online for GeoTIFF files that I wish to convert to EPSG:3857 for Mapbox. (Flood Data for the United States.) A sample file is here :
(http://flash.ou.edu/flash_web/maxthres/maxthres.20200712.012000.tif).
When I run gdalinfo on the above file, I get this output :
Driver: GTiff/GeoTIFF
Files: maxthres.20200712.012000.tif
Size is 7000, 3500
ERROR 1: PROJ: proj_create_from_database: datum not found
ERROR 1: PROJ: proj_create_from_database: datum not found
Coordinate System is:
GEOGCRS["NAD83",
    DATUM["unnamed",
        ELLIPSOID["unretrievable - using WGS84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4269]]
Data axis to CRS axis mapping: 2,1
Origin = (-130.000000000000000,55.000000000000000)
Pixel Size = (0.010000000000000,-0.010000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=EF5 v1.0
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-130.0000000,  55.0000000) (130d 0' 0.00"W, 55d 0' 0.00"N)
Lower Left  (-130.0000000,  20.0000000) (130d 0' 0.00"W, 20d 0' 0.00"N)
Upper Right ( -60.0000000,  55.0000000) ( 60d 0' 0.00"W, 55d 0' 0.00"N)
Lower Right ( -60.0000000,  20.0000000) ( 60d 0' 0.00"W, 20d 0' 0.00"N)
Center      ( -95.0000000,  37.5000000) ( 95d 0' 0.00"W, 37d30' 0.00"N)
Band 1 Block=7000x20 Type=Float32, ColorInterp=Gray
  NoData Value=-9999

I was able to find a training manual showing how the development team that creates QGIS to create these geoTIFs, and it looked like it was in EPSG:4326. With that, I tried the command :
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 maxthres.20200712.012000.tif maxthres.20200712.012000-projected.tif

This is very close, however some of the pixels are off. (I may be wrong about 4326 too). The following output is also presented:
Creating output file that is 6553P x 4279L.
Processing maxthres.20200712.012000.tif [1/1] : 0ERROR 1: PROJ: proj_create_from_database: datum not found
ERROR 1: PROJ: proj_create_from_database: datum not found

Is there some command I am missing, or do I need more information to be able to convert this properly?


texts: Slow performance with Add Join and Calculate Field?

texts:
I have a point layer and a table with approximately 100 fields. I need to join the table to the point layer based on a key field and create a new shapefile / feature class.

This script has been painfully slow (40+ minutes for 1 layer w/100k records) to run. The run time for this part of the overall script makes it unstable due to memory issues. Is there a better way to do this? The frustrating thing is that the same function takes less than 2 minutes when performed through the ArcMap user interface.

    #Join CLS shp to off system table based on KEYFIELD1
    arcpy.AddMessage("running ELSE subscript for: "+fc_name)
    arcpy.FeatureClassToGeodatabase_conversion(CLS_file,output)
    arcpy.AddMessage("finish import cls to gdb")
    CLS_name = os.path.basename(CLS_file)
    CLS_filename = CLS_name.split('.')
    CLS_fc = output+"\\"+CLS_filename[0]
    CLS_fc_lyr = CLS_fc+"_lyr"
    input_dbf_view = input_dbf+"_view"
    arcpy.MakeFeatureLayer_management(CLS_fc,CLS_fc_lyr)
    arcpy.MakeTableView_management(input_dbf,input_dbf_view,"","","")
    arcpy.AddJoin_management(CLS_fc_lyr,"KEYFIELD_1",input_dbf_view,"KEYFIELD1","KEEP_ALL")
    input = fc_name+"_2013_join"
    arcpy.AddMessage("finish create view and join")

    #Report the % of records that failed to join + create new shp
    arcpy.CopyFeatures_management(CLS_fc_lyr,output+"\\"+input,"#","0","0","0")
    ##arcpy.TableSelect_analysis(output+input,output+input+"_NoCoords", "ADD SQL STATEMENT")
    ##error_cnt = int(arcpy.GetCount_management(output+fc_name+"_NoCoords").getOutput(0))
    ##total_cnt = int(arcpy.GetCount_management(output+fc_name).getOutput(0))
    ##print "The error % for "+fc_name+ " is "+str((float(error_cnt)/float(total_cnt))*100)+'%'
    arcpy.AddMessage("finish copy features")

    #Calculate Lat / Long Cooridnates
    arcpy.AddXY_management(output+"\\"+input)
    arcpy.CalculateField_management(output+"\\"+input,"LONGITUDE","!POINT_X!","PYTHON")
    arcpy.CalculateField_management(output+"\\"+input,"LATITUDE","!POINT_Y!","PYTHON")
    arcpy.DeleteField_management(output+"\\"+input,["POINT_X","POINT_Y"])
    arcpy.AddMessage("finish add lat/long")


texts: Python pyqgis how to load vector layer from a list

texts:
How can I load vector layer from a list? Such as

latitude = [...]
longitude = [...]
some_data_belongs_to_the_coordinates = [...]

vector_layer = QgsVectorLayer(?)


texts: Arcade Attribute Rules on joined feature layers (ArcGIS Pro)

texts:
I have a series of attribute rules defined.  The various fields are referenced using $feature.<fieldname>.  I know my rules are valid and well defined because they work fine as long as the feature class is edited on its own. I can successfully edit both as individual features via the attribute editor and with the CalculateField tool.
As soon as I join the feature class to another layer and attempt to use the CalculateField tool I get a 'Field not found' error.  I am guessing this is related to the fact that the attribute can't find $feature.<fieldname> in a joined table.  If I manually attempt to calculate a field on the same joined table the arcade syntax switches to $feature['<tablename>.<fieldname>'] - however this syntax does not appear to work when the feature class is edited on its own.
Is there a syntax that will work both with a standalone feature class and when participating in a join ?


texts:
I currently have X amount of binary files that are in grid format.

I have the span of these files in the form of latitude/longitude, plus the geographic projection. How can I turn these into georeferenced geoTIFFs/HDF/netCDFs?

For reference the latitude span is -60 to 60, and the longitude span is 85 to 205.

Data is projected in EPSG:4326.


texts: Open Source counterpart to Euclidean Allocation

texts:
I'm looking for a function that can perform something similiar to the Euclidean Allocation process in arcGIS (http://resources.arcgis.com/en/help/main/10.1/index.html#/Euclidean_Allocation/009z0000001m000000/) but using open source tools-- ideally Qgis or Gdal.

Alternatively, if there is a function that can take a polygon layer (e.g. map of the world) and run an algorithm similiar to Thiessen polygons-- that is, fill in the gaps in the extent using the nearest polygon that would be most helpful. For this latter solution, I'd prefer if I didn't have to convert boundaries to points.


texts: How to query "not contains" in ArcGIS Server REST API?

texts:
I want to query data that is not in a specified polygon in ArcGIS Server Rest API.
How can I achieve this like below image?


texts: Extending line by specified distance in ArcGIS for Desktop?

texts:
I have a purely aesthetic layer which has arrow symbols. Some aren't showing up right because the line is too small. I have selected maybe 50 records where I need to extend this line by a given number (ex. 2 meters). The extend line tool only extends lines to a specified intersection, so this tool isn't what I'm looking for.

I've tried editing the shape length field but it won't let me.
Is there a simple way of doing this via field calculator or within the editor tool bar?


texts:
I have a directory that contains around 700 .shp files. I want to create a map in R that displays all of them. The easiest way i could think to do this was to bind them all together and then use ggplot2 to create a map. here is my code to read the files and create a list containing the files:

dir <- "./Basins/separate_basins/"
basins <- list.files(dir, pattern = "\\.shp$")
basins <- unique(file_path_sans_ext(basins))
cord <- "+init=epsg:2163"
templist <- list()

for (i in basins){

  thisbasin <- readOGR(dsn=dir, layer=i)
  thisbasin <- spTransform(thisbasin, cord)
  templist[[i]] <- data.frame(thisbasin)
}

I then try to use the bind function to bind all of the components in the list together:
    test_bind <- bind(templist)
Then I tried to test if this works using
    `plot(test_bind)
but it prints like this:

I want to display the polygons that are present in the original .shp files. how can I fix this?


texts: Calculating shortest distance from point to line using QGIS?

texts:
What is the best method to calculate the distance from the nearest road for a series of address points in QGIS, having a very big dataset? (QGIS 2.0.1, Win XP)

I read the answers to previous similar questions:

Measuring the distance between lines and points in QGIS

Drawing perpendicular lines in PyQGIS?

Calculating minimum distance between points and polygons in QGIS?

https://stackoverflow.com/questions/16448402/distance-of-point-feature-to-nearest-polygon-in-r

Nevertheless:

MMQGIS plugin find the disance to the nearest node of a line network, not the shortest (perpendicular) distance to a line.
R and Python proposed scripts correctly works for small datasets and can be implemented in the modeler, however for my job I have 530 points and a regional street layer with about 135'000 features, and cycling through all the features then taking the shortest distance for all points will take an enormous amount of time
I have problems in making GRASS work in QGIS 2.0.1 (GRASS in QGIS not working (windows XP)).


texts: tool to find the routing distance between two latitude/longitude points using OSM data

texts:
I need a tool to find the distance between two points whose longitude and latitude are given. I am making an application so it will require bulk queries, so i cannot use google maps distance matrix API that could tell me the distance via road, airways or any other.
So I turned up to OSM data, i downloaded the OSM data and then pgrouting tool. But I am unable to find the distance via road between two points from this data stored in postgres table.

Any help regarding how to fetch this data from the postgresql db or any other resource you know that can accomplish my goal will be greatly appreciated.


texts: How may I use Umlaute programmatically to show up correctly in PyQGIS messagebar?

texts:
So far I couldn't figure out how to write something like u"Einstellungen aus ArcView übernehmen" to show up correctly as the titel of my messagebar. At least I can save my source (although this seems to depend heavily on the editor) by marking it as

-- coding: iso-8859-1 -- .

But it doesn't give me the desired result for my messagebar. I also tried .encode('latin1'), which helped for writing the string correctly to a file, but not for displaying it correctly in the messagebar.
Any ideas how to do it?


texts: Intersection of two polyline geometries not returning expected coordinates?

texts:
T2B = arcpy.Array([arcpy.Point(1,2),arcpy.Point(-1, -2)])
L2R = arcpy.Array([arcpy.Point(0, 2), arcpy.Point(1, -2)])
L2R_line = arcpy.Polyline(L2R)
T2B_line = arcpy.Polyline(T2B)
cent = L2R_line.intersect(T2B_line, 1)
arcpy.AddMessage("X co_ordinate "+str(cent.X))

I am trying to find intersection of two polylines.  Here is my sample code
But this is not printing x coordinate of intersection point.  I would like to know how to get x and y coordinates from "cent" object.


texts:
I am trying to connect my GeoNetwork to QGIS using the CSW protocol. The GetCapabilities request is OK.

But when I try to search some metadata this exception is raised:

Search error: 'Raised exception while searching metadata :
java.lang.IllegalArgumentException: Error when parsing spatial filter (version: 1.1.0):
<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
  <ogc:BBOX>
    <ogc:PropertyName>ows:BoundingBox</ogc:PropertyName>
    <gml311:Envelope xmlns:gml311="http://www.opengis.net/gml" srsName="urn:ogc:def:crs:OGC:1.3:CRS84">
      <gml311:lowerCorner>-23 -90</gml311:lowerCorner>
      <gml311:upperCorner>23 90</gml311:upperCorner>
    </gml311:Envelope>
  </ogc:BBOX>
</ogc:Filter>.

Error is: java.lang.RuntimeException: Parsing failed for Envelope:
java.lang.RuntimeException: Could not create crs: urn:ogc:def:crs:OGC:1.3:CRS84'


texts: Insufficient RAM to run QGIS - yet I have 32GB available

texts:
I have a question about running QGIS 3.8 Zanzibar on a MacPro running Sierra OS 10.12. I am hoping to use the app for species distribution mapping and was teaching myself using an online tutorial. However, the app kept failing at the same point saying I required 513 MB to continue. My Mac has 4x 8GB RAM chips installed, all functioning, but there does not seem to be a way to allocate more memory to any particular application.
Does this mean I will be unable to use QGIS or is there a way to get the application to obtain more of the available RAM (32GB). At this point I have done a restart and have very little else memory hungry running, but still the same issue. This is frustrating all my plans for a very large project.


texts: Is there a comfortable way to generate efficient raster catalogs in QGIS?

texts:
Is there a comfortable way to generate efficient raster catalogs (such as in Arcgis), where polygon tiles are loaded as long as the scale does not hit a threshold value and images are only loaded when the scale is large enough?

I tried to compute a virtual raster in QGIS, but it did not work for 200 tiles and gave me error messages reproduced by other users as well. Loading all the single tiles into QGIS is not an option since these are rendered (extremely slowly) on any action I do within QGIS. Is there a perfomant solution embedding large raster data into QGIS projects without sacrificing too much performance?

QGIS can´t read and correctly open Esri raster catalogs.


texts: Why am I getting "Project file error"? XAMPP and QGIS Server

texts:
I managed to get QGIS Server running locally on Windows10 using XAMPP, following this guide.

I clicked http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities as the guide suggests, and got "Layer(s) not valid", but I figured that the server was working, since I wasn't getting a "not found" error of some sorts.

I then wanted to load my own map, but the guide doesn't really explain how to do it. After searching around on the internet, I tried:
http://localhost/cgi-bin/qgis_mapserv.fcgi.exe?MAP=/OSGeo4W64/apps/qgis/bin/MAPNAME.qgs, but got a "Project file error".

What am I missing?


texts: Using Ripley's K-Function and Kernel Density Estimation?

texts:
Having issues applying Ripley's K-Function and Kernel Density to apply to ACLED (Armed Conflict Location & Event Data) point pattern analysis.

The issues relates to how do you know the correct distance increment for Ripley's and output cell size and search radius for kernel density based on the data you are faced with.


texts: Query ArcGIS Server Map Service with 3857 SRID?

texts:
I am trying to return the polygon a point is within using a latitude and longitude. I set the Input Spatial Reference to 4326 since I am using a latitude and longitude but I am still getting no results. Is there something I am missing?

https://wildfire.cr.usgs.gov/arcgis/rest/services/geomac_dyn/MapServer/23/query?&geometry=-118.031616%2C34.222874&geometryType=esriGeometryPoint&inSR=4326&returnGeometry=false&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=3857&returnZ=false&returnM=false&f=html

That is the url I am using.


texts:
A script I made creates domains within a selected geodatabase.

I need it to not create the domains if they already exist because if they do, it gives an error saying invalid domain name.

What's the best way to check if a domain exists?


texts: Where can I find the data that QGIS used to custom its menus as well as it's toolbars?

texts: I have noted that QGIS uses some data to restore its menuitems as well as its toolbar items everytime  when it is started.And I have being trying to hunt them out in my hard disks but sadly failed...The whole registry has been searched with the key word "OSgeo" in the "find" dialog,but nothing found to be related with what i want,neither under the whole directory of "OSGeo4W". File types of ".ini",".xml","*.cfg" have been used for search. That has made me sad,you all know. So,any kindness to help ?


texts: CartoDB choropleth + bubble visualization with two independent value ranges in legend

texts:
I have been trying to display two values simultaneously i.e. one absolute value displayed by bubble size (as with layer wizard: bubble) and one relative value displayed by color (as with layer wizard: chloropleth). (Adding to longitude and latitude these are four dimensions I want to visualize.)

See http://visage.co/data-visualization-101-bubble-charts/ for examples of bubble maps.

I have now managed to combine the CSS for a bubble layer (marker-width) and the CSS for a chloropleth layer (marker-fill). However I fail to adjust the legend to indicate both dimensions.

Is there a way to do this with CartoDB?

The thumbnail for Torque Cat in the Map layer wizward shows bubbles unequal in size and color, but I don't see how this could serve my purpose.


texts:
What is the GeoPandas equivalent of ST_MakePolygon?

I want to convert a MULTIPOLYGON to one single polygon. Currently, using GeoPandas.boundary a MultiLineString is returned.

How can this one be converted to a regular POLYGON?
As far as I understand ST_MakePolygon would be the post GIS pendant.

edit

I am aware of https://stackoverflow.com/questions/2964751/how-to-convert-a-geos-multilinestring-to-polygon-using-python, i.e. convex_hull, but this only returns a rather distorted geometry.

I.e. in case one looks at the shape of Austria this is:

edit2

A df.explode() will convert the MULTIPOLYGON to POLYGON.

Some answers suggest to aggregate via dissolve, however:

type(df(by='column_with_constant_value', aggfunc='sum').iloc[0].geometry)

still is

shapely.geometry.multipolygon.MultiPolygon

The same for:

type(df.cascaded_union)
shapely.geometry.multipolygon.MultiPolygon

edit3

The Polygon I am looking at is the Austrian borders from GADM. https://biogeo.ucdavis.edu/data/gadm3.6/gpkg/gadm36_AUT_gpkg.zip is the download link.

unzip it
df = gp.read_file()

The distorted image looks like:


texts: Authenticating using only token in ArcGIS API for JavaScript?

texts:
I am trying to embed a web map from ou ArcGIS Online on our corporate web app.

I implemented a server-side token generator, since our app users do not have (and should not have) an ArcGIS Online account, only the server.

However, I could not find a way to authenticate through the JS API using only a token or setting a custom getToken function.

I tried to simply add the token to the layer url, but it opened the login popup nonetheless.

const testLayer = new FeatureLayer({
        url: "https://services3.arcgis.com/.../ArcGIS/rest/services/.../FeatureServer/0?token=..."
      });

Ideally I also want the Javascript to able to access an ArcGIS Online item using only the token, for example a web map JSON.

I'm using ArcGIS API for Javascript 4.


texts: Is there a way to add a 'home extent' to a cartodb map?

texts: Within the zoom bar, is there a way to include a home extent button that the user could click to zoom to the same extent as the initial extent of the map?


texts: ArcGIS API for JavaScript 3.19 - export map to PDF layers and attributes

texts:
Is it possible in an ArcGIS for JavaScript app, to export a map to PDF with layers and attributes, just like in ArcMap?
In the printtemplates export format in only see

pdf | png32 | png8 | jpg | gif | eps | svg | svgz

Edit: this is how it works in ArcMap, when you choose Export Map in the File Menu:


texts: Turning off colliding keyboard shortcuts in ArcMap

texts:
While working in the table environment in ArcMap, I am not able to write several Polish characters without causing some action. Trying to write ą (AltGr+a) results in a Select All, ń (AltGr+n) ends up with flashing the select object, Ż (Shift+AltGr+Z) results with undoing last action, while writing Ś (Shift+AltGr+S) ends up with displaying Advance Table Sorting window.

How do I turn off those shortcuts?

I can't really find them in Customize Mode.


texts: Finding State, County, Place, and Census Block FIPS codes from a coordinate (latitude and longitude)

texts:
I am looking for a way to gather the political areas (ie FIPS codes) from a latitude and longitude coordinate.

For example, if I were to input the coordinates: 47.66955, -122.31259 (Downtown Seattle), I would expect to get the following FIPS codes:

  Place (aka "City"):

  name: Seattle

  FIPS: 53_63000

  County:

  name: King County

  FIPS: 53_033

  State:

  name: Washington

  FIPS: 53

Ideally, this could also include the census block, but at the very least, I need to return the State, County, and Place (aka "city") for a given coordinate.

I also understand that not all coordinates will map to a city/county/state, and that's ok.

Is there a library or machine image that I can deploy to solve this problem? Or perhaps a service (paid or otherwise) that I can use in production?

I have been researching this problem, and I found some interesting solutions. But none of them are working for me, as described below:

ESRI service at tigerweb.geo.census.gov
It looks like there was a similar question here, but the service in the answer to the question is down. If I'm going to use a service, I'd like to see an SLA (I wouldn't mind paying for a service), or have a way to deploy it myself.

Data Science Toolkit
Also, I found this project which solves my problem: http://www.datasciencetoolkit.org/developerdocs#coordinates2politics but it doesn't appear to be working and is not maintained. First, it's not a service - I have to deploy it myself to rely on it in production. And when I tried deploying it myself using their recommended AMI, the coordinates2politics endpoint was down, as described in this issue: https://github.com/petewarden/dstk/issues/66

data.fcc.gov
I also found this service to be helpful:
http://data.fcc.gov/api/block/2010/find?format=json&latitude=47.66955&longitude=-122.31259&showall=true

but it doesn't have an SLA ,and it doesn't return any Places

Any other suggestions or tips would be very helpful.

NOTE: SLA stands for Service Level Agreement: https://en.wikipedia.org/wiki/Service-level_agreement (ie "99.5% uptime")


texts: Analyzing values of rasters with different resolution

texts:
I have two rasters, on precipitation and population, and the precipitation raster is at a larger resolution.

I am interested in creating a table which tells me all of the values of the smaller resolution population rasters for each larger resolution precipitation raster. For instance:

I know how to do this with ArcGIS (using Tabulate Intersect), but I do not have access to ArcGIS at the moment.

I am looking to conduct this analysis for 15 countries for about 240 months.

Here is a visual example for Austria:

Here, the precipitation data is at larger resolution in the purples and blues, for one month (Jan 1991)
The population data is in the smaller resolution in the greens


texts: Removing JSESSION value from URL of GeoNetwork

texts:
Context
GeoNetwork 3 installed on a local machine running under Jetty 9.4 Server. HttpOnly option set to false
Steps to reproduce:

Open a browser (Chrome) and type: http://localhost:8080/geonetwork
The GeoNetwork's main page appears

Expected
The browser redirects to: http://localhost:8080/geonetwork/srv/spa/catalog.search
Reality
The browser redirects to the same site but concatenates the Session Cookie: http://localhost:8080/geonetwork/srv/spa/catalog.search;jsessionid=node0v85mmlv58rjgazg4acdgnxb30.node0#/home
Objective
I want to remove ;jsessionid=node0v85mmlv58rjgazg4acdgnxb30.node0 from the URL.
Any idea?


texts:
my code that takes one particular transformation method ("DHDN_to_ETRS") but it doesnt give you opportunity to choose it from drop down list looks like this:

import arcpy
import os

arcpy.env.workspace = arcpy.GetParameterAsText(0)
arcpy.env.overwriteOutput = True
outWorkspace = arcpy.GetParameterAsText(1)

try:
    for infc in arcpy.ListFeatureClasses():

        dsc = arcpy.Describe(infc)

        if dsc.spatialReference.Name == "Unknown":
            print('Koordinate System unbekannt!: ' + infc)
        else:
            outfc = os.path.join(outWorkspace, infc)
            outCS = arcpy.SpatialReference(25832)
            transform_method="DHDN_to_ETRS"
            arcpy.Project_management(infc, outfc, outCS, transform_method)
            print(arcpy.GetMessages())

except arcpy.ExecuteError:
    print(arcpy.GetMessages(2))

except Exception as ex:
    print(ex.args[0])


texts: How can I extract the parallel/perpendicular distance between two lines with QGIS?

texts:
I'm preparing a guide of my town for firemans and is very important to know the width of the streets but I don't know how to process this.

I have polygons that represents the block (islands) houses that simultaneously create the streeets.


texts: Converting .multiple shapefiles into Image Pyramid?

texts:
I have multiple shapefiles, and I want to import them as image tiles.

How do I convert a shapefile into image format and then convert that single image file into an image pyramid such as that of OpenStreetMap (OSM)?

I have worked with GeoServer for backend and OpenLayers to render the map but now I am asking about converting my own custom shapefiles into OSM tile pyramid.


texts: Python script to iterate through geodatabase and assign field value based on spatial location

texts:
I have two types of spatial data. One is a geodatabase of parcel data by county in Florida, and the other is a shapefile containing wind_zones through out Florida (used for enforcing certain building codes). I want to create and populate two fields in the parcel data based on which county and which wind_zone they are in. I have found a post here: Assigning field value based on location without join that answers a similar question and is what I am basing my code off of, but I need to iterate this process through an entire geodatabase that contains ~70 counties. Is there a line of code that can be added to iterate through the geodatabase so I do not have to manually run the code for each county?

Here is the code I have so far:

import arcpy

# Set overwrite option
arcpy.env.overwriteOutput = True

# KEEPING ORIGINAL LAYERS (NOT CREATING ADDITIONAL "JOINED" LAYER)
# Create FeatureLayers
arcpy.MakeFeatureLayer_management("C:/Data/HAAS/parc_join/parcels_1.gdb",    "lyr_parcels")
arcpy.MakeFeatureLayer_management("C:/wind_zones/All_WZs.shp", "lyr_windZones")

# Add a "LOCATED_WZ" & "WZ_CAT" field
arcpy.AddField_management("lyr_parcels", "Located_WZ", "TEXT", "", "", "20")
arcpy.AddField_management("lyr_parcels", "WZ_Cat", "SHORT")

# Create a search cursor for the WZs
rows = arcpy.SearchCursor("lyr_windZones")
for row in rows:
# What you'll do is select each WZ one at a time, and then select all the parcels in that WZ and calculate the "LOCATED_WZ" field
# NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
    arcpy.SelectLayerByAttribute_management("lyr_windZones", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
    arcpy.SelectLayerByLocation_management("lyr_Parcels", "INTERSECT", "lyr_windZones", "", "NEW_SELECTION")
    arcpy.CalculateField_management("lyr_parcels", "Located_WZ", "'{0}'".format(str(row.getValue("COUNTY"))), "PYTHON_9.3", "")
    print "Finished processing " + str(row.getValue("COUNTY"))

# Create a search cursor for the WZs
rows = arcpy.SearchCursor("lyr_windZones")
for row in rows:
# What you'll do is select each WZ one at a time, and then select all the parcels in that WZ and calculate the fields
# NOTE: If you are using not using shapefiles, then you'll have to change the FID in the line below to OBJECTID (or similar)
    arcpy.SelectLayerByAttribute_management("lyr_windZones", "NEW_SELECTION", "\"FID\" = " + str(row.getValue("FID")))
    arcpy.SelectLayerByLocation_management("lyr_Parcels", "INTERSECT", "lyr_windZones", "", "NEW_SELECTION")
    arcpy.CalculateField_management("lyr_parcels", "WZ_Cat", "'{0}'".format(str(row.getValue("DESCRIPT"))), "PYTHON_9.3", "")
    print "Finished processing " + str(row.getValue("DESCRIPT"))

Because each county is it's own feature class, I suppose I don't need to add a new field to tell me county name, because I will already know which county it is in.


texts: WMS source, selecting and printing region on multiple pages

texts:
I am new with GIS software.

I've imported topographic WMS source which shows whole country in ArcGis 10.2.2 and QGIS.

I want to select one specific region and print it in narrow scale (say, 1:15 000) on multiple pages (which would be something like 30 pages).

How to do that?


texts: Splitting Roads at intersection for ArcGIS Network Analyst?

texts:
I have a file geodatabase Feature Class representing the roads of a city.
The Roads are not split at Intersections. ArcGIS needs roads to be split at the intersections, otherwise it assumes that the roads are not connected.
This Image shows the Connectivity model for ArcGIS Network Analyst, when the roads are not split.

How do I correct this problem and split the roads wherever they intersect?


texts: Differences Between the Clip Tool and the Intersect Tool

texts: How is the clip tool different from the intersect tool in terms of how it works or the result? I don't have a specific situation in mind. I was learning about each of these tools and was curious about how they differ since the results seem to be very similar at first glance.  Thanks!


texts: What is the proj4string of the Lambert azimuthal equal-area projection?

texts:
Shapefiles of European data are available on the Eurostat website. I have loaded the NUTS_2013_60M_SH Shapefile in R with the command;

europeRG <- readOGR(dsn = patheurope, layer = "NUTS_RG_60M_2013", stringsAsFactors = FALSE)

I have managed to transform the map of Europe to the web mercator projection system used by  OpenStreetMap and Google maps (EPSG:3857).

europewmercator <- spTransform(europeRG, CRS("+init=epsg:3857"))
plot(europewmercator)

But a 2001 JRC and EuroGeographics document entitled Map projections for Europe recommends using the Lambert Azimuthal Equal Area projection:

  [...] Based on the above, it is documented that the projection to be used for statistical mapping in the EU, is the Lambert Azimuthal Equal Area projection.
  The projection center should be positioned at 9° E and 53° N. [...]

Wikipedia page of the Lambert azimuthal equal-area projection.
This site says the epsg code of the "Lambert Azimuthal Equal Area" projection is 9820.
Unfortunately EPSG code 9820 is not available on https://epsg.io/?q=9820 or http://spatialreference.org
but it is available on http://www.epsg-registry.org/
The issue is that the R function sptransform() doesn't recognise epsg:9820.

plot(europeRG)
# Transform To Lambert Azimuthal Equal Area projection
europeLAEA <- spTransform(europeRG, CRS("+init=epsg:9820"))

Error in spTransform(europeRG, CRS("+init=epsg:9820")) :
  error in evaluating the argument 'CRSobj' in selecting a method for function 'spTransform': Error in CRS("+init=epsg:9820") : no options found in 'init' file

What proj4string could I use to perform the conversion?


texts: Use field calculator in QGIS so get number of points in polygone and store it in attribute table

texts:
I have to layers, one point one polygon (with 6 areas) stored as shapefiles.
Now I want to count the number of points in each area and write it into the attribute table of the polygons. A column sum_poi already exists.

Can I use the field calculator to populate the column in the polygon layer?

Thanks for your help.


texts: Determining distance between two points along line feature in one direction using ArcGIS Desktop?

texts:
I have three  layers, two different point layers and a line feature. The line feature is representing a river system, and one set of points shows monitoring sites along the river (purple markers) and the other shows dams and impoundments (red circles).

I would like to calculate the downstream distance between a dam and the closest monitoring station. This is a small subset of a larger river system. Ideally I would be able automate this process for ~100 dam to monitoring station distances and get a distance for each downstream station, not only the closest one (i.e. a situation in which there are two monitoring stations that are downstream of a dam). The points do not exactly intersect the river system, and most are within ~1000 m.


texts: Not possible to add feature to Shapefile

texts:
While in editing mode of an already existing layer, the "add feature" button is not enabled. However, it is possible to add a circular feature, move feature, etc.

What can I do to have the button enabled again?


texts:
I'm new in pgrouting.
Is it possible to get 10 (some) ordered shortest paths ordered by distance (cost) with PgRouting (ex: SELECT TOP 10 SHORTEST)?

Update: why "second shortest" needed?

"second shortest" or "10 shortest path", this is the because the shortest path is not the best, for example if all people go to work with the shortest path, it will become a traffic jam, so "the shorted path is not the best". This is can be applied too with "10 cheapest path". Is this impossible?

-aris setyawan


texts: Restricting values to domain for Enterprise Geodatabase?

texts:
For an Enterprise GeoDatabase, I'd like to restrict the values of a feature class' field to a range of values. Essentially, I want validation at the database level, so that regardless of how the data is being inserted, the constraints are persisted (e.g. ArcGIS for Server feature service, ArcPy, ArcMap, ArcObjects, etc.). I tested domains, however I noticed that from the feature server service endpoint I can insert a value outside of my specified range.

For example, if I specify a range to 0 to 1, I can set the field to 6 from the service endpoint and it will be inserted into the GeoDatabase successfully.

I think database triggers or possibly creating a server object extension may be the only option.


texts: How to get Average of Points in Polygon?

texts: There are many number of polygon in my study area and in each polygon I have 16 points available.Is there any way to get multiple value of points to a single point.


texts: When editing a polygon in QField, the polygon is highlighted in red. This makes it almost opaque and difficult to use the fill ring tool if digitising over a feature shown on the base map below. Is there a solution to this problem? I've tried playing with the rubberband settings in QGIS but this doesn't seem to make a difference in QField.


texts: Testing correlation between two polygon features

texts:
I am working on a project looking at county boundaries for the state of California. I have 2 polygon datasets with this same geographic boundary, and these datasets, when mapped individually, visually appear to have a spatial correlation. Polygons are mapped according to numerical attribute values, which range widely on both datasets.

Wondering what the best method might be to test for a correlation between these datasets?


texts: Exporting raster from QGIS to use it in other programs

texts:
I am having a hard time trying to set a raster in a way that I can use it in AutoCAD map, for example. So these are the steps that I am doing:
Below is the image from the google satellite with the vectors of the lanes (in UTM Zone 23S, as is the project)

I then exported the map as image, and as you can see in the print below, when I put the output image it in QGIS, it landed in the right place.

Well, with the little I know about projections, I thought I had successfully exported the image so I could use it in AutoCAD map. Then I just went and exported the vectors as dxf, in UTM zone 23s as well, but when I tried to import both the raster and new dxf they weren't in the same place at all!
To better explain all my steps, I'm inserting here what I'm doing when importing the data to AutoCAD map.
First I'm converting the vector to dxf in QGIS in Sirgas 2000 UTM 23S

Then I'm opening the dxf file in AutoCAD map and assigning the coordinate system (Sirgas 2000 UTM 23S)

And lastly I'm importing the raster created when I imported as image in QGIS

Here the result of the journey, the vector in a place and the raster far away from it

Could someone help me to understand what did I do wrong and how could I achieve the results that I'm seeking?


texts: Select Features grayed out in PythonAddins toolbar

texts:
I've created a custom toolbar using the pythonaddins module. For ease of use I'd like to add the Select Features tool to the toolbar.

The tool becomes grayed out as soon as I close the Customize window however.

Note that the Select Features tool on the Tools toolbar works fine.


texts: Compile error in creating QGIS QT Designer UI file

texts:
Working on my first attempt with QGIS QT Designer I am trying to compile my UI file.   I am using Python27 found inside my QGIS installation (purposely not using QGIS v3 at this point).  I am using the following command:
C:\Users\username\.qgis2\python\plugins\projectfolder>"C:\QGIS\bin\python.exe"
"C:\QGIS\apps\Python27\lib\site-packages\PyQt4\uic\pyuic.py" -x ui_testproject.ui -o ui_testproject.py

However, I am getting the following error message in the command prompt:
Traceback (most recent call last):

  File "C:\QGIS\apps\Python27\lib\site-packages\PyQt4\uic\pyuic.py", line 31, in <module>
    from PyQt4 import QtCore
ImportError: DLL load failed: The specified module could not be found.

I don’t see a QtCore.dll file (is there even one?), or is it just looking for QtCore?
So I take it this may me an environment variable issue?  So I have added:
PYTHONHOME  C:\QGIS\apps\Python27

PYTHONPATH  C:\QGIS\apps\Python27\Lib\site-packages\PyQt4; C:\QGIS\apps\Python27\Scripts

Even added .SPC and .PYD as part of PATHEXT
What is missing?  What am I doing incorrect?


texts:
right now i am working in "MapInfo 10.0" tool where i stuck in a situation. I have a "GPS track" and "existing roads". Now i just want to select those roads where GPS track is present. The problem is there is no common attributes between these features.

My idea : GPS buffer and select features within the desired radius or distance.

Is there a way or tool to do such task or it would be done by coding?


texts: Orthogonal correction of building footprints after extracted from a raster

texts:
I have ArcGIS 10.1 and am working on a project involving extraction of vegetation and building classes from LIDAR (unclassified) data and aerial imagery (resolution 30x30cm).

I have managed to extract each class after following many LIDAR processing and Spatial Analysis steps; however, I could not find an automated way of correcting building footprints according to possible orthogonal boundaries. I have applied Simplify algorithms (in Cartography toolbox) in ArcGIS 10.1 and the best result as far as I can get is given in the first image. My question is:

Is there any free/open source library/script (preferably in Python) to correct polygon edges to (almost) true building boundaries? More specifically, I am looking for a method similar to the one given in the second image below:

from Zhang etal (2006). Automatic Construction of Building Footprints from Airborne LIDAR Data.


texts: Identifying inserts, updates, deletes, and conflicts between a parent and child version in ArcSDE with ArcPy

texts:
I'm looking for a way to identify records that have changed between the present and the last time a version was reconciled.  It would also be great to be able to identify conflicting records.  Ideally I would be able to read in the attributes of the records that have changed or are conflicting.

I suspect that this is possible doing some sort of comparison with the delta tables, though I have very little experience with versioned features and ArcSDE itself.  I've looked around and haven't found any information on how to actually read or access those delta tables.

Is there a way to achieve this within Arcpy?


texts: Routing difference between GraphHopper and ORS

texts:
I know that ORS (openrouteservice) is developed under GraphHopper routing library, but comparing those two heavy truck queries, I found the same distance but with a difference of two hours driving!

https://maps.openrouteservice.org/directions?n1=48.4184&n2=2.284057&n3=7&a=50.632786,3.059524,43.603902,1.440694&b=4a&c=0&k1=en-US&k2=km

https://graphhopper.com/maps/?point=50.635635%2C3.063008&point=43.605637%2C1.448693&locale=en-US&vehicle=truck&weighting=fastest&elevation=true&use_miles=false&layer=OpenStreetMap

I know also that the travel time is calculated for each segment by using speed-limits for different waytypes [see Travel Time Calculation] (https://github.com/GIScience/openrouteservice-docs).

If that's the reason why both queries output a different driving time is that the ORS determines different speed-limits to the ways.

Is there a way to change those speed limits before depolying my docker image so they can fit OSM?


texts:
We are developing a system to disseminate information about roadkills. The source data comes from three different databases. We have set up a FME project which pulls geojson from the source databases, joins them according too a set of business rules and then exports them to a front-end database.

We are currently exporting our data to a PostGIS database. FME does this very nicely. Now my colleague suggests we use Elasticsearch since this will make the presentation of our data in leaflet less complicated. Elasticsearch evidently indexes geojson on the fly.

We are using FME version 15.0. How can I write the data to Elasticsearch from FME? Which writer should I use?


texts: Digitizing boundaries for the educational office according to its affiliated schools

texts:
I have a layer of administrative boundaries, and within these borders, there are a group of separate schools for boys and girls, and these schools are under the supervision of several educational offices.
My goal is to digitize boundaries (Polygon) for the educational office according to its affiliated schools, taking into account the administrative boundaries of the city.
How can I digitize these boundaries correctly given the lack of time
using ArcMap 10.8 or ArcGIS Pro?


texts: How to create dataset in python with creation options like in ogr2ogr with -dsco options

texts:
In ogr2ogr there is dsco option (data source creation option).
e.g. ogr2ogr -f DXF  -dsco trailer=trailer_buildings.dxf -dsco header=header_buildings.dxf out.dxf test_buldings_shp.shp. It works well, but how can I use the same dsco options in Python?

from osgeo import ogr,gdal
import os,sys

def create_dxf(dxf_out_path):
    header_path = u'c:\header.dxf'
    trailer_path = u'c:\trailer.dxf'
    dxf_driver = ogr.GetDriverByName('DXF')
    dsco = {'HEADER':header_path,'TRAILER':trailer_path}
    dxf_ds = dxf_driver.CreateDataSource(dxf_out_path,**dsco)
    gdal.SetConfigOption('DXF_INLINE_BLOCKS', 'TRUE')

This code throws an error: TypeError: Driver_CreateDataSource() takes at most 3 arguments (4 given)

If I keep only one argument in dsco dictionary I get the following error:
TypeError: 'header' is an invalid keyword argument for this function. I've tried UPPER and lower cases but nothing helps. Documentation tells to watch supported options for each format here. And I clearly see that in DXF format there are HEADER and TRAILER creation option. But Python errors tells me vise versa things. I think option name in ogr2ogr utility doesn't correspond to option name in python gdal library, and if it so, where i can find documentation with full list of options for each vector format and for each language and utility? (e.g. ogr2ogr: -dsco HEADER='dxf_path.dxf'; python: kwargs = {'dxf_header':'dxf_path.dxf'}, etc...)


texts:
I cannot find examples of OL3 integration with EXT, such as done with GeoExt.  It looks like GeoExt is waiting for funding to do this. What I want is a feature grid with vector data displayed in a table/grid and two-way interaction between selected features in the map and the grid. This is done in may OL 2 and GeoExt examples, but I haven't seen anything like this with OpenLayers 3.

Can anyone point me to such examples/documentation?


texts: Creating Conditional Statement in QGIS Field Calculator

texts:
I've looked up some of the other Q&As on this topic, however, none have addressed what I am seeking to do.

I have a large data set and I would like to search one field for any values over a certain amount (in this case 20) and if the value is larger than 20, I would like to return a value of 1 in another field, otherwise return a value of 0.

I have tried writing something with no luck.

I am using QGIS 2.8.


texts:
I am trying to export a map canvas in QGIS I've made composed of .SHP files to a GeoTIFF, GeoPDF or any other georeferenced image.

All I want to be able to do is set a desired output resolution. I am unable to find anything on how do this. I keep getting pointed towards the Print Composer but that is for paper maps (DPI, etc.).

Where can I find this simple feature?


texts: QGIS Field Calculator update existing field using python console

texts:
How do I update existing field with a common number using Field Calculator Python console in QGIS 2.18.10?

I have tried:

processing.runalg('qgis:fieldcalculator', layername, "fieldname", 2, 20, 0, False, "3", path+"mynewlayer")

But my concern is each time I want to update a column I would have to create a new output layer. Is there a way of updating field without creating new or am I missing something?


texts:
I would like to know if there is a way to identify lakes based on an elevation raster (like SRTM or ASTER)? As far as I can understand, it would mean identifying any cell below the elevation of neighbouring cells, but I can't find anything related so I suppose it's not that simple...

The closest thing I figured out was to use r.watershed to identify small basins (10 to 100 threshold), and then rank them according to their accumulation. Some basins show a quite higher level than the surrounding ones, so I figured out it could mean a lake or something...

Any idea of how to do it better, in GRASS or QGIS? Is it even possible?


texts: Creating symbol with size and color based on two variables (one text, one numeric) in ArcMap?

texts:
I have a shapefile in which I need the symbology to be both size graduated and to represent another variable (say for example, "yes" or "no"), and another variable that is a numeric amount (say, from 1 to 10).

I'm attaching a diagram to illustrate what I want.


texts:
I'm looking for a package to develop a simple offline GIS desktop application that works on a tablet/laptop for using it on the terrain (namely for annotating existing maps).

The program has following requirements:

Able to show/load vector and image data offline (formats doesn't matter much)
Able to edit basic data. So able to save vector data.
Desktop or Android
We want to create our own routines for editing the data.

I have looked a bit to ArcGIS SDKs, but offline is in beta there so I was wondering if there are other (maybe open source) toolkits that some of you have experience with.

I'm not interested in things like buffering, geoprocessing,... We want to be able to show vector and image data and provide our interface on top of that.


texts: Calculating area of different land uses within boundaries

texts:
I am looking to calculate the percentage cover of different land classes in a number of districts in southern Malawi. The land cover map (raster) can be found here - http://www.masdap.mw/layers/geonode%3Amalawi_landcover_2010_schema_2, and the boundary data can be found here - http://geoportal.rcmrd.org/layers/servir%3Amalawi_adm2

I have used "tabulate area" in spatial analyst tools which seems to be what I need - I now have a table with values for each land cover type in every district. However, the values are very small (most are 0.0001 for example)... can anyone tell me what units this is in? My map is 30x30m resolution.

Also, once I have figured out what units the area of each land cover is, can anyone tell me how to find out the total area of each district? This will allow me to figure out % cover.


texts: Cut polygon shapefile with 2 parallel lines to derive new inside polygon?

texts: Using Python scripting: How can I clip a polygon shapefile with parallel lines (polyline shapefile) to derive new inside polygon?


texts:
Do any one knows online free tutorial in Introduction to GIS, where i can find the below topics.
1.Coordinate sysem
2.Georelational Vector Data model
3.Raster Model
4.Data Input
5.Geometric transformation
6.Spatial Data Editing


texts: Using angle attribute to split voronoi polygan to three small polygon in QGIS?

texts:
I have a point data set, that include [id] and [angle] attribute, like 0/120/240。I made voronoi layer use point, i want to use angle attribute split voronoi polygan to three small polygan, how to do in QGIS or GRASS ?

I try to calculate max distance of one point and associate polygan, using distance and angle draw line,  next find interaction point of line and polygan edge,  last using interaction point split voronoi polygan.

the way is true ?
How to do?


texts: How to export shapefile with colours from QGIS

texts: I want to export shapefiles from QGIS. The shapefiles have been colored by category, but when i export them, the exported files does not have colors, but only the lines that divide the categories. I want to mantain the colors while exporting instead.


texts: QGIS gdal_contour not respecting scale_factor/offset for netCDF

texts:
I am working with some netCDF files and want to import netCDF parameter's data as a  Raster and build a contour layer for it. I am using gdal_contour for this.

When I import the netCDF and choose a parameter (water_temp) in QGIS, the raster is loaded into the map with no problem and displays values in the range of roughly 4 degC to 31.25 degC.

However, when I use gdal_contour to make a contour layer for it, the values are in the range of -15944 to 11250. It certainly doesn't help that among other issues, it takes forever to generate the layer because I'm specifying an interval of 1.0 and the value range is far larger than the expected temperatures for Celsius.

From what I can tell, it looks like perhaps gdal_contour either isn't respecting the raster band's offset and scale_factor or has no knowledge of it. I understand that the netCDF is storing the temperature values as integers instead of floats to optimize file size, but I'm a bit confused by why QGIS can understand the offset when reading the netCDF into a raster layer, but not when generating a contour layer.

Am I missing something, or is there perhaps a caveat to using gdal_contour of which I'm unaware?

The command I am using to generate the conotur layer is:

gdal_contour -b 1 -a water_temp -i 1.0 -snodata -30000.0 -f "ESRI Shapefile" NETCDF:"C:/path/to/input/netcdf/INPUT.nc":water_temp C:/path/to/output/layer/OUTPUT.shp

The scale_factor, offset, and associated metadata for the band are:

add_offset=20
missing_value=-30000
NETCDF_VARNAME=water_temp
scale_factor=0.001
STATISTICS_MAXIMUM=11250
STATISTICS_MEAN=5475.011083141
STATISTICS_MINIMUM=-15944
STATISTICS_STDDEV=5863.9957197805
units=degC
_FillValue=-30000


texts:
I'm using QGIS 3.12 on Win10 x64.
I'm trying to configure QField for 1:n photos. I have a point layer representing trees, and I want to be able to add multiple photos to that pt.
I have set up as follow.
Layer ptTree has a unique field 'uuid'. Layer photo has a unique field also 'uuid' and a foreign key field 'fkid'.
I have created a relation between the two, as in the image below:

I have then configured the photo form as in the images below: first image shows the settings for the field 'uuid', second for the field 'fkid' and third for the field 'path'.

On this last image for the 'path' field, I think I have set 'relative paths' correctly, but I can't find any guidance on 'display resource path', or if you scroll down on the 'widget type' settings, there is also settings for 'Integrated Document Viewer' - this is hidden in the picture, but I have this set to 'no content' - I've tried 'image' but this does not seem to work.
In the form settings for the ptTree layer, I have the cardinality for the ptTree-photo widget set as 'Many-to-one' as in the image below.

In QField, when editing the ptTree layer, I get an add button to add a picture. I can then click the camera icon to bring up the camera. I can then take a picture, and press OK to accept and a record is created. But no path information is entered, and I cannot see a thumbnail of the picture. This was previously working on the same hardware (Galaxy Tab S4) so it's not a hardware issue.
Can anyone see what I may be doing wrong?


texts: Display tile grid borders with leaflet (visual debugging)

texts:
I'd like to have the tile grid borders displayed on a leaflet WMS Layer in order to see which areas are requested as separate maps. On top of that it should be possible to overlay a border to the tiles actually received from the server. That would allow checking server behavior at or near tile borders.

The only solution I've found so far is to enable server side watermarking on the WMS, which is influenced by tile generation parameters (e.g. meta tiling with GeoWebCache) and can thus be also used for comparison with the other border grids.

Is this feasible (even with a client other than leaflet) or are there any other workarounds to get the same functionality?


texts:
I am trying to put an image inside the geoJSON file.

I know, that Geojson works as string and I cannot input the image directly, however  I found some stuff, like this below:

http://www.gistechsolutions.com/leaflet/DEMO/basic/popupimage.htm

which could be working, when code would be written properly.
I believe, that problem lies somewhere in the quotes.

https://teamtreehouse.com/community/somebody-please-help-me-put-images-in-json

Anyhow my code looks like:

"type": "Feature",
             "properties": {
             "Title": "Stretham house",
                 "Head": "Flat to leave",
                 "Description": "<a href='https://en.wikipedia.org/wiki/Cambridge#/media/File:KingsCollegeChapelWest.jpg'File directory"
             "Pict":"image.png"
        },
            "geometry": {
            "type": "Point",
            "coordinates": [
            0.211212,
            52.344149
                ]
                  }
                    }];

And then...

var openreach = L.geoJSON(openpoor, {
pointToLayer: function (feature, latlng) {
feature.properties.myKey = feature.properties.Title + ', ' + feature.properties.Head
    return L.circleMarker(latlng, btMarkerOptions);
},
    onEachFeature: function (feature, layer) {
        layer.bindPopup('
'+feature.properties.Title+'

Address: '+feature.properties.Head+'

'+feature.properties.Description+'
'+feature.properties.URL+'
');
}

    }).addTo(map);

As you can see I tried to put an image by 2 ways. Link opens me the image URL, however it's not visible at the pop up.

Another option fell through due some code mismatch.

I would like to have an effect similar to Google MyMaps

I am aware, that GeoJSON does not support a way such this.


texts: qgis-processing: how to specify optional outputs?

texts:
Suppose a script has a variable number of outputs, which may or may not be produced depending on settings:

##Input_layer=vector
##Output_1=output vector
##Produce_output_2=boolean
##Output_2=output vector

from qgis.core import *

def create_my_output(path):
    # whatever...

create_my_output(Output_1)

if Produce_output_2:
    create_my_output(Output_2)

If this script is run with "Produce output 2" unchecked then QGIS displays a warning dialog saying Output 2 was not correctly produced, even though not producing Output 2 was clearly intended behaviour.

Is there a way to tell QGIS that an output is optional?


texts:
in a project I got some raster files with the following content (header and data are in one txt file):

HOUR 17
MINUTES 30
YEAR 2012
MONTH 5
DAY 4
CITY Emilia;Bologna
NROWS 99
NCOLS 156
ULXMAP 11.2087224106
ULYMAP 44.5539582648
XDIM 0.0016666668
YDIM 0.0012500001
CELLE_TRUE 2999
CELLE_FALSE 1
0.15478748  0.15478748  0.15478748  0.15478748  0.15478748
0.15478748  0.15478748  0.15478748  0.15478748  0.15478748
0.15478748  0.15478748  0.077375025 0.018038042 0.01576243
[...]

I tried to identify the format but did not get far. It seems to be a "kind of" ESRI ASCII grid, but it features non-square pixels (as you can see), which is not supported by ESRI ASCII grid (AFAIK). It might be an Arc Grid format (I found this link: http://landcover.usgs.gov/faq_importing_arc.php) but then I am lost. GDAL does not recognize the format, QGIS and ArcGIS 10 don't like the format either. Does anybody know this particular format (coming from an Italian source)? Or give me a hint how to process it to an more "readable" format?


texts: Converting lines to polylines with PostGIS/QGIS?

texts:
I have a PostGIS table containing linestrings.  Some of these are parts of larger logical constructs (tracks, roads or boundaries).  I have been going though the data using snapping tools to make sure that that coordinates match up properly.

I now want to change my "tracks" layer to have polyline features, one for each logical track so I can label them on the maps and export them to GPS as a single named entity, or convert them to a route.

I am pretty sure that the answer will involve ST_MakeLine but I cannot see how construct the query?

The line table has an attribute 'part_of'  which has the name of the logical entity and a sequence attribute too.

Pointers to a good basic tutorial on working with various geometric constructs would be useful.


texts: Assigning numbers to grid of polygons using ArcPy?

texts: I would like to number the shown grid of polygons.  The numbers should go left to right top to bottom, starting at 1.  I'm using ArcGIS Desktop with a Basic license, and would prefer to do this with ArcPy.


texts: Applying function to each element of SpatialLinesDataFrame

texts:
How do I apply functions to each element of a SpatialLinesDataFrame, in particular using other data also stored as an attribute (e.g. coordinates or elevation per line vertex)?
I've a number of SpatialLine objects which are grouped together in a single SpatialLinesDataFrame (SLDF) object with an "Id" name. Using the SLDF, I can extract elevations to the points along each line, and calculate things like the length of each line as in the example below.
library(raster)
library(sp)
library(rgeos)

# made up elevations somewhere in England
dtm <- raster(xmn=412115, xmx=439324, ymn=352217, ymx=391748)
dtm <- setValues(dtm,runif(ncell(dtm),min=1,max=150))
projection(dtm)<-CRS("+init=epsg:27700")

# create some lines as SLDFs
x1 <- spLines(rbind(c(420000,360000), c(420300,365000), c(420600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=1))
x2 <- spLines(rbind(c(425000,360000), c(425300,365000), c(430600,370000)),crs=("+init=epsg:27700"),attr=data.frame(id=2))
list_sldf=c(x1,x2) # in the main code I'm using, there are 100s of thesefirst added to a list
merged.lines <- do.call(rbind, list_sldf)  # combine to single SLDF

# extract elevation to each line
merged.lines$elev = raster::extract(dtm, merged.lines)

# calculate path length of each line
merged.lines$path_length=rgeos::gLength(merged.lines, byid=TRUE)

I now want to calculate additional metrics for each line, initally cumulative length per vertex and then eventually elevation specific derivatives e.g. along line curvature extract.
As an example, I've tried the following for calcualting cumulative distance (giving me an x axis to my "y" axis elevations):
If I try looping (just to start - ideally I want to vectorise this) the spatial dataframe to access each row, then I can't add new variables. For example, the function I have below to calculate the cumulative length between the vertices of each line is based on the coordinates of a line (e.g. cordinates[merged.lines[1,]) but I can't work out how to apply this sequentially across the SLDF.
#' Calcualte euclidean distance between 2 points
#' @example euc.dist(c(0,0),c(2,2))
euc.dist <- function(x1, x2){
  sqrt(sum((x1 - x2) ^ 2))
}

#' @param xys - a matrix of 2 columns (x/y)
#' @example
#' xys=merged.lines@lines[[1]]@Lines[[1]]@coords # a matrix
#' cumDist(xys)
cumulative.Dist<-function(xys,x='x',y='y', longlat=FALSE){

  colnames(xys)<-c(x,y)

  Dist <- 0
  for(i in 1:length(xys[,x])){
    if (i==1){
      print("first point")
      Dist[i]=0
      }
    else{
      print("other points")
      Dist[i] = euc.dist(c(xys[i,x],xys[i,y]), c(xys[i-1,x],xys[i-1,y]))
    }
  }

  DistCumulative <- 0
  for(i in 2:length(xys[,x])) {
    DistCumulative[i] = Dist[i] + DistCumulative[i-1]
  }

  #xys_cumulativeDist <- cbind(xys, Dist, DistCumulative)
  return(DistCumulative)
}

This works to get the values on an individual basis:
xy_df_=as.data.frame(coordinates(merged.lines[1,]))
cumulative.Dist(xy_df_)

How do I get it to run over each line in the SLDF?
Can I vectorise this?
for (i in 1:length(merged.lines)){
  merged.lines[i,]$cumulative_length=cumulative.Dist(as.data.frame(coordinates(merged.lines[i,])), x='X1', y='X2')
}

Perhaps I am setting the variable incorrectly.


texts: Find Feature class into another feature data

texts:
I am trying to find the feature classes from one feature data into another feature classes. But I'm getting nothing at the end as a result of the matc statement. The data in the arc environment is XYZ.LandCover while in the other environment data naming format is pqr_landcover

One environment is SDE and other is postgis
Below is a script I'm using

import arcpy
import os, subprocess, fnmatch, string, re

arc = ("Connections\\arc.sde")

post = ("Connections\\Post.sde")

arcpy.env.workspace = arc

listfeatures = arcpy.ListDatasets()
for lf in listfeatures:
    print lf

    fclist = arcpy.ListFeatureClasses("","",lf)
    for fc in fclist:
         print fc
    arcpy.env.workspace = post
    postgisdata = arcpy.ListDatasets()
    for postg in postgisdata:
      print(postg)
    postlist = arcpy.ListFeatureClasses("", "",postg)
    for postfc in postlist:
        print postfc

        if postfc in fc:
                matches = matches+1
        if matches > 0:

            print str(matches)+" match found for "+postfc

        else:
            print "No matches found for "+postfc