texts: Proper use of Modulus "%" operator in ArcMap

texts:
Objective:
Gather pour points to define a specific-sized watershed. In this case, 1000 pixels. (for this thought exercise, watershed is to mean any point that flows through a pour point, regardless of size, not the NHD definition)
Layers:
DEM, flow direction, flow accumulation, the like.
Flow accumulation symbolized with 1/2 standard deviations on top of the DEM:
My mathematician brother suggested the use of the modulus operator, which returns the remainder of a division problem, thus allowing something like `Con(("flow_accum" % 1000) < 10, 1, 0) in the Raster Calculator.
Theoretically, this should return a blip in the stream every thousand pixels of flow accumulation, but it doesn't.
What just the modulus ( just "flow_accum" % 1000) returns is this:

I do not know how to interpret this. It should have low values for all pixels that are not near a X,000 but it appears to have just done some sort of rescaling to 1000 or something.
I did a manual selection of the first 12 watersheds with:
Con(("flow_accum" == 1000), 1, Con(("flow_accum" == 2000),2, Con(("flow_accum" == 3000)
,3,Con(("flow_accum" == 4000),4,Con(("flow_accum" == 5000),5,Con(("flow_accum" == 6000)
,6,Con(("flow_accum" == 7000),7,Con(("flow_accum" == 8000),8,Con(("flow_accum" == 9000)
,9,Con(("flow_accum" == 10000),10,Con(("flow_accum" == 11000),11,Con(("flow_accum" == 12000)
,12,Con(("flow_accum" == 13000),13,Con(("flow_accum" == 14000),14,0))))))))))))))

Which after running the watershed tool, returned this:

What I need help with is making sense of the output of the modulus function, and defining arbitrarily-sized watersheds.


texts: Saving borough boundaries as MultiLineStrings or MultiPolygons?

texts:
I am creating a database in PostGIS and I have a column that will save the boundaries of boroughs.

The app will later query the database, asking what points, lines or polygons are inside those boundaries.

The borough boundaries are closed linestrings, but I dont know if this affects the query speed or the efficiency of functions like ST_Within or ST_Intersects.

In that context, what is the correct way to save boundaries , as closed linestrings or as polygons ?


texts:
I am brand new to GRASS GIS (and GIS in general) and I am trying to use this for some real estate research.  However, before even going down that path, I'm attempting to follow the Quick wxGUI tutorial (https://grasswiki.osgeo.org/wiki/Quick_wxGUI_tutorial) just to see if I can open a map.  I am having difficultly even doing this.  Here is what I'm attempting to do:

I have installed Grass GIS 7.4.0 in Mint Linux
I downloaded the 'GRASS 7 full data' in zip format (nc_spm, 145MB)
I have created a GIS Data Directory i '/home/me/grassdata'
in the 'grassdata' folder, I have 3 directories (landsat, PERMANENT, user1) and 3 files (CREDITS.txt, HISTORY.txt, VERSION.txt)
When I open GRASS GIS, I navigate to the 'grassdata' folder, it indicates 'No GRASS location found in /home/me/grassdata'. Create a new location or choose different GRASS database directory'

Which directory do I need to browse to for the 'Grass GIS database directory'?

Also, I don't see 'nc_spm_08' or 'spearfish60'.


texts: Using only Selected Layer in mMapLayerComboBox

texts:
I am building a plugin and so far I have tried simplifying the tasks. But my greatest challenge is making sure that the selected layer in the "mMapLayerComboBox" is what is exactly used in the operation shown below:
def run(self):
    mains = self.dlg.mMapLayerComboBox.currentLayer()       ### vector layer
    laterals = self.dlg.mMapLayerComboBox_2.currentLayer()  ### vector layer
    ovalay = self.dlg.mMapLayerComboBox_3.currentLayer()    ### Raster layer

    for layer in checkedLayers:
        if layer.type() == QgsMapLayer.VectorLayer and self.dlg.checkBox_3.isChecked():
           # get extension about the raster
           filename, file_extension = os.path.splitext(layer.source())
           output1 = self.folderName + "/Profiles/m-length_" + layer.name() + ".shp"
           processing.run("qgis:exportaddgeometrycolumns", {"INPUT": mains, "CALC_METHOD": 0, "OUTPUT": output1})

In the "INPUT", I had specified the layer to use. However, while the task works, the main problem is that the task executes for all loaded vector layers in QGIS, even if I have a basemap like from google satellite it just tries to run the task on all available vector layers instead of the one I specified as:
              mains = self.dlg.mMapLayerComboBox.currentLayer()       ### vector layer

What I get as an output is a result for all the loaded layers including the specified layer in the "mMapLayerComboBox". I just want one output and not for all loaded vector or raster layers as the case may be.
How do I resolve this problem when using "mMapLayerComboBox" from QGIS built-in custom widget?
I want my task to only perform operations on the selected layer, that's all.


texts: combine raster with another raster mask file - gdal

texts:
I have 2 rasters in GeoTiff format. One has all of the land and water and the other is a mask file that shows either 0 for water or 1 for land. I basically want my output to just have the land values and all the areas that were water in the first image to just be 0,0,0. Is there a way to do this with any of the gdal command line tools?

Basically I want to do what gdal_rasterize does with a shp file mask but with a raster instead.


texts: Photogrammetry route flown using Litchi misses photos at some waypoints

texts:
I have recently started using Litchi to fly waypoint-based photogrammetry flights with a DJI Phantom 4 Pro. The routes were created in DJI Flight Planner, then edited, uploaded, and executed in the field using Litchi. The photos collected on these flights will be used to create aerial photo mosaics and high-resolution DEMs.

However, I’ve encountered an unexpected problem. Every flight so far has a number of waypoints where photos were not collected by the drone. The location of these missed photos varies between routes, and appears to be random. For example, in the photogrammetry route below, a photo was not collected at Waypoint 39 the first flight. The second time I flew the exact same route, photos were not collected at Waypoints 5, 6, and 14.

A different route with 99 waypoints in another location ended up with 9 missing photos. Other users have encountered similar issues, but I can’t seem to find any suggestions on what might fix the problem.

Does Litchi (or the Phantom 4 itself) have some sort of built in tolerance or distance threshold beyond which it will not take a photo? In other words, if the drone is too far from a waypoint once it reaches that point in the route, will it skip that photo? For example, if the drone reaches a waypoint and then a gust of wind blows it to one side could that cause it to skip a photo? This is just a guess, as the 99-waypoint route that missed more photos was flown on a windier day.


texts: How to create list from spatial join table using ArcGIS Desktop and ArcPy?

texts:
I have performed a spatial join on some data and now I would like to take that data and create a list. I'm trying to do this dynamically so that it can be easily repeatable via a script or model.

I have two fields in my data that I care about: Station Name and Page Number. The overall purpose of my project is to take a 632 page mapbook and divide it into smaller books based on the Station Name Boundary. Exporting to pdf only the pages that are relevant to each station.

My spatial join table has given me One_To_Many list of each page and what station boundary it falls within.

Now I'd like to create a list that is something like:

Station1
Pages: 1,2,3,4,5,6...

Station2
Pages: 8,9,10,11,12...

and so on.

For what its worth I'm using 9.3.1 and the maplogic extension for the mapbooks.


texts: GeoServer GeoJSON WMS layer request in Mapbox

texts:
I'm creating a Mapbox GL JS web map which pulls GeoJSON layers from GeoServer. I've used the documentation to find my problem lies in the 'GetFeatureInfo' get request URL which follows the following format (though in one line of course):

http://localhost:8080/geoserver/wms?
INFO_FORMAT=application/json&
REQUEST=GetFeatureInfo&
EXCEPTIONS=application/vnd.ogc.se_xml&
SERVICE=WMS&
VERSION=1.1.1&
WIDTH=970&
HEIGHT=485&
X=486&
Y=165&
BBOX=-180,-90,180,90&
LAYERS=CampusMap:dwt-campus-json&
QUERY_LAYERS=CampusMap:dwt-campus-json&
TYPENAME=CampusMap:dwt-campus-json&

My problem is that this link does not work, and when I run it in a browser I receive various WMS error files that are not very helpful. I really just want to be able to get this link to work so I can have this building footprints file load on my basemap and be query-able for future use.


texts:
I created a Shapefile from Google Earth (kmz) consisting of a set of points which I had used as a layer for Heatmap analysis.

I opened the kmz file in QGIS and utilized the heatmap plug-in and it has 8 layer units as the default radius.

Can anybody enlighten me what does this 8 layer units mean?


texts: Reversed values from slope analysis of DEM in QGIS?

texts:
I have a DEM as a Geotiff and I have conducted slope analysis in QGIS. I've done this many times and everything was ok but this time I received reversed values: 0 for max slope and 89 for min slope. I have no idea what is wrong. Geotiff was generated from .asc files and seems to be ok.

I suppose there is something wrong with the first Geotiff file.
I took "native" .asc file and made slope analysis without a problem.
After that I generated Geotiff from the Slope.asc and is also fine.


texts: Keeping alias field names on CSV export using ArcPy

texts:
I'm exporting my feature class to a CSV file. I figured how to keep the fields I need on export.
Is there a way to keep the alias field name on export?
I'm also open to writing out the correct headers but can't figure out how.
#Export out as CSV
out_folder_path = arcpy.GetParameterAsText(0)
outputCSV = os.path.join(out_folder_path, "Points" + ".csv")
fields = arcpy.ListFields(fc)
fieldsToKeep =['REPORTER','RECOWNER','SCINAME','COMNAME','RECBASIS','GlobalID','OCCSTATUS','MANAGESTAT',
               'POPSTAT','OBSDATE','COUNTRY','STATE','COUNTY','PRIVATE','LATITUDE','LONGITUDE',
               'ABUNDANCE','INFESTAREA','INFESTUNIT','PCTCOVER','VISITTYPE','COMMENTS']

field_names = [field.name for field in fields if field.name in fieldsToKeep]

with open(outputCSV,'w', newline='') as f:
    w = csv.writer(f)
    w.writerow(field_names)
    for row in arcpy.SearchCursor(fc):
        field_vals = [row.getValue(field.name) for field in fields if field.name in fieldsToKeep]
        w.writerow(field_vals)
    del row


texts: How can I add a second Toolbar Button and Dialog to a QGIS plugin built with Plugin Builder?

texts:
I want to create a toolbar with 2 buttons, each button opening a different form.

The steps I followed are:

Created a plugin using Plugin builder
Created a toolbar with 2 buttons by doing the necessary modifications in the implementation file.

Until here, everything works fine, so when first button is clicked, the corresponding form opens correctly. This form is the one created by the Plugin builder.

Then I want to have a second form which will open when the second button is clicked.
So:

Created a new form in Qt designer (form2.ui),
Translated the ui file to py file (make) and imported the class from the translated py file into the implementation file, got an instance of the class (self.dlg2 = Ui_form2())

and in the run() method of the implementation file I wrote:

def run2(self):
    self.dlg2.show()
    result = self.dlg2.exec_()

But then I get an error message that this new form has no .show() method. Should it inherit something in order to implement the .show() method?
But how can I modify the code of the second form in order to inherit something since all changes in the translated py file from the ui file cannot be saved?


texts: Identifying the cause of NOT ST_IsValid()

texts:
I have a table of raw GeoJSON in a Postgres database, with PostGIS installed.  Some of this GeoJSON is unfortunately invalid, for varying reasons.

SELECT
    id, name, geometry
FROM shapes
WHERE
    NOT ST_IsValid(ST_GeomFromGeoJSON(geometry));

Sample log output:

NOTICE:  Self-intersection at or near point -98.12763909659509 40.524842899558408
NOTICE:  Ring Self-intersection at or near point -89.500197172165002 45.280231548720003
NOTICE:  Nested shells at or near point -80.15264234 34.224868579999999
NOTICE:  Duplicate Rings at or near point -80.155240000000006 34.217115999999997
NOTICE:  IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4
NOTICE:  IllegalArgumentException: point array must contain 0 or >1 elements

Self-intersection is by far the most common error.  I'd like to edit out the self-intersection completely, such as via ST_SimplifyPreserveTopology.  However, when I attempt to do so:

SELECT
    id,
    name,
    geometry,
    ST_AsGeoJSON(ST_SimplifyPreserveTopology(ST_GeomFromGeoJSON(geometry), .00000001)) as new_geometry
FROM shapes
WHERE
    NOT ST_IsValid(ST_GeomFromGeoJSON(geometry));

Those IllegalArgumentExceptions kill the query:

ERROR:  First argument geometry could not be converted to GEOS: IllegalArgumentException: Invalid number of points in LinearRing found 3 - must be 0 or >= 4

I tried using ST_IsSimple to filter those geometries out, but that function also triggers an IllegalArgumentException.

How can I more accurately identify which validity error a GeoJSON string suffers from so that I can filter out those nastier geometries?


texts:
I wrote a specific function to calculte metrics on a catalog but this function (myMetrics) work with a LAS file but don't work with a catalog?

myMetrics = function(las, ...)
{
  UseMethod("myMetrics", las)
}

myMetrics.LAS <- function(las){
  metrics <- list(
    min_x <- min(las@data$X),
    max_x <- max(las@data$X),
    min_y <- min(las@data$Y),
    max_y <- max(las@data$Y),
    min_z <- min(las@data$Z),
    max_z <- max(las@data$Z),
    min_classification <- min(las@data$Classification),
    max_classification <- max(las@data$Classification)
  )
  names(metrics) <- c(
    "min_x", "max_x", "min_y", "max_y", "min_z", "max_z", "min_classification", "max_classification"
  )
  return(metrics)
}

myMetrics.LAScluster <- function(cluster)
{
  las <- readLAS(cluster)
  if (is.empty(las)) return(NULL)

  metrics <- lasmetrics(las, ~myMetrics(las))
  return(metrics)
}

opt_chunk_buffer(ctg) <- 0
opt_chunk_size(ctg) <- 0
opt_output_files(ctg) <- ""
plots_metrics <- catalog_apply(ctg, lasmetrics, func = .myMetrics(ctg))

  Error in catalog_apply(las, lasmetrics, func = .myMetrics) :
    object '.myMetrics' not found


texts: How to parse and display KML on Bing maps under SSL

texts:
This may seem like a simple question at first, but let me explain.  We have several mapping solutions, however one is strictly Bing maps (v6) based.  "Use OpenLayers" will not be a possible solution for this question.

Our product runs behind a login and SSL and we use the HTTPS hosted version of Bing Maps. However... Bing maps under SSL does NOT support adding a KML file.  It does under HTTP. The way it works is a URL is passed in, then the MS Server downloads the file, parses into Bing shapes sends back to client and the client renders it.  This is all automated by just adding the KML URL. However this is not supported under SSL. There's a 3-year old post on MS forums about this already.  Also, running Bing maps on HTTP on our page will not work as it generates the wonderful "mixed content" warning and is generally not good practice.  So that's not a solution either.

So I've been thinking about consuming the KML file through our own server-side proxy (since URL is from a different domain) and parsing it manually into Bing-friendly shapes.

Any ideas about parsing KML and displaying on Bing maps?  Essentially what OpenLayers (OL) does already, but it's all in an OL format.  Is it possible to reuse some of that KML code?  Or do you know of another library that would allow us to parse KML and display on Bing maps v6?

Here's a link to the MS Forum discussion where they basically acknowledge the problem under SSL


texts: Splitting self-overlapping lines with points using PostGIS

texts:
I have a line geometry that partially overlaps itself (think about a line going A -> B -> C -> B -> D) and a set of points (not on the line but close to it) that I want to use to split the line in segments.

I tried to project the points on the line but all the solutions here involving ST_ClosestPoint and ST_Snap don't work, so I created new lines by "extending" the lines between the points and ClosestPoints (considering sin and cos for correct heading) and collected these lines to use them in ST_Split.

When I use ST_Split(line, geom_collection), the line is correctly cut into segments based on the points but only where it doesn't self-overlap; where it does self-overlap, it is splitted into segments at every single vertex, not only at the points of my geometry_collection.

Do you have any idea?


texts: Complicated expression for Calculate Field for arcpy

texts:
I have a working formula that I can put in the Field Calculator which I'd like to use in my python code:

  .3-.001*Log ( [Field1]+1  )-.02*Log ( [Field2]+1 )-.04*Log (
  [Field3]+1 )

My code:

import math
expression =".3-.001*Log(!Field1!+1)-.02*Log(!Field2!+1)-.04*Log(!Field3!+1)"
arcpy.CalculateField_management("Layer","Field4",expression)

  ExecuteError: ERROR 999999: Error executing function. Syntax error
  Failed to execute (CalculateField).

How can I write this expression?


texts: OpenLayers3 vs Leaflet for styling multifeature geojson

texts:
I'm experimenting with various options for displaying geographic information to the public in html format. The ultimate aim is to build an online application where anyone can add the route they cycle to work and points of danger, for a project comparing perceived risk with actual accident data.

Anyway, getting to the point, some colleagues and I have added our routes to a geojson file hosted by GitHub for a 'proof of concept' map. We used the nifty geojson.io interface to create this, and it can be seen on our GitHub page.

The question, now that we have some example data online is how best to visualise it, styling the points for different levels of severity and risk. A colleague has demonstrated how the points can be styled based on attributes, but the code is very long compared with what can be done quickly in CartoCSS:

 function initialize()  {
    var map = L.map('map').setView([53.851,-1.61], 12);
    L.tileLayer('http://tile.openstreetmap.org/{z}/{x}/{y}.png', {
        attribution: 'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © CloudMade',
        maxZoom: 18
    }).addTo(map);

    L.geoJson(geojsonFeature, {
        onEachFeature: onEachFeature,
        pointToLayer: pointToLayer,
        style : myStyle,
    }).addTo(map);
}

function onEachFeature(feature, layer) {
    var popupContent = feature.properties.Person + " " + feature.properties.Place;
    layer.bindPopup(popupContent);
}

function pointToLayer(feature, latlng) {
    return L.circleMarker(latlng, myStyle);
}

function myStyle(feature, layer)    {

    switch (feature.properties.Place)   {

        case 'Home' :
            return style4;
        break;

        case 'Work' :
            return style4;
        break;

        case 'Danger':
            if (feature.properties.Severity == 1)   {
                return style1;
            }

            else if (feature.properties.Severity == 2)  {
                return style2;
            }

            else return style3;

        }
}

var style1 = {
    "color": "#00FF00",
    "radius" : 8,
    "weight": 1,
    "opacity": 1,
    "fillOpacity": 0.5
};

var style2 = {
    "color": "#FFFF00",
    "radius" : 8,
    "weight": 1,
    "opacity": 1,
    "fillOpacity": 0.5
};

var style3 = {
    "color": "#FF0000",
    "radius" : 8,
    "weight": 1,
    "opacity": 1,
    "fillOpacity": 0.5
};

var style4 = {
    "color": "#0000FF",
    "radius" : 8,
    "weight": 1,
    "opacity": 1,
    "fillOpacity": 0.5
};

I have looked at the geojson documentation for both Leaflet and OL3, but cannot find out how to do this properly. So really I have two questions, one specific and one more general:

Would OpenLayers3 be able to do this better?
Are there any other options for quickly styling geojson data for webmaps?


texts: Couldn't load SIP module - python disabled

texts:
Recently my QGIS started refusing to load SIP module, from day to day without any significant changes in the system (I'm on Windows7 x64). Message thrown at startup:

Nie można wczytać modułu SIP. (Couldn't load SIP module)
Wsparcie dla Pythona zostanie wyłączone. (Python support will be disabled)

Traceback (most recent call last):
  File "", line 1, in
ImportError: cannot import name wrapinstance

Wersja Pythona: (Python version)
2.7.3 (default, Apr 10 2012, 23:31:26) [MSC v.1500 32 bit (Intel)]

Wersja QGIS: (QGIS version)
1.8.0-Lisboa 'Lisboa', f440b60

Lokalizacja Pythona: (Python location)
['C:/PROGRAMS/Geodesy/OSGeo4Win/apps/qgis/./python', 'C:/Users/Muciek/.qgis//python', 'C:/Users/Muciek/.qgis//python/plugins', 'C:/PROGRAMS/Geodesy/OSGeo4Win/apps/qgis/./python/plugins', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\qgis\\python', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib', 'C:\\Windows\\system32\\python27.zip', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\DLLs', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\plat-win', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\lib-tk', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\qgis\\bin', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages\\PIL', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages\\win32', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages\\win32\\lib', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages\\Pythonwin', 'C:\\PROGRAMS\\Geodesy\\OSGeo4Win\\apps\\python27\\lib\\site-packages\\wx-2.8-msw-unicode']

I tried reinstalling various versions of QGIS (1.7 - 1.9) and manipulated with system Path and PYTHONPATH both with no result.

QGIS starts perfectly on different user account (without administration rights). It means that it may be a problem with sone software installed but still I have no idea what could it be.

Any help or suggestion would be appreciated, thanks in advance.


texts: Group geometry points and create a linestring for each group

texts:
I have a table of points at a Postgres database that have two main columns (geom geometry, time timestamp).
I need a query to retrieve as linestrings or as array of the geometries the points that are close to each other by 1 m aprox and follow the sequence of time.
The points table:

Example:

Within this set of points, how can i get in a query those two sequences of points as two different linestrings? Also excluding points that are isolated.
Assuming that this is a particular case within millions like those two cases, which would be the more efficient way to find all the linestrings?


texts: Merging hundreds of Aerial Photos into one big file?

texts:
I have tried to merge my 405 Aerial Photos into one. I have FME and ArcInfo 10.0 at my disposal on Windows 7 64bit.

I have tried the following: create a mosaic; export all images from a rastercatalog; create new empty raster and then loaded the rasters, and in FME I have tried the RasterMosaiker.  Everything starts ok, but then due to the file size creates an error after many hours.

Does anyone have any experience creating such large single raster data files and can offer some tips on how I could achieve this,

My only restriction is that I don´t wan´t to decrease resolution.


texts: Wireline internet (cable, dsl) availability data source

texts: I'm looking for data on spatial availability of cable, dsl, etc. internet service in Washington State, USA.  Mainly interested in the big telco providers Comcast, ATT, etc.


texts:
I have calculated curvature using ArcGIS (curvature) and SAGA (Slope, Aspect, Curvature; based on 9 parameter 2nd order polynom (Zevenbergen & Thorne 1987) method), however, the output rasters are not the same.

I've done some research and now I understand how curvature works in ArcGIS:
http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/how-curvature-works.htm

But it's still a mystery to me how SAGA calculates curvature.

Which formula is used?


texts: Google Maps API: Can't load GeoJSON when inside a function

texts:
I use the sample provided on the Google Maps API page (https://developers.google.com/maps/documentation/javascript/datalayer).

I put the sample inside a function like following

function loadCableGeom(){
        _this.map.data.loadGeoJson('https://storage.googleapis.com/mapsdevsite/json/google.json');
        _this.map.data.setStyle({
          fillColor: 'green',
          strokeWeight: 1
        });
        console.log("LOAD");
}

function removeCableGeom(){
        _this.map.data.setMap(null);
}

I am sure the call to loadCableGeom() is done as the console.log("LOAD") works well.

If the loadGeoJson is put outside the function just above the function definition (in this case I remove the underscore), the GeoJSON is well loaded. And in this case the GeoJSON is removed when the function removeCableGeom() is called.


texts: ArcGIS 10: create perpendicular transects to stream at specified intervals

texts:
I am trying to do what my title says. The goal is to input a line, segment it in an evenly distributed way by a specified distance, convert the vertices to points, then draw lines perpendicular to the segment at a specified distance. These lines will then be converted to point in a similar fashion and return the maximum elevation value both east and west of each point. This whole process is meant to extract samples of an approximation of bank height for a river. I have found tools and processes that do the line segmentation (even though annoyingly it seems one cannot access the nifty Editor toolbar processes in python), but I am having trouble finding a script that makes the perpendicular transect lines.

Has anyone done anything like this before?

Note: I have found the script here, but it has too many malfunctions and I gave up trying to fix it (http://forums.arcgis.com/threads/49206-Perpendicular-transects-at-regular-intervals). Thank you!


texts: How to select values with decimals > 0 in a field?

texts:
I'm trying to write a script that would select the features which have the decimal value in the field 'Length' greater than 0.

Here's what I could think about it:

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface
import re

layer = iface.activeLayer()

selected_fid = []
r = re.compile(r"^\s*(?=.*[1-9])\d*(?:\.\d{1,2})?\s*$")

for feat in layer.getFeatures():
    if r.match(feat['Length']):
        selected_fid.append(feat.id())
    break

layer.select(selected_fid)

Traceback (most recent call last):
  File "", line 12, in
TypeError: expected string or buffer

I've read about using re for this purpose, here: https://stackoverflow.com/questions/8609714/regex-greater-than-zero-with-2-decimal-places but I'm not sure how to integrate this with the PyQGIS selection function.


texts: Is it possible to use GeoMesa Accumulo as a Geoserver Data store for real time data? i.e. Can Geoserver be configured to frequently read the data from Accumulo and update the layer?


texts: Create points inside polygon based on attribute in QGIS

texts:
I have a polygon layer with population and I need to create a new point layer based on that polygon layer.

The polygon layer has an attribute with the specific number of people living in each polygon.

How do I create X points (where X is the absolute number of people based on the attribute information) inside the boundaries of each polygon?


texts: Is there a a way to add a topology dataset to an mxd with arcpy?

texts:
I have a script that identifies certain layers from a folder containing several geodatabases.  Part of the scripts purposes, or intended purpose is to have it add a topology dataset to the dataframe. I have tried both of these methods below with no success.  It's as if the topology cannot be treated as a layer and therefore cannot be added.  Is there a workaround to allow for topology data to be added to the map?  If anyone has a solution for this I would be very grateful.

 elif name == "Topo_" + str(fileNameFull) + "_Topology":
        tempTopoPath = os.path.join(tempData,str(fileNameFull) + "_topoTemp")
        topologyPath = os.path.join(r'N:\29000s\29910\GIS\GI02\QC_Processing',fileNameFull,"FW_" + str(fileNameFull) + "_" + str(gdbDate) + ".gdb","Layers",str(name))
        topoLayer = arcpy.MakeFeatureLayer_management(topologyPath, tempTopoPath)
        topologyLayer = arcpy.mapping.Layer(tempTopoPath)
        arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, topologyLayer, "TOP")
        arcpy.RefreshActiveView()

Here is another method I tried.

elif name == "Topo_" + str(fileNameFull) + "_Topology":
        topologyPath = os.path.join(r'N:\29000s\29910\GIS\GI02\QC_Processing',fileNameFull,"FW_" + str(fileNameFull) + "_" + str(gdbDate) + ".gdb","Layers",str(name))
        topologyLayer = arcpy.mapping.Layer(topologyPath)
        arcpy.mapping.AddLayerToGroup(df, targetGroupLayer, topologyLayer, "TOP")
        arcpy.RefreshActiveView()


texts: CartoDB, Torque map sending events to JS when a market is rendered

texts:
I have a Torque map like this one and at the same time the bubbles are showing up I would like to render a growing graph with the number of bubbles accumulated. I don't expect CartoDB to do this for me but if CartoDB is able to send events to my JS functions I can do it.

So, the question is, is it CartoDB able to send events to JS functions when a new market is rendered?


texts:
When transforming the map from raster to vector, I get the following warning:
WARNING: Concurrent mapset locking is not supported on Windows and the result: The following layers were not correctly generated.
• C:/Users/bened/AppData/Local/Temp/processing_FnViAx/ca4332f52a844a92a3e192bcbd0d33c3/outpu


texts: OpenLayers 2 and D3.js:Set the right size of the SVG overlay

texts:
I am following an example of D3.js used in conjunction with OpenLayers 2 to display GeoJSON on a map.

http://bl.ocks.org/mbertrand/raw/5218300/quakes.html

In this example, the author set the bounds of the vector collection to be the size of SVG overlay:

var bounds = d3.geo.bounds(collection),
      path = d3.geo.path().projection(project);

var bottomLeft = project(bounds[0]),
      topRight = project(bounds[1]);

svg.attr("width", topRight[0] - bottomLeft[0])
      .attr("height", Math.abs(topRight[1] - bottomLeft[1]))

but in my codepen, I got a negative value in topRight[0] - bottomLeft[0].

http://codepen.io/shevawen/pen/ojqpRM?editors=001

An error shows in Chrome like this:

Error: Invalid negative value for <svg> attribute width="-134.70037333332988"

after set the svg.attr to a fixed value, its OK.

svg.attr("width", 2000).attr("height", 2000)


texts:
Is there a way to publish a JSON URL that displays a live feed of current bus vehicle locations in ArcGIS Online?

Here's the URL: https://transitdata.phoenix.gov/api/vehiclepositions?format=json


texts: Don't show empty categories in classified style

texts:
I have a geojson file containing line and polygon features.

I want to render them according to the field Feature Type which contains different types. When I use categorized symbology, it creates also categories which contain no elements. I think this is because some categories are only contained in the line feature layer, and some are only in the polygon layer.

Is there a way I can use categorized symbology, but leave out these empty categories?


texts: Set seed value for random points in vector grid

texts: I am using QGIS 2.18.11 and I have overlaid a 50m x 50m grid on a map (through Vector Grid) and have 98 grid cells in total. What I need to do is to generate 12 points randomly per grid cell that is “reproducible”. Setting the seed of “RandomPointspolygonfixed.py” didn’t work because the result shows the random points have similar pattern across the cells as shown in a snapshot below. Is there any way to do it?


texts: Using a projection with units in meters, but plotting with degrees on axes?

texts:
I've reprojected a map from WGS84 into into Transverse Mercator in an attempt to better represent the proportions of the region, but I would still like to have the axes in degrees, instead of meters.

This is the definition I'm using:

new.proj<-"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

and the map is plotted via the plot() function, producing the following figure.

I'm not sure if I'm grossly misunderstanding the point of using different CRS - I realize that Transverse Mercator uses meters as the units, but could it not be plotted as degrees if necessary?


texts: I'm starting to use QGIS desktop and many other tools like Python, Geoserver or PostGIS as a Geo-information's master student and I want to know if can I get a better performance for these in Linux, refering to speed boost (based in benchmarks), difficulty to work with and difficulty to install packages.


texts: Where are Esri LasDatasetBuilder examples?

texts:
Does anyone know of any working examples leveraging the new ESRI ArcGIS 10.1 Lidar classes i.e LasDataSet or LasBuilder class?

I've looked around on http://resources.arcgis.com/en/help/arcobjects-net/ but have not found anything other than the brief descriptions of individual class names and parameters they take in.  However, nothing in respect to working examples and how to tie everything together.


texts: Changing color of pixel using Google Earth Engine?

texts:
Is there a way to change the color of a single pixel in Earth Engine?

Generally, can run I run an algorithm for each pixel in earth engine in a region that analyzes each pixel's time-series and outputs a number that can be used to color each pixel in the region as a result?


texts:
I have a vector layer of 108374 lines  I want to split the vector layer into field names of the layer, which is 17. Now, I want to create a new vector layer containing one of the field names e.g. "OBJECTTYPE". From this field name I want to calculate the distance from a certain category to another vector layer.

I have tried to use "split vector layer" in QGIS, but it gives an error: "TypeError: the sip.simplewrapper type cannot be instantiated or sub-classed".

I'm using QGIS 1.8 Lisboa for Mac.

Thanks in advance.


texts: How to obtain timestamp and user id from .osm file in PostGIS enabled Postgres database

texts:
At present, osm2pgsql seems to
ignore the timestamp and user id of the elements, a serious drawback
for some applications. According to a
bug report this is not
solved by the addition of the --extra-attributes tag.
I've just tested this in osm2pgsql 0.83.0 and am wondering if there is a good way around it?

My plan B is to do pre-filtering with osmosis and then load the data into R using osmar. Still, it would be nice to have the PostGIS option working also.


texts: Which Custom Geographic Transformation Method to use for Nad83 to Nad 83 (CSRS)

texts:
I couldn't find a NTv2 file for transforming GCS NAD83 to GCS NAD83(CSRS) so I'm trying to create a custom transformation using the 'Create Custom Geographic Transformation' tool in ArcMap. The tool seems pretty straightforward but I don't know what transformation method I should use.

I'm using ArcMap 10.3 and am working in UTM zone 11N (in BC) if that is useful at all.


texts: Showing an error in the code (.eeObject.name is not a function)


texts: Loading a png image as a map with ggplot2/ggmap

texts:
I am trying to find a workaround for ggmap's missing support of world maps (i.e. you can't create any maps that show latitudes > 80°, due to idiosyncrasies in the mapproj package).
To a limited extend however, it seems possible to create empty world maps and save them as an image (png etc.), even if you can't use the ggmap object directly as one normally would in ggmap(get_map(...)).

That's why I'd like to load a png (ideally, one I created with ggmap) into ggplot2 and use that as a map instead. How exactly can I do that?

I am aware that you can load background images in ggplot2 (see this stackoverflow question). But I'd also like to plot points on my map, so it's important that the latitude/longitude values are mapped correctly.


texts: PDF File doesn't show with other Shapefiles

texts:
I have data as PDF file about the spatial categories of the German sea.

I sucessfully managed to import files about Windparks and military zones into QGIS by simply dragging them into the window and they matched with my other maps. My goal is to redraw the zones to get some vector files.

However when i tried the exact same thing with the PDF files of the protected sites ( Schutz und Eignungsgebiete) and other files, it only showed me empty layers.

I tried it with different CBS but i didnt seem to find the correct one even tho it says WSG 84.

When i import the files into an empty project they will show, however other files will then show empty.

Here's the dataset:

http://www.bsh.de/de/Meeresnutzung/Wirtschaft/CONTIS-Informationssystem/ContisKarten/NordseeSchutzUndEignungsgebiete.pdf


texts: MapServer in WMS mode for geotiff folder

texts:
I spent a few days trying to configure MapServer to work as a WMS server. I have a folder full of geotiff files (454 GB total, 62 MB per file). Each raster has embedded pyramid and I generated .vrt for that set.
This works in Geoserver, but performance is very poor. So now I am trying to do the same in MapServer (hopefully it will be much faster).

Serving shapefiles in MapServer works, but when I try to publish tif files, in QGIS I get this message:

  Invalid Layer: WMS provider Cannot calculate extent Raster layer
  Provider is not valid (provider: wms, URI:
  contextualWMSLegend=0&crs=EPSG:3908&dpiMode=7&featureCount=10&format=image/png&layers=Urbane%20zone%202012&styles=&url=http://localhost:8181/cgi-bin/mapserv.exe?map%3DC:%5C%5Cwmstif.map

Here is my map file content:

MAP
   EXTENT 6304500 4709990 6635000 5027700 #Tile Cache
   TEMPLATEPATTERN "."
   SIZE 512 512
   MAXSIZE 8192
   STATUS ON
   UNITS METERS
   NAME "Test map name"
   IMAGETYPE AGG_PNG
   OUTPUTFORMAT
    NAME "AGG_PNG"
    DRIVER "AGG/PNG"
    MIMETYPE "image/png"
    IMAGEMODE RGBA
    TRANSPARENT ON
    EXTENSION "png"
   END
   CONFIG "PROJ_LIB" "C:\ms4w\proj\nad"
   PROJECTION
      "init=epsg:3908"
   END
   WEB
      METADATA
         ows_title "Zoran example"
         ows_enable_request "*"
         ows_srs "EPSG:3908" #Recommended
         ows_abstract "This is an example WMS server from MapServer"
      END
      LOG "C:/ms.log"
   END
   LAYER
      NAME "Urbane zone 2012"
      STATUS ON
      TYPE RASTER
      UNITS METERS
      DATA "F:\geoserver_data\data\urbane_zone\index.vrt"
      METADATA
         ows_title "Urbane zone 2012"
         ows_abstract "rezolucija 25 cm"
         wms_extent "6304500 4709990 6635000 5027700" #Can use the map level extent values
      END
      TEMPLATE "blank.html"
   END
END

In QGIS I tried "Ignore axis orientation" options, but it doesn't help.

I am running MapServer on Windows Server 2012 R2, SSD disk, CPU i7-4790K, RAM 32GB. Raster files are on WD 3TB Green HDD.


texts:
I'm looking for an open source project or library to parse French addresses, i've been into github and i have tried a lot of those address parser but with no luck the result is not satisfying .
what i need is something like : https://www.address-parser.net/

which is able to identify extra info from an address .


texts: Vertices only showing up in a small portion of data view when editing

texts:
When I'm trying to edit features, or even digitize new features, the vertices or new lines only show up in the bottom right corner of the data view.  I've tried restarting ArcMap and making sure only certain features are selectable, but the problem persists.  I basically have to hover over vertices and then place them blindly.  Once I save edits, the vertices show up wherever I moved them.  Same thing goes for new features.  I basically have to click blindly and once I double click it all shows up.  However, a small square in the bottom right portion of my screen works as its supposed to.  It is within the program though, because if I shrink the window down and move it around it does not affect where vertices are seen.  I saw a similar post that was answered through changing some properties, but I'm not sure I understand, and their issue was for the entire screen.

Edited in response to below question:
Thank you for the response.  The link to the similar problem is below.                                              Vertices not showing in Edit session (ArcMap 10.5)                                                                                I'm editing vertices, but this also happens when I'm creating new features using the editing toolbar.  When you are editing a feature and you click on the edit vertices button, usually all of the vertices pop up and you can move them around, add more, delete them, etc.  I can find a vertice by hovering over certain places on the line (the vertices are actually there) I just can't see them.  I can even move them and edit them or delete them if I guess where they are correctly, I just cannot see any of this until after I save edits.  Hopefully that answers your questions.


texts: Model Builder ArcGis Iterate file and rename files

texts:
I want to go through a folder and subfolders to copy some files in a gdb. At the moment I have tried to do it using the iterator file and the Copy tool. The problem is all the files have the same name, then I would need to rename and give them a only name before to copy them in the gdb.

I am not too sure how I could rename them automatically, any advice?


texts: Text codes keep coming in as numbers, stripping out leading zeros

texts:
I'm trying to merge census tract data with a census tract shapefile in CartoDB.

Despite saving the tract ID codes as text/string in both .csv and .xlsx (including reimporting the .csv as text into Excel and resaving) formats, CartoDB keeps stripping out the leading zeros and treating the IDs as numbers when I import the files.

When I convert it back to a string in CartoDB doesn't put the leading zeros back in.

This is the only column the two datasets have in common, unlike some other county data I've worked with where I can merge on the Name of the polygon.

Any other ideas?


texts:
I'm using the Natural Earth shapefile to generate maps of ethnographic data. For every location where I have evidence of (for example) whether men or women had the job of brewing beer, I want to show a coloured marker.

My problem is that Mapnik insists on not letting the dots overlap, or even touch. This obscures a lot of the data, since many of the locations are very close to each other.

I've tried creating the markers this way:

ps = mapnik.PointSymbolizer()
ps.file = svgfile
ps.allow_overlap = True
ps.ignore_placement = True

and like this:

ps = mapnik.MarkersSymbolizer()
ps.file = svgfile
ps.allow_overlap = True
ps.ignore_placement = True
ps.spacing = 1
ps.placement = 'point'
ps.multi_policy = 'each'

but nothing seems to have any effect. I've tried adding and removing these properties in all possible combinations, to no effect.

For each type of marker I build a layer, style, rule etc like this:

r = mapnik.Rule()
r.symbols.append(ps)

s = mapnik.Style()
s.rules.append(r)

m.append_style('symbol_' + symbol.get_id(), s)

ds = mapnik.MemoryDatasource()
symbol_layers[symbol.get_id()] = ds

l = mapnik.Layer('%s_layer' % symbol.get_id())
l.datasource = ds
l.styles.append('symbol_' + symbol.get_id())

m.layers.append(l)

The actual points get added into the data source for their layer:

f = mapnik.Feature.from_geojson(json.dumps({
  "type": "Feature",
  "geometry": {
    "type": "Point",
    "coordinates": [float(marker.get_longitude()), float(marker.get_latitude())]
  }
}), ctx)
symbol_layers[marker.get_symbol().get_id()].add_feature(f)

How can I work around this? This is a major issue for me. I can't really use a visualization that hides half my data.

Edit: Setting clear_label_cache on the layer made mapnik ignore collisions from one layer to another. It's still avoiding collisions within layers, however.

Edit2: I was able to work around the problem by adding a separate layer and data source for each marker and using clear_label_cache. Not exactly elegant, but at least I get the output I need.


texts:
I'm trying to figure out why I can't copy a featureclass from one file GDB to another. I can copy the featureclass to within the source GDB. There's no name conflict in the target GDB. I don't think I've set any GDB constraints that would stop me. I'm running 10.0 ArcGIS. Any thoughts?

Here's what I've done to explore a bit and the error that I get.

[Dbg]>>> arcpy.env.workspace
u'C:\\Workspace\\wk_NHDPlus/NHDPlus03_Final.gdb/PRMS'
[Dbg]>>> arcpy.env.overwriteOutput
True
[Dbg]>>> arcpy.Exists('nsegment_flowlinesFL')
True
[Dbg]>>> arcpy.GetCount_management('nsegment_flowlinesFL')
<Result '5530'>
[Dbg]>>> desc = arcpy.Describe('nsegment_flowlinesFL')
[Dbg]>>> desc.dataElement.catalogPath
u'C:\\Workspace\\wk_NHDPlus\\NHDPlus03_wk.gdb\\Output\\nsegment_flowlines'
[Dbg]>>> arcpy.CopyFeatures_management('nsegment_flowlinesFL', 'C:\\Workspace\\wk_NHDPlus\\NHDPlus03_wk.gdb\\Output\\junk')
<Result 'C:\\Workspace\\wk_NHDPlus\\NHDPlus03_wk.gdb\\Output\\junk'>
[Dbg]>>> arcpy.CopyFeatures_management('nsegment_flowlinesFL', 'C:\\Workspace\\wk_NHDPlus\\NHDPlus03_Final.gdb\\PRMS\\junk')
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
  File "C:\ArcGIS\Desktop10.0\arcpy\arcpy\management.py", line 1943, in CopyFeatures
    raise e
ExecuteError: ERROR 000210: Cannot create output C:\Workspace\wk_NHDPlus\NHDPlus03_Final.gdb\PRMS\junk
Failed to execute (CopyFeatures).


texts: How to import or reproject SRTM 30 sec DEM

texts: I downloaded the global SRTM DEM from here. When I add the DEM as a raster layer (in my case topo30.grs file) the left/western edge of the file is at Greenwich and the right/eastern edge is at +360 degrees. This makes the layer run off the right side of the required area. The western hemisphere needs to be negative up to -180 degrees. I have attached an image which shows what I mean. Is there an easy way to import it with the right projection or reproject the SRTM. I don't understand why they created the DEM like this in the first place.


texts: Field format lost when exporting events layer to shapefile - long field type converted to double

texts:
I have a MS Access 2010 database with about 80 000 records in it. It is point data for animal locations. I would like to bring this into GIS and project the data into lambert conformal conic. I am experiencing some difficulties.

First I tried exporting the database in MS Access 2010 as a .csv before importing to ArcGIS 10.1, but there is apparently a known bug where Access will cut off your double fields to two decimal places. This doesn't work for me: I need those lat longs to be precise!

Then I tried exporting the data as an Excel .xlsx file type. Here all of my long integer fields from Access were converted to doubles (no good - for future processing I need those fields to be longs or integers in GIS).

Then I tried to directly import the x,y data using an OLE DB Connection. So far so good... the events file from the x,y data is imported as NAD 83 and the integer fields from Access are shown as being the long data type in Arc GIS.

Crossing my fingers I tried to then export this events file (NAD 83) as a shapefile (using the background data frame projection, lambert conformal conic, as to properly project the data for visualization). All of my fields that were previously classed as "long" in GIS have become doubles - and now they don't work in any of my post processing programs (namely, GME the geospatial modeling environment).

This is my first question - not sure if I can upload some sample data but I'm happy to provide it if I can figure out how to do it!

I have considered just bulldozing through them and making all new fields in GIS using the field calculator to change the field types, but I will be performing these kinds of operations regularly and since my data has about 50 fields I'd rather get this solved!

I am trying to export these data to a shapefile so that UDID is in long format.


texts: Make biodiversity heatmap in QGIS out of polygons and not points

texts: I would like to make a biodiversity heatmap in QGIS out of polygons and not points. I have a shape file made up of many overlapping polygons. I know how to make a heatmap with a set of points by using centroids, or by using a QGIS plugin called the Tombio FSC Biological Records Tool, but how do I do this with polygons? An example shape file of multiple sub-polygons can be downloaded from the IUCN (iucnredlist.org/resources/spatial-data-download), an animal conservation group, by creating a free account with them and downloading any of their shape files (e.g. marine mammals, mammalia, amphibia, reptiles). I want my map to look like this one, for example.


texts: Is licensed file geodatabase really secure?

texts:
I am planning to deliver our geo data in Licensed File Geodatabase format. This option was quite impressive as customer can have good look of the data with limitation and expiry date before purchasing.

Before using this Licensed File Geodatabase, I have performed some tests how other software handle this Licensed File Geodatabase. Here are my observations:

With in ESRI environment (tested Arcgis 10 and Arcgis 10.6) everything is perfect. You cannot open the data without installing license file.
In Safe FME also you cannot open data.
Most interestingly, Mr. QGIS 3.0  can open this secure data with possibility of export this data to any format :)

Here is the workflow to open this data in QGIS:

You should use QGIS OpenFileGDB (default driver) to open this data. If you have OGR FileGDB Driver installed on the machine then first uninstall it otherwise, you cannot open the data (Tested with QGIS version 3.0.0-Girona).

Here you can find the Licensed File Geodatabase for your test.


texts:
I'm trying to create a 30 meter, Int16 raster covering the entire continental US using Python/GDAL. This raster will ultimately be used to do raster algebra in ArcMap by downstream users.

Ordinarily I would use GTiff for this, but I'm getting the error message

ERROR 6: A 141014 pixels x 145929 lines x 1 bands Int16 image would be larger than 4GB but this is the largest size a TIFF can be, and BigTIFF is unavailable.  Creation failed.

So my question is twofold: (1) is BigTIFF the best format for creating a Arc-compatible raster, and (2) if so, how do I make it available?


texts:
I'm trying Geomesa Kafka and the GeoMesa Kafka Quick Start tutorial.
The producer writes features into Kafka data store and I can see new records in my Kafka topic.
However, my consumer
...
SimpleFeatureSource consumerFS = consumer.getFeatureSource(sft.getTypeName());
...
try (SimpleFeatureIterator iterator = consumerFS.getFeatures().features()) {
                System.out.println("Current consumer state:");
                while (iterator.hasNext()) {
                    System.out.println(DataUtilities.encodeFeature(iterator.next()));
                }
            }

every time returns empty iterator ( iterator.hasNext()=null ).
What could be the reason?
P.S. Writing records to the producer's datastore, causes following WARN message:
o.l.g.f.k.s.KryoUserDataSerialization$   : Skipping serialization of entry: USE_PROVIDED_FID -> true


texts: How to check point intersection with bbox made from two points?

texts:
I have two points which represents lower left and upper right point of BBOX.
I have one point which SHOULD be somewhere inside and I want to test it.

I have all of them as 'geometry(POINT, 4326)' data type.

When I try functions:

ST_MakeBox2D( geom4326ur, geom4326ll )
ST_Envelope ( ST_MakeBox2D( geom4326ur, geom4326ll )
Box2D

I get only linestring geometry connecting these two points.
I don't get 4 points of box which I need to transform them to polygon and then check with ST_Intersect if all points are inside their bounding boxes.

Can you help me solve this? I would like to have most effective solution - I am running this on millions of points.


texts: What is default ArcGIS raster calculator angular unit?

texts: Is the default arcgis raster calculator angular unit degrees or radians?


texts:
I am fairly new to the FME and don't have much hands on with all the available transformers.

I have a task to save two different sets of clipped data (via clipper transformer) into the multiple GDBs. I will try to explain in brief below:

I have a dump of road network data of two different timelines for one of the state in India. Also I have a polygon shapefile containing administrative areas of the same respective state.

I need the output containing clipped data of both the road shape into one GDB based on all administrative areas present in the polygon shape of the admin. In FME, I don't have any idea about how to create a writer that can write the output in multiple GDBs.

In short, I seek help from anyone who can help to set up my workbench which will help me to write the clipped data coming from two different timeline into the multiple GDBs (respective administrative areas).


texts: QGIS 3.14 with GRASS 7 error on batch dxf export

texts:
I tried to do a batch export of 19 geopackage layers to dxf and got an error I don't understand (I don't speak Python very well). First error of 19 pasted in below from my log, all of them the same.

Processing algorithm 1/19…
Algorithm v.out.dxf starting…
Input parameters:
{'GRASS_MIN_AREA_PARAMETER': 0.0001,
'GRASS_REGION_PARAMETER': None,
'GRASS_SNAP_TOLERANCE_PARAMETER': -1.0,
'input': '1_634370ab_e31a_4d68_91b9_bdbacd5dc470',
'output': 'H:/Data/Vector/map/layers/2.dxf'}

g.proj -c proj4="+proj=longlat +datum=NAD83 +no_defs"
v.in.ogr min_area=0.0001 snap=-1.0 input="H:\Data\Vector\map\layers\1.gpkg" layer="1" output="vector_5fb0b2ffecf5f2" --overwrite -o
g.region n=42.1083199772876 s=42.0816710148285 e=-74.4838672516787 w=-74.5197998630435
v.out.dxf input=vector_5fb0b2ffecf5f2 output="H:/Data/Vector/map/layers/2.dxf" --overwrite

Starting GRASS GIS...
Traceback (most recent call last):
File "C:/PROGRA~1/QGIS3~1.14/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 434, in processAlgorithm
Grass7Utils.executeGrass(self.commands, feedback, self.outputCommands)
File "C:/PROGRA~1/QGIS3~1.14/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 401, in executeGrass
feedback.reportError(line.strip())
TypeError: reportError() missing 1 required positional argument: 'fatalError'

Algorithm v.out.dxf failed…
Execution failed after 0.73 seconds


texts: How to create a raster delineating all drainage basins using pysheds python

texts:
I am going through the tutorial of "Pysheds" available in the following website
https://www.hatarilabs.com/ih-en/elevation-model-conditioning-and-stream-network-delimitation-with-python-and-pysheds-tutorial
Corresponding to all the steps mentioned in this tutorial, similar steps can be done in ArcMap (under Spatial Analyst --> Hydrology Tool).
Now the issue is I want to know is there any function available in Python Pysheds library to delineate all the drainage basins available in the DEM (after generation of flow direction map from DEM) as it is there in ArcMap's "Basin" tool under Hydrology.


texts:
I need to do some coordinate transformations on the browser and I cannot figure why I am getting wrong results with Openlayers 3.

Please see the below jsfiddle where I did a sample that defines the EPSG:2321 projection (using the proj4 definition from spatialreference.org) and then transforms a coordinate in that projection to EPSG:4326. Here's the fiddle:

https://jsfiddle.net/ak882y4b/

The result should be 32.93843, 39.96738 (5 digit display precision) whereas ol.proj.transform returns 32.938839053463106, 39.96748782887031. That is quite a large error margin.

What is wrong with this picture?


texts:
I want to change SLD of WMS layer. I think everything is well but SLD does not apply to layer.
Here is the code:

var SLD = '<?xml version="1.0" encoding="UTF-8"?>'+
'<StyledLayerDescriptor version="1.0.0" xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"'+
    'xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"'+
    'xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd">'+
    '<NamedLayer>'+
    '<Name>'+ layer.Name + '</Name>'+
    '<UserStyle>'+
        '<Name>' + layer.Name + '</Name>' +
        '<Title>A small red flag</Title>'+
        '<Abstract>A sample of how to use an SVG based symbolizer</Abstract>'+
        '<FeatureTypeStyle>'+
        '<Rule>'+
            '<Title>Red flag</Title>'+
            '<PointSymbolizer>'+
            '<Graphic>'+
                '<ExternalGraphic>'+
                '<OnlineResource xlink:type="simple" xlink:href="burg02.svg" />'+
                '<Format>image/svg+xml</Format>'+
                '</ExternalGraphic>'+
                '<Size>'+
                '<ogc:Literal>20</ogc:Literal>'+
                '</Size>'+
            '</Graphic>'+
            '</PointSymbolizer>'+
        '</Rule>'+
        '</FeatureTypeStyle>'+
    '</UserStyle>'+
    '</NamedLayer>'+
    '</StyledLayerDescriptor>';

    SLD = encodeURI(SLD);

    var source = layer.WMSLayer.getSource();
    source.updateParams({ STYLES: undefined, SLD_BODY: SLD });
    this.map.render();


texts: Select layer by attribute using expression qtcreator/pyqgis

texts:
I'm trying to make a function like ArcGIS's select using Python 3.6 , QGIS 2.8.6 and QtDesigner, so this is my code:

def execute(self):
        layers = self.iface.legendInterface().layers()
        selectedLayerIndex = self.dlg.layers.currentText()
        layer=None
        for lyr in QgsMapLayerRegistry.instance().mapLayers().values():
            if lyr.name() == selectedLayerIndex:
                layer = lyr
                qgis.utils.iface.setActiveLayer(layer)
                expr = QgsExpression(" \"BORONAME\" = 'The Bronx'")
                it = layer.getFeatures( QgsFeatureRequest( expr ) )
                ids = [i.id() for i in it]
                layer.setSelectedFeatures( ids )
                self.dlg.expression.setText(" ")
                self.dlg.lineEdit.setText(" ")
            break

I'm choosing a layer that has an attribute called "BORONAME" with the value 'The Bronx' when I try the code on Python console plugin it works fine, the results shown selected, but I got nothing when I use the plugin that I created.


texts: Adding blank rasters to mosaic in ArcGIS Desktop?

texts:
I've created a raster mosaic from OS 5 contour data to then process a number of viewsheds.

The issue is the area im working with is an island surrounded by sea so my mosaic is has no data in the areas of sea leaving holes when I make the viewsheds.

Is there a way to add blank/dummy rasters for these areas?


texts: Refreshing the GeoCache and Boundaries in OpenGeo suite

texts:
Currently if I add a new site i can query the layer and see where the point should be however it does not always display in Geoexplorer or OpenLayers preview. However sometimes it does show in Openlayers but not Geoexplorer. I think this is down to the GeoCache, and the Bounding Box around the data. I clear my internet history each time too.

Edit:

I have tried deleting the Geoserver cache, reseeding the layer at zoom levels 0-15. The odd thing is at zoom level 2, all points appear and then when I zoom in a progressively loose points on the map. However the data is still visible, if I query the layer?

Any suggestions

Thanks


texts:
I've a problem with GDALs TMS mini driver. I've created an XML file, which I can open in QGIS. I am using local tiles on my computer which I pre-rendered with Mapnik, the map covers entire country of Austria and the XML looks like this:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>file:///f:/austria/map/${z}/${x}/${y}.png</ServerUrl>
        <SRS>EPSG:3857</SRS>
        <ImageFormat>image/png</ImageFormat>
    </Service>
    <DataWindow>
        <UpperLeftX>-20037508.34</UpperLeftX>
        <UpperLeftY>20037508.34</UpperLeftY>
        <LowerRightX>20037508.34</LowerRightX>
        <LowerRightY>-20037508.34</LowerRightY>
        <TileLevel>17</TileLevel>
        <TileCountX>1</TileCountX>
        <TileCountY>1</TileCountY>
        <YOrigin>top</YOrigin>
    </DataWindow>
    <Projection>EPSG:3857</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
    <ZeroBlockHttpCodes>204,303,400,404,500,501</ZeroBlockHttpCodes>
    <ZeroBlockOnServerException>true</ZeroBlockOnServerException>
</GDAL_WMS>

The problem is: Whenever I zoom/pan into such window extents that is not completely covered with my tiles, the result is a scrambled image like this:

When I zoom/pan to an area covered with my tiles, the result is OK:

What I think is the issue: GDAL cannot merge and display the map correctly if there are missing tiles. If I was loading tiles from an http source, the server would probably return http request codes like 404 or 204. When configuring the TMS service in the XML, you have an option to define what to do when a tile is not returned/can't be found using these 2 tags:

<ZeroBlockHttpCodes>204,303,400,404,500,501</ZeroBlockHttpCodes>
<ZeroBlockOnServerException>true</ZeroBlockOnServerException>

The problem is, I don't know what is the response when accessing locally stored tiles and some are not found, or how do you specify what to do in such a case. I was also thinking of limiting the extent of the layer, but there are no options in the XML syntax for doing so.


texts: Open Source online software for making maps and displaying data?

texts:
I have noticed a number of online maps with data, however I wonder if there are any similar to umap that allows integration of data as in codap and offer features like openwebgis, ideally with an installation guide.

I have seen cartodb however installation seems to be a challenging task.
I need to set it up to show my boss so we can perhaps use it in my organization while not infringing on any licenses.


texts: why are all polygons of geojson-layer in OpenLayers 3 black?

texts:
I want to color each polygon of my geojson according to its Hex-Value at "DN". I tried this:

var map = new ol.Map({
    target: 'map',
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM()
        }),
        new ol.layer.Vector({
            title: 'added Layer',
            source: new ol.source.GeoJSON({
                projection : 'EPSG:3857',
                url: 'json2.json'
            }),
            style: new ol.style.Style({
                fill: new ol.style.Fill({color : '${DN}'})
            })
        })
    ],
    view: new ol.View(
        center:[52.5243700 , 13.4105300],
        zoom:2
    })
});

My geojsonfile (json2.json) has this pattern:

{
    "type": "FeatureCollection",
    "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

    "features": [
        { "type": "Feature", "properties": { "DN": "#FFFFA5" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 13.282823309930697, 52.660705106390886 ], [ 13.285901346780852, 52.660201769236203 ], [ 13.282823309930697, 52.660705106390886 ] ] ] } }
    ]
}

The polygons are all displayed but they are all black although there is not the HEX-code of black in any "DN"..


texts: OpenLayers : intersect between polygon and circle not accurate

texts:
When I look for intersect between circle and polygon. He finds intersect between these two, who are not visually intersected.
Everything should have the same projection ('EPSG:3857'). But I don't know how I can check the projection of each feature.
My code :
var circle = new Circle(transformedCoordinates, radius);
var circleFeature = new Feature(circle);
var circleSource = new VectorSource();
circleSource.addFeature(circleFeature);
var circleLayer = new VectorLayer({
  source: circleSource
});

var features = this.neighbourhoodLayer.getSource().getFeatures();
var nisCodesInCircle = [];
var circleExtent = circleFeature.getGeometry().getExtent();
for(var i = 0; i < features.length; i ++) {
  if(intersects(circleExtent, features[i].getGeometry().getExtent())) {
    nisCodesInCircle.push(features[i].get('id'));
  }
}


texts: Converting LLA (lat, long, alt) to ENU (east, north, up) in ArcMap?

texts:
I converted LLA (lat, long, alt) to ENU (east, north, up) in MATLAB using mapping toolbox but failed to view in the ArcMap. In ArcMap, it is viewing as single point.

Is there any other way to convert and view in ArcMap?


texts: The difference of transforming raster file in different cordinates between export and projection

texts:
I am using ArcGIS Desktop version 10.4.1.
I have a raster file with Geographic coordinate system WGS84. For my analysis, I needed to transform it into UTM projection. So I used the export function:
http://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/exporting-a-raster-in-arcmap.htm

I used The Data Frame (Current) option t export the raster dataset using the spatial referencing specifications of the data frame. So the coordinate of the current data frame (UTM) is different from  the coordinate of the raster file. It works very well.

However, when I use the project function of ArcMap to re-project from geographic coordinate to UTM and conducted the analysis, I noticed that the results are somehow different.

I googled all possible documents explaining the difference of export and project, I could not find a document explain the difference of those function

Can anybody  explain the difference of two functions?

It seems that export is doing the same thing as project when the coordinates of the source and the current data frame are different. But I want to make sure.


texts: Google Maps Javascript API, MarkerWithLabel z-index

texts:
First time poster, long time listener.

I am developing a Google Maps application and am using MarkerWithLabel to display the number of parking spots available in various parking lots around our town.

The z-index of the MarkerWithLabel appears to be messed up and no matter what I do to change z-indexes nothing seems to help.

As you can see from the screen shot.

From what I can ascertain (via developer tools in Chrome) is that the (P) icon has a z-index of 1, the label also has a z-index of 1 and the current position marker has a z-index of 10000. So I am unsure why the current position marker sits underneath the label but yet on top of the (P) icon.

Ideally I want the current position icon (in red) to sit over top of both the label and the (P) icons.

Can anyone help me with how z-indexes work with the Google Maps Javascript API as I am obviously missing something here.


texts: Reclassify raster in QGIS with GDAL raster calc

texts:
I'm trying to reclassify some categorical raster data in QGIS using GDAL raster calculator. My input Geotiff has four values, [26, 27, 31, 32] which I want reclassified to [1.25, 2.5, 3.75, 5] respectively.

I want to use GDAL raster calc so I can use GDAL's lossless compression creation options(the input raster is quite large).

The logic I've tried is the following, where inputraster.tif = A

(A=26)*1.25 + (A=27)*2.5 + (A=31)*3.75 + (A=32)*5

Where am I going wrong? I'm unfamiliar with Numpy but from the examples I've found in the docs and online this logic seems right, but the process crashes immediately.


texts: Setting the scale text for a spefice map inside a label

texts:
I am trying to archive the following:
I have a print setup with 3 map canvases.  I want to set  the scale text for a specific map canvas (lets say its ID is "main_view").

In QGIS 2.14. I can use the normal scale bar and set its style to "numeric" and link it to my desired map canvas. However I have dynamic text in which the scale text will be integrated like this:

Further Infos:
[%"open_info"%]
Created: [% day( now())%].[% month(  now() )%].[% year(  now() )%]
Scale: [%$scale%]

So I can not use the transformed version from the scale bar. The command [%$scale%] yields just a 0 because I dont know how to link to my map specific canvas.


texts: Joining multiple tables to multiple feature classes using ModelBuilder / arcpy

texts:
I have a list of tables (1,2,3) and would like to join the tables, tables by table, to a list of feature classes (1X,2X,3X). So that table 1 combines with feature class 1X, table 2 with feature class 2X and so on. The first table of the correspondent workspace should combine with the first feature class of the correspondent workspace, the second with the second and so forth.
I think it is similar to this problem where a solution is provided but in arcpy which I am not yet into! (https://community.esri.com/t5/python-questions/join-multiple-tables-to-multiple-feature-classes-arcpy/td-p/181549/page/2)
I experimented with lists, batch function and iterators but could not find a solution. The batch function seemed suitable, but it seems only to work with one list and the counterpart must be a single feature class/table as far as I could see.
I will show pictures of what from my point of view might be the closest solution but isn't working anyways. I thought I might iterate through the tables and the feature classes (one iteration in a submodel).
Or will this model try to combine every table with every feature class in my workspaces?

Edit: I tried a solution in arcpy. I runned it in ArcMap 10.7 python console. So I think it is using Python 2.7. (this is why I changed the f' strings from the solution in the link above [only available from Pyhton 3?] to "%" for inline variable substitution).The script runs but joins all tables to the last feature class in the list:
import arcpy

#Define Geodatabase Path
featureGDB = r"C:\C_M2\arcpy\list\list_Events_SBA_wp.gdb"
tableGDB = r"C:\C_M2\arcpy\list\list_distance.gdb"

#Create Feature List Variable
arcpy.env.workspace = featureGDB
featureList = arcpy.ListFeatureClasses()

#Create Table List Variable
arcpy.env.workspace = tableGDB
tableList = arcpy.ListTables()

#Set table Names to featureList Names equivalent
for table in tableList:
    tableNum = table [0:-5]
    #Find Table and Feature from Lists with same Name and join them
    for feature in featureList:
        if feature == "%tableNum%":
            inData = featureGDB + "%\feature%"
            inField = "OBJECTID"
            joinTable = tableGDB + "%\table%"
            joinField = "OBJECTID"
            arcpy.JoinField(inData, inField, joinTable, joinField)

The table names are "WP_20200723_001_dist" till "WP_20201010_005_dist". The feature class names are the same less the "_dist" anex: "WP_20200723_001" till "WP_20201010_005". Both lists have the same number of items.


texts: QGIS 3.0 not able to load processing plugin

texts:
I solved issues in the past opening my vector menu options with my processing plugin (in version 2.16), but recently upgraded my computer to a Mac OS 10.13 and my QGIS to version 3.0.1. My vector menu once again wasn't showing, but when I went to Plugins->Processing I got an error that I didn't understand which is shown below.

Can anyone help me through this?

  Couldn't load plugin 'processing'

  ModuleNotFoundError: No module named 'osgeo'  Traceback (most recent
  call last):   File
  "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 311, in loadPlugin
      import(packageName)   File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 664, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level)   File
  "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/init.py",
  line 28, in
      from processing.tools.dataobjects import *          # NOQA   File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 664, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level)   File
  "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/tools/dataobjects.py",
  line 51, in
      from processing.algs.gdal.GdalUtils import GdalUtils   File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 664, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level)   File
  "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/plugins/processing/algs/gdal/GdalUtils.py",
  line 35, in
      from osgeo import gdal   File "/Applications/QGIS3.app/Contents/MacOS/../Resources/python/qgis/utils.py",
  line 664, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level) ModuleNotFoundError: No module named 'osgeo'

  Python version: 3.6.5 (v3.6.5:f59c0932b4, Mar 28 2018, 05:52:31) [GCC
  4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)]  QGIS version: 3.0.1-Girona Girona, exported


texts: Reading raster datasets from file geodatabase with FME?

texts:
I have a file geodatabase (10.2) with geodatabase feature classes as well as geodatabase raster datasets.

I am using the Esri Geodb API reader in FME Desktop 2014 to read feature types from this geodatabase.

However none of the raster datasets stored in the geodatabase are being acknowledged by FME.

How can I read raster datasets from an Esri file geodatabase with FME?


texts:
I have tried to figure out why the layer is not showing up on the map. My code is as follows.

var markersLayerSource = {
    user_name: 'someaccount',
    type: 'cartodb',
        sublayers: [{
            sql: baseSql,
            https: 'force https',
            cartocss: cartoCss,
            interactivity: 'cartodb_id'
            }]
        }

        cartodb.createLayer(map, markersLayerSource, {https: true})
            .addTo(map)
            .on('done', function (layer) {
            console.log(layer);
            layer.getSubLayer(0).setInteraction(true);
            layer.getSubLayer(0).on('featureOver', function(e, pos, pixel, data) {
                // print data to console log
                console.log("Event #" + data.cartodb_id + ", name " + data.name + ", max population: " + data.pop_max);
            });
            //layer.getSubLayer(0).show();
            //cdb.vis.Vis.addInfowindow(map, layer.getSubLayer(0), ['cartodb_id']);
        });

Any one with pointers on why i specifically need to have cdb.viv.Vis.addInfowindow inorder for the layer to show up on the map? Am I doing something wrong if one needs to render without the addinfowindow?


texts: How to execute pgsql2shp from batch file

texts:
I am trying to execute pgsql2shp for a remote PostgreSQL db using batch file using the below code

for /f %%i in ('psql -U ihadmin -d ih_gis_production -h <xxxxx> -c "select f_table_name from geometry_columns where f_table_schema='mls_dataset';"') do (
set file_name=D:/Processing/DB_Dump/airdrie/mls_datatset/%%i
set table_name=%%i
pgsql2shp -h <xxxxx> -f %file_name% -u ihadmin -P xxxxx ih_gis_production mls_dataset.%table_name%
)

This will get table names from from the above select query and loop and export those tables to shapefile.
When I run this from batch file, I am seeing pgsql2shp help options in command prompt. Any correct way to call pgsql2shp via batch?


texts:
I am trying to create a line feature as the result of another line in PyQGIS. The new line should start, where the old "ends", face the "same direction" and should have a defined length of X.

Unfortunately, the direction of both lines is not the same using the following code. What is wrong?:

p1 = QgsPointXY(14.18064, 50.63963)
p2 = QgsPointXY(14.15027, 50.66042)

d = QgsDistanceArea()
d.setEllipsoid('WGS84')
degrees = p1.azimuth(p2)
rad = math.radians(degrees)
# distance is 1000m
d_upfront = d.computeSpheroidProject(p2,1000,rad)
line_layer = QgsVectorLayer("Linestring?crs=4326", "Measurements AUTO", "memory")
QgsProject.instance().addMapLayer(line_layer)
line = QgsGeometry.fromPolylineXY([p1, p2])
feat = QgsFeature()
feat.setGeometry(line)
line2 = QgsGeometry.fromPolylineXY([p2, d_upfront])
feat2 = QgsFeature()
feat2.setGeometry(line2)
provider = line_layer.dataProvider()
provider.addFeatures([feat,feat2])

provider.forceReload()
line_layer.reload()
line_layer.triggerRepaint()

The example shows the difference:


texts: How do I remove TIR bands from current stacked ASTER

texts:
How do I remove TIR from already stacked ASTR that has VIS and SWIR bands. I want to do it in ENVI 5.0. Is there a workaround to do this ?

The bands I stacked has VIS, SWIR, and TIR.


texts: Weighted-area averages using dissolve for SINGLE PART features

texts:
I'm trying to compute area-weighted averages of a field when dissolving a feature layer. I thought the solution would be to:

create a new field that is the product of the original field value and the original shape area: new_field = !original_field!*!Shape_Area!
Perform the dissolve on field my_class with the SUM statistic for new_field: SUM_new_field = Sum(new_field)
Divide new_field by the dissolved feature area: aw_field = !SUM_new_field!/!Shape_Area!.

I additionally want to use single-part features only, so I uncheck the "multipart features" option in Dissolve. However, the result is that the values of aw_field is identical for every feature of a given value of my_class. It appears that ArcGIS first dissolves into multi-part features, does the statistics computation on the multi-part features, and then splits the multi-part features into single-part features and simply assigns the multi-part feature value of the statistics field to each resulting single-part feature.

What is the most efficient work-around to get an area-weighted average for dissolved single-part features? I can see a procedure where you first dissolve into single-part features with no statistics, use the resulting layer as an identity feature or for zonal statistics on the original layer, and then join those results back to the dissolved feature... but that seems obnoxious and needlessly computer-intensive. Does anyone have a better idea?


texts: How to add map into browser using OpenLayers 3 and GeoServer?

texts:
I am new to OpenLayers 3, no error is showing in console, but the map is not showing in browser.

My code is:

    <title>mymap</title>

     <link rel='stylesheet' href='css/ol.css'>

Please help


texts: Subtracting two raster layers based on condition in QGIS

texts:
I have 3 raster layers in QGIS. I want to subtract two from each other based on a condition in a third layer.
'Layer1': has no data in some parts and I used Close gaps to generate 'Layer1_fill'. What I want to do is to subtract the filled values in 'Layer1_fill' from 'Layer2' whenever 'Layer1' has no data.
If 'Layer1' has no data and 'Layer2' has data, then 'Layer2' - 'Layer1_fill'
To do so, I converted nodata in 'Layer1' to 0 and called it 'layer1_null0'.
I tried this expression in Raster Calculator but it didn't work
("layer1_null0@1" = 0 AND Layer2 != 0) * ("Layer2@1" - "Layer1_fill@1")

Any suggestions will be appreciated.


texts: Georeferencing of PNG and Convert to TIFF using Python

texts:
I have a PNG image which I am trying to georeference using four corners of coordinate information of the defined input image.
When I am trying to convert PNG to TIFF, I am getting 3 bands (R,G,B) raster. Then I have copied the metadata from the original raster and tried to save it as GeoTIFF along with copied metadata using rasterio library.
In the final output, I am getting only one e band as R band. but I want all the bands in my GeoTIFF output.
Below, I have attached, input PNG and output GeoTIFF images. While converting PNG to GeoTIFF, the blue and black features are missing. only yellow is capturing.
Below is my code,
import rasterio as rio
import numpy as np

#Input png image, to convert as geotiff
img = rio.open('C:/Users/Documents/Study_Area_2.png')
img = img.read([1,2,3])
img = img.astype('uint16')

# Input image for coordinate reference
with rio.open('C:/Users/Study_Area_2_RGB.tif') as src:
    naip_data_ras = src.read()
    naip_meta = src.profile

#output images
with rio.open('C:/Users/Georeferenced.tif', 'w', **naip_meta) as dst:
    dst.write(img,[1,2,3])

Input Png image

Output GeoTIFF


texts: Is there a way to classify simple non-spectral images in ArcGIS?

texts: I was wondering whether there's a way to classify simple non-spectral TIFF images in ArcGIS (maybe based on similar colour values of image pixels like graphic software applications such as Photoshop can do)?


texts: I am exporting my QGIS project with categorized information but for some unknown reasons, the layer only shows when you zoom in to a larger extent and not when you zoom out or zoom to the whole extent. Other non categorized layers works perfectly fine both symbols and text. I am using OpenLayers which I would like to stick to. The categorized symbols works fine with leaflet but I do not want to use reaflet for other reasons-reaflet also give text issues.


texts: Cell Size and nonsensical results from GRASS r.walk

texts:
I've been wrestling with two main problems the past six hours. I am using 30m resolution aster/srtm merged dem. When I use the r.walk function to calculate a walking time surface from a point, the output looks good, but when I check the values, they make no sense. That is, the time it takes to travel is way underestimated, often by a factor of 3 or more. I am assuming that the results are in seconds, so I divide that by 3600 to get the value in hours. Actually, given the terrain (Andes), the times are impossible because the velocity implied by the results is way more than 5km/hr (about 20km an hour). I tried using just the dem and setting lambda to 0 for the friction. I tried using a reclassed slope as the friction surface. Both yield similar results.

The second issue is that when I zoom in to check the cell size, it's about 1.5km on each side! I don't know how that's determined because all the raster inputs are 30m resolution. The extent of the GRASS region isn't that big, either. Maybe 150 by 150km. I can't find any resources that address this issue.


texts:
I am using the new valhalla plugin on QGIS 3.10. I am trying to create isochrones, but I receive the following error:
KeyError 'metric'

I have checked the properties of the project (metric) and the routing option works perfectly.
I am using the Mapbox API.
What could I be doing wrong?
The following shows up in the log:
The error log: {"error_code":110,"error":"Insufficiently specified required parameter 'locations'","status_code":400,"status":"Bad Request"}
2020-12-18T06:57:16 WARNING Traceback (most recent call last): File "C:/Users/Andrew Rickard/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\valhalla\gui\ValhallaDialog.py", line 320, in run_gui_control for feat in isochrones.get_features(str(i), isochrones_ui.costing_options): File "C:/Users/Andrew Rickard/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\valhalla\common\isochrones_core.py", line 138, in get_features metric = isochrone['properties']['metric'] KeyError: 'metric' – Andrew Rickard


texts: Draw a skyplot with calculated ECEF values (in meters)

texts:
I have ECEF values that I calculated from ephemeris subframes 2/3.  Using this information, I'd like to place these values on a skyplot based on my current location.  I'm using a C# picture box for my skyplot image, and a custom control to show the dot and SV PRN.

I am struggling in these areas:

How to merge the ECEF information with my current location to provide a useful x,y for my skyplot.
How to determine the bounds of my drawing.  Currently i estimate that from edge to center of my drawing is 90 pixels.  So i attempted to do some scaling 90 / 20,000,000 (max distance for satellite, guessed), then multiplied that by my ECEF.x, and ECEF.y values.

I'd appreciate any help.
Thanks,
Jason


texts:
I currently have a coordinate based map of the US and it's states for topojson (https://github.com/mbostock/topojson).

The coordinates are a JSON array (as topojson works), data looks like this:

[-123.233256,42.006186], [-122.378853,42.011663]...

Now what I'd like to do is combine states into regions, outside of doing this by hand (and that idea makes me feel sick), are there any tools which I can use to accomplish this?

I have QGIS but have no clue where to even start (let alone how to combine states in QGIS).


texts: How to convert a polygon to a grid of polygons using QGIS?

texts: For example, the city of Rotterdam represented as a hexagonal grid.


texts: How to convert polygons to centriods without losing attribute names?

texts:
I have a GeoJSON file of polygons and I need to convert it to centroids. Finding the centroid is very easy in QGIS.

Problem

My attribute names are very long and have a colon between them, e.g. building:structure and shape:elevation. Converting to centroids results in a shapefile format. It gives the geometry but changes the attribute names, truncating them and converting colons to underscores, e.g. building_l and shape_elev.

I need the attribute names in the centroids to be the same as in the original polygon.


texts: Explaining scale decrease inside parallels and outside increase for Lambert conformal conic projection?

texts:
For a Lambert conformal conic projection, when projecting the cone on the globe, two parallels are formed and the scale decreases between the two parallels and increases outside them.

Why does this phenomenon of decreasing scale occur inside the parallels and increasing scale occurs outside them?


texts: Layers from SpatialLite Database are not displayed correctly in QGIS

texts:
While working a lot with layers from spatialLite databases I face a returning bug in which layers from databases aren’t displayed correctly in the project.
The error affects persisting layers as well as virtual layers, which are created by sql query.
Layers are not displayed. Although related tables are shown correctly in the db-manager, related attribute tables are displayed empty, while the number of Objects is shown correctly.

The error suddenly appears at any time, while working on a project. Thus all layers from the database are affected, the only opportunity for fixing is to start a new project.
While searching for solutions on the internet, zooming in and out may sometimes fix any displaying errors in virtual layers. While trying out still the problem remains.


texts: Running two versions of ArcGIS Desktop on same machine?

texts:
I have the full installation Media for 10 (9.3/SQL05 already loaded). And two editor licenses (one loaded and a new license file). What manner of wizardry is involved in getting around the conflict detection where I can run 9.3/SQL05 and 10.0/SQL08 on the same machine?

I have a gdb that I sync with at a far remote local that can't be helped at this time. But I need some of the new functionality in 10 for other projects.

I found this post:
http://blogs.esri.com/Dev/blogs/arcgis/archive/2009/10/23/Installing-ArcGIS-9.4-_2800_and-can-you-leave-9.3.1-on-the-machine_2900_.aspx

but there must be a way to get both working on the same machine.

I find no real help at all on the net. Most people just upgrade to 10.


texts: How to determine the version of an SDF file?

texts:
I'm looking for a tool or a method to reliably determine which SDF version a given file  is implementing.

So far I've been using Autodesk's SDF Loader with the /info flag, but it only supports SDFv2.

Unfortunately I have a large number of SDF files of various version where i need determine which ones need to be upgraded to SDFv3. As FDO Toolbox only supports SDFv3 and we are upgrading from MapGuide 6.5 which only supports SDFv2.

So in summary I'm looking for a reliable way to determine the major and minor version number of the SDF file hopefully in a single tool. (I'm aware of FME, leaving it as a last resort due to cost)


texts: Why is my network route performing oddly

texts:
I'm running Arcmap 10.2.2 and am messing with network analyst for the first time.  I used QGIS to download an OSM layer that I imported into an enterprise geodatabase.  I created a network with all the roads and set connectivity to Any vertex.  I selected two points which can be seen in the attached image.  I ran a route with just these two points and it goes way out of the way.  The distance should about 20 miles and it comes out as 118 miles.

Any thoughts how I can correct this issue?


texts:
Try to do this

and you'll get an error message. However you can easily get statistics on [FID] using Summary Statistics tool. Somebody has a good explanation of this? Also note that from within the table [FID] is treated like non-numeric 'field'.

Extended question:

Statistics on [FID] can be very quick way to complete many tasks. For example find duplicates stored in other field. In example below field [COORD] stores string equivalent of point coordinates:

If I apply Summary Statistics tool, case field [COORD], statistics field [FID] and join resulting table back to original points table using [COORD]. I can find duplicate points using query

"pnts.FID" <> "st.FIRST_FID"

There are very many other tasks that can be accomplished using sequential record number.

My question is WHY one cannot perform statistics from TABLE VIEW(window), which is very fast way to perform statistics on the field, and at the same time this can be done using Summary Statistics tool. That is not so quick, because it is a long way down the tree System Toolboxes-Analysis Tools-Statistics.


texts: Defining custom projection using baseline survey coordinates

texts:
I have CAD data that I need to incorporate into a GIS analysis. The CAD data does not have a defined projection, but the drawings include a table of survey baseline points. I need to know if I can build a custom projection based on this table:

Is this enough information to construct a custom projection string, or do I need more data?

The documentation references a "Borough of Queens Coordinate System", which turns up nothing in Google. My hope was that it's some obscure projection only used in New York City but that it could be recreated using the above table.


texts: Applying SLD to vector layer in OpenLayers 2?

texts:
I'm trying to apply an SLD style to a vector layer in OpenLayers. I already took a look on google and here and I got a lot of information, but still I don't manage to achieve the result I'm looking for. Here is the code:

var map;

function init() {

    OpenLayers.ProxyHost = 'proxy.php?url=';

    map = new OpenLayers.Map("map", {});
    var park = new OpenLayers.Layer.WMS('Park',
                       'http://localhost:8080/geoserver/PK/wms',
                       {
                           layers: 'PK:park'
                       },
                       {
                           isBaseLayer: true,
                       });

    map.addLayers([park]);

    var roads = new OpenLayers.Layer.Vector('Roads', {
        protocol: new OpenLayers.Protocol.HTTP({
            url: 'roads.xml',
            format: new OpenLayers.Format.GML({
                extractAttributes: true
            })
        }),
        strategies: [new OpenLayers.Strategy.Fixed()]
    });

    map.addLayers([roads]);

    //SLD
    var format = new OpenLayers.Format.SLD();
    OpenLayers.Request.GET({
        url: 'style/roadsStyle.xml',
        success: function (req) {
            sld = format.read(req.responseXML || req.responseText);
            styles = sld.namedLayers.interpreted.userStyles;
            roads.styleMap.styles.default = styles;
        }
    });

    LayerSwitcher = new OpenLayers.Control.LayerSwitcher();
    MousePosition = new OpenLayers.Control.MousePosition();
    PanZoomBar = new OpenLayers.Control.PanZoomBar();
    map.addControls([LayerSwitcher, MousePosition, PanZoomBar]);

    if (!map.getCenter()) {
        map.zoomToMaxExtent()
    };

};

Pretty simple one, but, with this, the console returns error and the style is not applied and I don't know what I did wrong.


texts: OpenLayers - remove all layers of map apart from base map layer?

texts:
What i need is to delete all layers from the map without removing the base one (in my case OSM)

Here is some code

    new Map({
      layers: [
        new TileLayer({ source: new OSM() }) // <== the one i dont want to be removed
      ],
      controls: [
        new Zoom({ zoomInLabel: '', zoomOutLabel: '' })
      ],
      view: new View({
        center: proj.fromLonLat([center.lon, center.lat]),
        zoom
      }),
      target: div
    });

I can iterate over them with a forEach and applying map.removeLayer(layer); but i can't distinguish it with the other layers.


texts: Is it possible to access the building textures from the Google Earth API?

texts:
I'm trying to get building textures from Google Earth and apply them to my own models. My overall goal is to be able to walk around in virtual reality.

So far, I've found some information from the samples provided by google. In one example, the textures are held in texture folder.

http://kml-samples.googlecode.com/svn/trunk/resources/textures/

I doubt it's this simple, I imagine you need to either use KML or call Google Earth's LAYER_BUILDINGS.


texts: Landsat Observation Zenith and Azimuth Angles

texts: I am doing a project with Landsat which needs to find out 2 Landsat images of the same area for different time, have both the same observation zenith and azimuth angles. But I can not find any information about satellite observation zenith or azimuth angle of Landsat images. Any one has an idea where I can find that information?


texts:
I was wondering if anybody would be interested in collaborating?

My Idea is using QGIS to analyze and quantify Immunofluorescence raster images of cells from human serum.


texts: Calculate attribute table field using a near function (or a similar function)

texts:
I'm having an attribute table with a field called 'name' which represents the streetnames of an area. The streetnames are derived from OSM, but unfortunately not all street (or parts of streets) have a streetname related to them. I believe that the reason for this is due to the difference of streettype (primary, secundary etc). What I want to do is to fill these missing streetnames with the streetname of the street that is closest to the street with the missing name.

I believe the best option to do this is to calculate the field using a near function, but I'm not that experienced in Python or other programming language. Therefore, I'm not really sure if my theory is plausible and/or the right one.

I have added a screenshot of the attributetable to give you an idea:


texts: Leaflet not registering touch events on static map?

texts:
I am attempting to create a tool that will present users with a web map, which cannot be zoomed, panned, etc. Just a static map, dynamically presented, ie at a semi-random location for the purpose of testing various aspects of map design/legibility.

I want users to be able to trace a line on the map with their finger using a touchscreen interface. I've been using leaflet, and things were going fine while I was developing on my laptop (click -> mousemove -> click events to draw a line), but when I went to test the application on my tablet, (touchstart, touchmove, touchend events??) nothing seems to happen.

Is leaflet somehow disabling these events? I am trying to add them to the map in the same way I had added click events, with the .on method of the map object. No errors show up, it's just that nothing happens.

The (hopefully) relevant part of my code is below.

$m = L.map('map',{
    'center': $mapCenter,
    'zoom': 17,
    'zoomControl': false,
    'attributionControl':false,
    'animate':true
});
// disable all map movement
$m.dragging.disable();
$m.touchZoom.disable();
$m.doubleClickZoom.disable();
$m.scrollWheelZoom.disable();
$m.keyboard.disable();
// use OSM tiles for testing
L.tileLayer($tilesource).addTo($m);
// create and then add A and B
L.marker($A,{icon:$Aicon}).addTo($m);
L.marker($B,{icon:$Bicon}).addTo($m);
// add event listeners for user interaction
$m.on('touchstart',trackCursor);
$m.on('touchend',closeTrack);
// monitor cursor position
$m.on('touchmove',function(e){$cursorLocation = e.latlng});

I am using leaflet 1.0.2 and the Firefox 50.1.0 on Android.


texts:
I am trying to follow this tutorial that has buttons to toggle (all points are on the same layer)

CartoDB Toggle map tutorial

I try to follow as close I can for my own purpose but I always get the error under line 71 "uncaught reference Error: layer is not defined" (see attached screen shot)

Can anyone tell me what I did wrong? I tried messing around with all the Layers, Sublayers, layer, sublayer...but non worked.

Here is the script

beginning of script tag

var layerUrl = 'my cartodb url';
    var sublayers = [];

  function main() {
    cartodb.createVis('map', layerUrl, {
        shareable: true,
        title: true,
        description: true,
        search: true,
        tiles_loader: true,
        center_lat: 22.273904,
        center_lon: 114.174510,
        zoom: 14
    })

    .done(function(vis, layers) {
      layers[1].setInteraction(true);
      layers[1].on('featureOver', function(e, pos, latlng, data) {
        cartodb.log.log(e, pos, latlng, data);

      });

      map = vis.getNativeMap();

        var subLayerOptions = {
  sql: "SELECT * FROM final_landreclamation_all" }
        var sublayer = layer.getSubLayer(1);
            sublayer.set(subLayerOptions);
            sublayers.push(sublayer);

    })
    .error(function(err) {
      console.log(err);
    });

    var LayerActions = {
all: function(){
  sublayers[1].setSQL("SELECT * FROM final_landreclamation_all");
  return true;
},
Saiwan: function(){
  sublayers[1].setSQL("SELECT * FROM final_landreclamation_all where district_en = 'Saiwan'");
  return true;
},
Sheungwan: function() {
  sublayers[1].set({
    sql: "SELECT * FROM final_landreclamation_all where district_en = 'Sheungwan'",
  });
  return true;
}

}
    $('.button').click(function() {
    $('.button').removeClass('selected');
    $(this).addClass('selected');

LayerActions[$(this).attr('id')]();

});

  }

  window.onload = main;

end of script tag