texts: Condition statement not returning correct value

texts: I may be missing something obvious here, but I cannot figure it out. In my attribute table under the area_ft_su column, I have a bunch of NULL entries. I want to change those NULL entries to zero without having to manually change each one. The output preview looks correct, but when I click "Ok", nothing happens. Nothing is updated. Am I missing something or is my syntax wrong?


texts: Adjust leaflet routing machine draw (animate) route speed

texts:
I use leaflet routing machine and it's working fine with me.

Here is my code:

<!DOCTYPE html>

    <meta name="viewport" content="width=device-width" />
    <title>Test</title>
    @Scripts.Render("~/bundles/jquery")
    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.4/dist/leaflet.css"
          integrity="sha512-puBpdR0798OZvTTbP4A8Ix/l+A4dHDD0DGqYW6RQ+9jxkRFclaxxQb/SJAWZfWAkuyeQUytO7+7N4QKrDh+drA=="
          crossorigin="" />

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.2.0/dist/leaflet.css" />
    <link rel="stylesheet" href="https://unpkg.com/leaflet-routing-machine@latest/dist/leaflet-routing-machine.css" />

It draws the line immediately after I load the map, but I want it to draw the line slowly like this sample for example


texts:
Does anyone know of a list of countries with their associated official language codes? I haven't found anything by Googling yet.

What I'm looking for is something like this:

Country name/code              Language code
[...]
Bulgaria (BG)                  bg
Belgium (BE)                   nl, fr
[...]

Some further clarifications:

How is there not a correlation between a country and its official language(s)?

In my application this is relevant because I want to redirect someone who accesses it from the US to ${path}/us/en/, so I need to get the official language tag. I get the country code from a geolocation web service and I need the associated (official) language code.

For example, I have everything on my computer set to english, but I live in Romania. When I access it I don't want to get redirected to the UK version, ro/en does not exist, so I want ro/ro.

I wanted to see if there was an already aggregated list so that I wouldn't have to do it manually. If there isn't, I'll build one :).


texts: Reclass a landuse shapefile based on another shapefile using R

texts:
I have a landuse shapefile that I want to modify based on another shapefile.

In this example I have three classes in landuse shapefile A.

Say I'd know that everything within one of the small shapefiles, called B, outlined in black belonged to class four.

Is there a way to create a new shapefile C based on the landuse shapefile A and a new shapefile B that assigns class 4 to everything that falls in shapefile B's geometry?

I tried workarounds for that, substracting B from A and then combining them back together, but I couldn't find a method that adds both polygon fields together AND the attribute table without creating things like area.1, Class.1, area.2, Class.2 and so on.

I'm really stuck with this but I believe there must be a not-so-hard solution to this.

Here are two example shapefiles, landuse and industry.
I want to produce a shapefile where all that is covered by industry is classified as industry and the rest of the landuse stays the same


texts: Create repeating pattern of raster or polygons in ArcGIS

texts:
I want to create a repeating pattern of raster OR polygon features based on an existing sample, so that the sample fills up a larger area in a repeating pattern. For example, with a sample of 10x10 units I want to fill an area of 100x100 units with 100 samples.

The goal is to create a mask of repeating pattern, similating a specific landscape type or pattern. This mask is then to be clipped by polygons to fit in a map, used for landscape change and connectivity analysis.

I have browsed around existing tools, but the only way I have found so far is to manually copy+paste the feature x times, risking both inconsistency and topological errors.


texts: Converting MDB File With Unicode Characters

texts:
I'm using C# wrapper of GDAL and I'm having a problem converting MDB file that contains unicode characters.

Lets take for example conversion to KML file,
non-English characters are replaced by '?'

How do I solve it?


texts: Clip raster using multiple country shapefiles

texts:
I currently use gdalwarp to clip a GeoTIFF file using a single shapefile. However, I would like to have the option of creating a clipped GeoTIFF file using multiple countries together.

Is using GDAL the proper way to approach this?


texts: Is there a way to access the 'data' from the cartodb.createLayer function?

texts:
Can I access the data portion of the SQL queries used in the cartodb.createLayer function without the use of special event handlers?

I am making a choropleth map, and categorizing the data with Jenks natural breaks. So, in order to determine the breaks needed to set the right values in CartoCSS, I need to access all values in the data, run the Jenks process, and set the proper CartoCSS. Currently I do this via an independent SQL query prior to running createLayer:

var sql = new cartodb.SQL({ user: 'user' });
var sqlQuery = "SELECT * FROM table;"
sql.execute(sqlQuery)

This works fine, but I would rather not run the independent SQL query, if I can somehow extract that information from the createLayer function.

I set up the layer source with four sublayers:

var layerSource = {
  user_name: 'user',
  type: 'cartodb',
  sublayers: [
    {
      sql: q1,
      interactivity: 'cartodb_id',
      cartocss: s1
    },
    {
      sql: q2,
      interactivity: 'cartodb_id',
      cartocss: s2
    },
    {

      sql: q3,
      interactivity: 'cartodb_id',
      cartocss: s3
    },
    {

      sql: q4,
      interactivity: 'cartodb_id',
      cartocss: s4
    }
  ]
}

Then I create the layers:

cartodb.createLayer(map, layerSource, {cartodb_logo: false})
.addTo(map)
.done(function(layer) {
  for (var i = 0; i < layer.getSubLayerCount(); i++) {
    subLayers[i] = layer.getSubLayer(i);
  }

etc...

I know of no way to extract the attribute values from the four SQL queries that are called to the make the layer - from neither the returned 'layer' nor any of the elements in subLayers.

I can access the 'data' from a featureOver or featureClick event handler on a sublayer, e.g.

 subLayer[0].on('featureOver', function(e,pos,latlng,data){
    $('.leaflet-container').css('cursor','pointer');
    if (data.cartodb_id != highlight.cartodb_id) {
        highlightArea(data);
 }

but I need the data before that - to set up CartoCSS and create a proper legend.

Any ideas?


texts:
PostgreSQL's ST_Distance_Sphere function calculates the distances on earth:

float ST_Distance_Sphere(geometry pointlonlatA, geometry pointlonlatB);

In its documentation, there is a peculiar note:

  This function currently does not look at the SRID of a point geometry and will always assume its in WGS 80 long lat.

My geometries are stored in WGS 84, which seems to be the de-facto standard for storing geospatial entities in PostgreSQL, an the reference coordinate system used by the GPS system.

Why does ST_Distance_Sphere uses WGS 80? Should I expect any significant error due to the coordinate system difference?


texts: Conditional drop-down menu list in a python script tool

texts:
I have created a script tool from some python code, and in a previous question, I learned how to enable a drop-down list in this tool by changing the Validation code - see link and image below:

Setting up Drop-Down list in Parameters of Python script tool?

Now I would like to have the second user-input in the tool have its own drop-down list whose members change depending upon the choice in the first input - see image:

For example, if the user chose the town "Winooski", I would like to store that input as a variable, maybe called 'name', which might then be included in a block of code like this:

name = self.params[0]
roads = "Z:\\OPS\\TechnicalServices\\Culverts\\GetCulverts\\GetCulverts.gdb\\%s" % name
rows = arcpy.SearchCursor(roads)
self.params[1].filter.list = sorted(list(set(r.getValue('CTCODE') for r in rows)))
del rows

In the 'roads' directory defined in this code, I have a feature class for every town name, which only includes the roads for that town - so the tool will have the user choose a town, and they will then be presented with a list of state roads in that town to choose from. I tried sticking the above code within the Validation code just beneath the block of code that solved my first problem (in the first image above). I'm assuming that the Validation code only runs when you first open the tool. Any clues as to where I should be sticking some code to do this (maybe in the code for the tool itself..?), or how I should be modifying the small block I wrote?


texts: QGIS credentials storage - security issue?

texts:
QGIS is the main GIS application used at my workplace. Access to the plugin repository has recently been removed due to a perceived security risk relating to the way QGIS stores proxy settings, including user passwords, as below (from the QGIS configuration page).

Is there a way to stop QGIS from saving and retaining records of these details, and instead requiring users to enter them each time they access plugin repositories?


texts: QGIS: attribute table in print-composer, format numbers and sort according to same field

texts:
In the print layout of QGIS (I use 3.12.2 on Win) you can add an attribute table. I have a list of towns I want to display ordered according to their size (population). That's quite easy: in the feature properties of the table just select main properties and than attributes. There, you can select which fields from the attribute table of a certain layer you want to show - even custom-generated fields with expressions are possible. Beneath, you can choose one of the fields defined above for the sort order. Perfect! It orders my towns in the correct order.

However, I would like to display the no. of inhabitants in a formatted style with a delimiter - thus 1.000.000 instead of 1000000. Works perfectly with the expression format_number(population,0). The problem: formatted this way, the sort-order seems to interpret the numbers as strings, thus ordering like this:

15,000

150,000

20,000

200,000

1,500,000

2,000,000

It is not possible to select something else for the sorting-order than the fields that are visible in the table - thus using an additional field (hidden in the table) for the definition of the sorting order is not possible.

I tried to create a virtual layer as explained here: https://gis.stackexchange.com/a/339276/88814.

However, that is not a solution in my case (see why here: QGIS: no geometry in virtual layer). So how can I implement a correct sort ordering that respects the desired formatted numbers?


texts: ArcGIS merge/unite two rasters with same extent with RasterCalculator

texts:
I've a problem with the ArcGIS Raster Calculator:

I have two 2-Bit Rasters with the same extent, the only data stored in it are 1-values and NoData.
Now I want to combine this two Rasters: If in layer one and/or two the value is 1, I want a 1 in my new raster. If the two layers are 0, it should be NoData.

I tried:

((Raster1 == 1) | (Raster2 == 1)) *1

Doesn't work, neither I don't know if this should be correct..

Con(((Raster1 == 1) | (Raster2 == 1)),1)

And the last one:

Con((Raster1 == 1), 1, Raster2)

All the time it returns me a new Raster which is only 1 if Raster1 AND Raster2 is 1. How should I do this? What's the mistake?


texts: Setting up a QGIS Plugin Development Environment in Visual Studio and PTVS

texts:
I would like to create a plugin using Visual Studio 2015 and Python tools (PTVS). I have a basic hello world plugin written outside visual Studio and It loads fine in QGIS 2.10 Pisa. I have installed the Visual Studio and PTVS and I am looking for a way to work using the Visual Studio Environment.

I have set up a Python environment within visual studio however I cannot get the REPL prompt to work. Has anyone done anything similar?


texts: How to create arrows lines automatically on ArcGIS, for a thematic map of directed flux?

texts:
I need to create a map with the Import and Export Flux between countries. I have the needed data in my table. But I wish to know if there is any tool or extension which will allow me to create those lines or arrows automatically (similar to how we do with labels, but those create words, which organize themselves in the Data Frame).

For that all I should do is to indicate lineament of origin and destination, representing the seaports.

I need to do something like this, but with arrows:


texts: I'm working with QGIS and I want to convert a line into many points, does someone know how?


texts: Arcpy crashes in loop geoprocessing (general behaviour)

texts:
I have several standalone scripts that create a list of shapefiles to process (or coverages in one of the scripts).  Sometimes my script crashes after it processes three files, sometimes ten, sometimes with just one to go. In some cases, it's quite complex operations, but sometimes, it's just a simple arcpy.Union_analysis(inShp, outShp) or with CalculateField_management().

It seems totally random because I just delete the outputs before restarting my scripts and it may or may not work the next time. Eventually it works, sometimes on the first try, but other times, I could have to retry it 2 or 5 or 10 times.

Any ideas what could cause this?
Should I clean or delete something from memory that could cause a crash?

Thanks!


texts: How to add popups to MapImageLayer (JavaScript API 4.0 esri)

texts:
I have to create a web map application with webservices (WMS or WFS).

I use the API for JavaScript of ESRI (4.0) and MapImageLayers for my data.
But I didn't find how to add popups (or "popupTemplate"= to a MapImageLayer.

I use this sample for my code (https://developers.arcgis.com/javascript/latest/sample-code/sandbox/sandbox.html?sample=layers-mapimagelayer-sublayers)


texts: When I try to draw a 2.5 mile buffer around a point in QGIS, the entire screen is filled with the buffer.  I have changed the projection of the shapefile and the project, as you can see from the screenshots.  This is using QGIS 3.8.0


texts: Extracting map properties (extent, metadata, coordinate system, transformation) using arcpy.mp?

texts:
I am trying to fetch map properties using arcpy.mp (since I am using ArcGIS Pro trial version) from mxd files.

I am able to get the extent info for individual layers but not for the dataframe. Sample code

fc = "\demo.gdb\dist"
spatial_ref = arcpy.Describe(fc).spatialReference
print("Co ordinate : ", spatial_ref.name)
extent = arcpy.Describe(fc).extent

There are code samples available with arcpy.mapping module (for desktop vsersion), but I am using ArcGIS Pro, which supports arcpy.mp module.


texts: Intersecting three polygon layers in ArcGIS?

texts:
I have three data sets of San Francisco, TerraSAR-X, ALOS-PALSAR, Radarsat2. I have extracted KML files using Envi Sarscape.

In order to find the intersection of these areas, I have converted these kml's into polygon shapefiles using:

  Conversion Tools ---> From KML ---> KML To Layer
  Conversion Tools ---> To Shapefile ---> Feature Class To Shapefile
  (multiple)

Each of those shapefiles have the following coordinate system when I used the properties from ArcCatalog

But when I used the following toolbox to extract the intersection shapefile

  Analysis Tools ---> Overlay ---> Intersect

The extracted shapefile will be as follows!!

What's the problem?

Edit to answer radouxju
Union would give the following result


texts: Disable greeting message from spatialite cmdline

texts:
I can emit plain values in SQLite3:
$ echo 'select 1;' | sqlite3
1

However, SpatiaLite always emits a greeting:
$ echo 'select asgeojson(st_point(1,2));' | spatialite
SpatiaLite version ..: 5.0.0    Supported Extensions:
        - 'VirtualShape'        [direct Shapefile access]
        - 'VirtualDbf'          [direct DBF access]
        - 'VirtualText'         [direct CSV/TXT access]
        - 'VirtualGeoJSON'              [direct GeoJSON access]
        - 'VirtualXL'           [direct XLS access]
        - 'VirtualNetwork'      [Dijkstra shortest path - obsolete]
        - 'RTree'               [Spatial Index - R*Tree]
        - 'MbrCache'            [Spatial Index - MBR cache]
        - 'VirtualFDO'          [FDO-OGR interoperability]
        - 'VirtualBBox'         [BoundingBox tables]
        - 'VirtualSpatialIndex' [R*Tree metahandler]
        - 'VirtualElementary'   [ElemGeoms metahandler]
        - 'VirtualRouting'      [Dijkstra shortest path - advanced]
        - 'VirtualKNN'  [K-Nearest Neighbors metahandler]
        - 'VirtualGPKG' [OGC GeoPackage interoperability]
        - 'SpatiaLite'          [Spatial SQL - OGC]
        - 'VirtualXPath'        [XML Path Language - XPath]
PROJ version ........: Rel. 7.2.0, November 1st, 2020
GEOS version ........: 3.8.1-CAPI-1.13.3
RTTOPO version ......: 1.1.0
TARGET CPU ..........: x86_64-linux-gnu
the SPATIAL_REF_SYS table already contains some row(s)
{"type":"Point","coordinates":[1,2]}
$

Both the greeting and my desired GeoJSON are going to stdout:
$ echo 'select asgeojson(st_point(1,2));' | spatialite >/dev/null
the SPATIAL_REF_SYS table already contains some row(s)
$

How can I run spatialite command-line tool and get it to emit only the stuff that I need (in this example, the GeoJSON)?


texts: How to display local file's geometries as rubber bands in QGIS

texts:
I'd like to get geometries from a local file (e.g. GML) displayed as rubber bands (QgsRubberBands) in QGIS. However, I wouldn't like to load the layer into the map because in the real scenario there could be a large number of files to get geometries from.

I've been attempting it with the following code (you can copy it directly to the QGIS Python Console), but the problem is that only the last feature is displayed:

from qgis.core import QgsFeature, QgsVectorLayer
from qgis.gui import QgsRubberBand
def loadGMLAsRubberBands(fileName):
    featList = []
    rb = QgsRubberBand(qgis.utils.iface.mapCanvas(), True)
    vlayer = QgsVectorLayer(fileName, "layer", "ogr")
    allAttrs = vlayer.dataProvider().attributeIndexes()
    vlayer.select(allAttrs)
    f = QgsFeature()
    while vlayer.dataProvider().nextFeature(f):
        featList.append(f)
        rb.addGeometry(featList[-1].geometry(), None)

loadGMLAsRubberBands("/tmp/gml/sample_three_features.gml")

A sample GML file can be found here.

I've already attempted several things with no success, namely to store layers (vlayer), features, geometries and/or rubber bands in a list (as suggested in some QGIS mailing-list threads).

I'm using QGIS 1.8.0-Lisboa (rev. a1255fc), Python v.2.7.2, GDAL/OGR v.1.7.3 on Ubuntu/Linux.

Any help would be appreciated.


texts: Getting an ordered list of ids of selected objects using PyQGIS

texts:
I am trying to get the ids of the objects of my selection. The list should be ordered from my first selected object to the last selected object.
I do not have much experience with using SIGNAL and currently, I am using selectionChanged.connect which prints out any changes in the selection and shows the orders I wish. However, I cannot assign these ids in order to further use them (the ids).
The user should ideally select the objects one by one and if multiple objects are selected in one action, it should throw an error.
Is there any way I could do this?


texts: Conflict with leaflet-markercluster-rails and another gem?

texts:
I'm wondering if anyone is aware of some sort of conflict between leaflet-markercluster-rails and another gem that would prevent the markers from displaying properly?

It's the only thing I can think after almost a week for messing around with code.  NOTHING seems to work.  See my original post for my code.

I'm running rails 4.2.5 on Windows box if that mattes.  ANY help would be wonderful. :)


texts: I need to draw Rectangle with Simple Fill Symbol

texts:
I need to draw this rectangle on layoutview :

I try this but it doesn't work. How can I draw this rectangle  ?

        IPoint p1 = new PointClass();

        p1.X = 9.9936; p1.Y = 5.8777;

        IPoint p2 = new PointClass();

        p2.X = p1.X; p2.Y = p1.Y + 2.62;

        IPoint p3 = new PointClass();

        p3.X = p1.X + 0.152 ; p3.Y = p1.Y;

        IPoint p4 = new PointClass();

        p4.X = p3.X; p4.Y = p3.Y + 2.62;

        IPointCollection ring = new RingClass();
        ring.AddPoint(p1);
        ring.AddPoint(p2);
        ring.AddPoint(p3);
        ring.AddPoint(p4);

        IRgbColor color = new RgbColorClass();
        color.Red = 255; color.Blue = 28; color.Green = 20;

        ISimpleFillSymbol simpleFillSymbol = new SimpleFillSymbolClass();
        simpleFillSymbol.Color = color;

        IElement element = new PolygonElement();
        element.Geometry = ring as IGeometry;  // I try this.
        element. ???  // I need to apply symbol on this element. but it hasn't symbol property

        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
        IPageLayout pageLayout = new PageLayout();
        pageLayout = mxdoc.PageLayout;
        IGraphicsContainer graphicsContainer = pageLayout as IGraphicsContainer;
        graphicsContainer.AddElement(element, 0);

        IActiveView activeView = pageLayout as IActiveView;

        activeView.Refresh();

This solution has problem at :

    IElement element = new PolygonElement();
    element.Geometry = ring as IGeometry;
    ((IFillShapeElement)element).Symbol = (IFillSymbol)simpleFillSymbol;

@Berend's solution is working..


texts: Points along geometry & QChainage unexpected result in QGIS (missing points)

texts:
I am using Windows 10 enterprise and QGIS version 3.10.1 A Coruna. I have found an unexpected result when I use QGIS to generate points every 160 meters along a road layer (i.e. points evenly spaced along each road in my layer). I find a pattern of what I believe to be aberrantly missing points. I expected the resulting points layer to have points spaced every 160 meters along each and every road.  What results instead is sections of road without points - while there are some points spaced 160 meters apart there are then large swaths (along the same roads) that have no points.

I have used both "points along geometry" tool in the processing toolbox as well as the QChainage plugin - both yield identical results.

I have tried re-projecting my layer with no change in results.  The current projection is EPSG 3857.

Here's a screenshot of the roads with patchy/missing points:


texts: Including metadata when importing features into GeoPackage using ArcGIS Desktop?

texts:
After much frustration reading vague documentation and getting obscure errors I managed to create a GeoPackage using ArcMap v10.6.1 and import features into it.

I was very happy until I discovered that no metadata is being imported along with the features. The GeoPackage is meant to be a modern replacement for shapefiles, but without metadata I can't use it.

Does anyone know how to get the metadata to transfer along with the features?

The procedure I use at the moment is:

Create an empty .gpkg file using ArcToolbox > Data Management Tools > Workspace > Create SQLite Database, and selecting Spatial Type = GEOPACKAGE. This only needs to be done once as a copy of the empty .gpkg file can be made outside of ArcGIS whenever one is needed.
Open the layer you want to convert to a GeoPackage
In the catalog window, navigate to the .gpkg file, right-click, and select Import > Feature Class (single)..., which opens a dialog box.
In the Input Features field, select the layer you want to convert. The Output Location field should already be populated with the still-empty .gpkg file. Enter an arbitrary name in the Output Feature Class field. Then click OK. A new window should appear showing the progress of the copy operation.
If you get an unhelpful error saying it didn't succeed, check that you have at least standard (ArcEditor) licence activated. It doesn't work with the basic (ArcView) licence.
Click the Close button once the operation is complete.
Note that ArcMap keeps a lock on the .gpkg file meaning you can't move or do many other things with it until you close the ArcGIS application completely.

Now when I right-click the freshly created layer under the .gpgk file and select Item Description, I get an error saying "The item's XML contains errors" where the metadata should be.


texts: Create a raster from file geodatabase featureclass based on a field from a related table

texts:
I have a USDA Statsgo soil polygon featureclass stored in a file geodatabase along with an attribute table; I have a one-to-one relationship class defined based on the soil type field. When I bring the featureclass into ArcMap I can view the attributes in the related table when I use the identity tool.

What I would like to do is create a raster from this featureclass using values from the related table. This seems simple enough, but when I used the Polygon to Raster conversion tool  the only fields available to use are those which are stored in the featureclass itself. If I define a join on the same featureclass and table I can then select the appropriate field in the joined table, which makes me wonder why I bothered creating the relationship class in the first place (resulting raster values are numeric - I lose my classification codes of A through D).

Can anyone explain 1) why I can't use the relationship class to access fields in the related table, and 2) how to keep the alphabetic soil classification codes?


texts: Selecting touching polygons by PostGIS for QGIS atlas?

texts:
I‘ve a PostGIS polygon layer where each feature represents a page of an atlas in the QGIS print composer.

CREATE TABLE public.composer (
gid serial NOT NULL,
geom geometry(Polygon, 31256),
page_number varchar(8),
page_north varchar(8),
page_west varchar(8),
page_south varchar(8),
page_east varchar(8),
CONSTRAINT composer_pkey PRIMARY KEY (gid));
;

INSERT INTO public.composer (geom, page_number) VALUES
(ST_GeometryFromText('Polygon((0 0, 300 0, 300 200, 0 200, 0 0))', 31256), 'A-A'),
(ST_GeometryFromText('Polygon((300 0, 600 0, 600 200, 300 200, 300 0))', 31256), 'A-B'),
(ST_GeometryFromText('Polygon((0 0, 300 0, 300 -200, 0 -200, 0 0))', 31256), 'A-C'),
(ST_GeometryFromText('Polygon((300 0, 600 0, 600 -200, 300 -200, 300 0))', 31256), 'A-D')
;

What I want to do is insert the numbers of the touching features into the attributes of each feature.

Can anyone help me with the SQL query?

Doing this manually is quite boring since my atlas contains around 100 features.


texts: How to edit AutoCad coordinate drawing boundaries

texts: How can I edit the Autodesk-AutoCad x and y coordinate surface? I would like it to have a larger coordinate boundary than the default...


texts: Forcing decimal numbers to use comma in ArcGIS Field Calculator?

texts:
Can I add the comma "," to force the decimals to remove from the attribute and replace with the number?

For example, I have an Acres field column that says 506320.0750 and I want to replace the decimal with the comma, say 506,320? I want that number to be permanment.

Is that possible? I'm using ArcGIS Desktop 9.3.1


texts: How to bind to coordinates , in sld GeoServer?

texts:
I have coordinates and there is an image (link to png image), I need to bind it to the coordinates (onlineresource tag) in the GeoServer. An example of how to google maps here https://developers.google.com/maps/documentation/javascript/examples/overlay-simple

    <ExternalGraphic>
        <OnlineResource xlink:type="simple"
            xlink:href="http://mysite.com/tn_${strToLowerCase(STATE_ABBR)}.jpg">
        <Format>image/jpeg</Format>
    </ExternalGraphic>

I want OnlineResource set as geo-referenced


texts: qgis snaps to outside instead of boundary

texts:
I have polygon layers that I need to be covered by other polygon layers, and many of my inner layers have a vertex on the border of an outer layer or in some cases they share an edge.  Both ST_Covers from Postgis and the topology checker ('must not overlap') inside Qgis confirm that some vertices are snapping outside of the outer polygons instead of on the boundary.

Is there any way to ensure that the inner polygons snap to the boundary instead of outside the outer polygons?


texts: Merging rasters using QGIS raster calculator?

texts:
I have 102 .tif files to merge. I want to use the raster calculator on these, but it said there was not enough memory when I tried using the raster calculator on the virtual raster of all 102 files. When I try raster -> misc -> merge, there is an error message saying:

  Error 5: Pixel type must be float32

.

Any recommendations?


texts: Mosaic rasters with same suffix (python)

texts:
I have a lot of rasters which have next names structure:

2017-08-02_1011_image.tiff
2017-08-13_5erf_image.tiff
2017-05-02_1555_image.tiff
2017-08-02_1001_image.tiff
2017-08-02_1001_image.tiff
2017-08-06_abc3_image.tiff
2017-08-06_45ei_image.tiff
2017-08-06_0005_image.tiff

I want to mosaic files with same suffix in names (4 symbols before «_image») and set resulted mosaics names with this four symbols – like «1001_mosaic».
I try this code:

import arcpy
from arcpy import env
from arcpy.sa import*

env.workspace = "H:\test\Image_names"
arcpy.env.overwriteOutput = True
print "Processing mosaic"

for i in range(1,5):
   ImgList = arcpy.ListRasters("*"+str(i)+"_image.tiff")
   arcpy.MosaicToNewRaster_management(ImgList, env.workspace, "image"+str(i)+"_mosaic.tif",  pixel_type="16_BIT_UNSIGNED", cellsize="", number_of_bands="4", mosaic_method="", mosaic_colormap_mode="MATCH")
   print "Completed!"

But it mosaics all files with same single symbol before «_image» from 1 to 5:2017-08-06_0005_ image.tiff with 2017-05-02_1555_ image.tiff; 2017-08-02_1011_image.tiff with 2017-08-02_1001_ image.tiff and so on.

I’m new with python.


texts:
I have a series of tram stops on a Carto map and would like to indicate the user's location.

Is this doable in Carto?

https://jamieparker.carto.com/builder/161cab96-859b-433d-8a72-5ae5653b1b8e/embed

I've played around unsuccessfully Does Cartodb have a function to change center location based on user location?

Probably user error on my behalf, although Carto has changed considerably since that post.


texts:
I am using Jquery Google Map V3 with markerclusterplus. I want to animate a cluster icon when hover the icon. My Clusterhover event code is,

var mc = new MarkerClusterer(map, markersarray,mcOptions);
    google.maps.event.addListener(mc, "mouseover", function (c) {
          console.log(this);
          $(this).effect( 'bounce', { times: 3 }, 'slow');
          //$("div.cluster").effect( 'bounce', { times: 3 }, 'slow');

    });

But I am facing following error,

Uncaught TypeError: $(...).effect is not a function(anonymous function) @ groupmap.html:309S.trigger @ main.js:20(anonymous function) @ markercluster.js:186

How do I achieve this?


texts: Rescaling NDVI from 0-255 to -1-1 in ArcGIS Pro?

texts:
My NDVI raster ranges from 0-255. My goal is to get values from -1 to 1.

If using the Raster Calculator, what would be the correct equation to achieve this?

I'm working in ArcGIS Pro.


texts: Layer 'save as' file displaying wrong/original name in TOC

texts:
I'm running ArcMap 10.3 and trying to create a new copy of a layer in my map document, which I can do, but when I load the new copy back into my map, the TOC displays its name as the old file from which I copied it, not the new filename.

For example, say I want to create a new version of Layer X that's currently in my map and call the new version Layer Y. I simply right click Layer X, select "Save As Layer File", specify folder destination, and click Save. When I click "Add Data" to bring it into my map, the Add Data dialog box displays the correct/new name (Layer Y in this case), but when I select it to bring it into my map, my TOC will now display 2 Layer X's (the original + the new one I just brought in) instead of the original Layer X and the new Layer Y. I have successfully done this before, so I know it's possible - especially given that it's not a particularly advanced geoprocessing step.

Just curious if anyone else has run into this and how I can get around it?


texts: Does anyone know some good tutorials to handle spatial data in R? I know that there are some packages like rgdal and maptools but I haven't actually found usefull tutorials to get in touch with the basic functions (read, plot, change symbology, join data, etc.) and/or to get a broader overview what are the acutal capabilites of R for spatial analysis and manipulation of Spatial data. Also recomendations for books, papers, blogs, etc. would be warmly welcome.


texts:
I created a polyline feature and then I created a shapefile. Then overlayed another map and I wanted to analyse the distribution of the distance from the polyline, but the extents of the shapefile of the polyline are limited by the polyline lengths, cutting out part of the area I would be interested to (as you can see in the pic below).
Is it possible to edit the size of the shapefile, so I can cover all the area?


texts: GeoJSON data visible in browser but not showing on map

texts:
I am using OpenLayers, GeoServer and Shapefiles.

I have published the shapefile to GeoServer and under layer preview I can see JSON data in browser when I select GeoJSON option from the dropdown for that file.

Problem is I cannot get the data to display on the map in OpenLayers.

This is the link I followed

And here is my code:

<!doctype html>

    <title>OpenLayers example</title>

    <link rel="stylesheet" href="https://openlayers.org/en/v4.6.5/css/ol.css" type="text/css">

My Map

I also read somewhere that I need to load features from URL dynamically using AJAX calls. If so, how would I do that? I didn't find anything beginner friendly on the subject.


texts: getFeatures() delivers some fields more multiple times

texts:
I am extracting data from selected fields using .getFeatures() in PyQGIS (see code below). The result is correct, but for some reason python delivers some features multiple times (see below; Lamu power project three times and Vipingo twice...). Any ideas why?

vectorLayer = QgsVectorLayer('C:/AIRPOLIM-ES/buf2.gpkg', 'buf2', "ogr")
QgsProject.instance().addMapLayer(vectorLayer)

rasterLayer = popKE
QgsProject.instance().addMapLayer(rasterLayer)

zoneStat = QgsZonalStatistics (vectorLayer, rasterLayer , 'buf2_pop_', 1,QgsZonalStatistics.Sum)
zoneStat.calculateStatistics(None)

##extract results

for feature in vectorLayer.getFeatures():
    print(feature['Plant'],',',feature['buf2_pop_sum'])

Results:

Kaloleni Cement Plant , 22221449.12682629
Lamu Power Project , 22570960.33722942
Lamu Power Project , 22570960.33722942
Lamu Power Project , 22570960.33722942
Pokot power station , 40265891.87767924
Vipingo power station , 4560350.673526992
Vipingo power station , 4560350.673526992


texts: Creating sub sections (groups) of layers in QGIS?

texts:
I am working on a project regarding Fish species and fishing methods with many layers.

This creates confusion.

I would like to create layer sub sections (groups) like the following:

Fish Species

species 1
species 2
species 3

Fishing Methods

method 1
method 2
method 3

and so on.

Possibly even a 2nd level sub section (layer group).

I have no idea how to tackle this.

I'm using QGIS 2.4.

To clarify an image of the project is below.


texts: Viewing a point geojson layer as markers in cesium

texts:
I have a point geojson layer and i would like to add and view it using cesium. I tried with the below code

var dataSource1 = new Cesium.GeoJsonDataSource();
var roads = dataSource1.load('Apps/SampleData/layout&area.geojson');
roads.then(function(dataSource1) {
    viewer.dataSources.add(dataSource1);
    viewer.zoomTo(dataSource1);

    //Get the array of entities
    var entities1 = dataSource1.entities.values;

    var colorHash = {};
    for (var i = 0; i < entities1.length; i++) {
        //For each entity, create a random color based on the state name.
        //Some states have multiple entities, so we store the color in a
        //hash so that we use the same color for the entire state.
        var entity = entities1[i];
        var name = entity.Text;
        var color = colorHash[name];
        if (!color) {
            color = Cesium.Color.BLACK;
            colorHash[name] = color;
        }
        //Set the polygon material to our random color.
        entity.point.material = color;
        //Remove the outlines.
        entity.point.outline = false;

    }
}).otherwise(function(error){
    //Display any errrors encountered while loading.
    window.alert(error);
});

As a result, i get all the points viewed as markers. How should i make changes to view it as a point?


texts: Estimating Area of farm under cultivation using NDVI?

texts:
I have  shapefiles for more than 1000 farms.

I want to use Hydrological model to model the water balance in these farms. Each farm is devided in to seperate parts .

My problem is that the whole area of these farms  is not under cultivation and in each season farmers decide based on available water and their budget to put how much of their land under cultivation

In other words, I do not know what the percentage of farms is being used by each farmer.

I have estimated NDVI from landsat8. These NDVI values can help me to find the area of each farm which is being used by farmer.

How can i find the area of my shapefiles which  have specific NDVI? For example, if I have shapefile for 200 Hac farm , what area of this shapefile has NDVI between 0.1 and 1?


texts: Calculate value from an equation using Field Calculator Python/VB Script

texts:
I wanted to calculate a result value from 3 fields using an equation.

VBScript or Python codes are also accepted..:-)

the equation is RESULT= [(R^2/3)(S^1/2)]/n..


texts: Handling Errors with Esri Leaflet Dynamic Map Layers

texts: I have an application that is using ESRI for leaflet to pull in MapServices.  I can successfully pull in map services and render in my application.  However, sometimes the map services return an error (400, 500, etc).  When this occurs I can no longer click on other data I am showing on the map (geojson data coming from a MongoDB).  Is there a way to handle these errors so they do not interfere with the user experience?


texts:
I'm clipping several rasters with different polygons and I've realized that some polygons fall on top of areas with no date values of the raster, returning me as a result completely empty clipps with NAN values. Is there a way for my code not to write raster without data?
import rasterio
import os

root = 'F:/LCCD/tile/'
shpfiles='F:/LCCD/bf2.shp'
outputPath='F:/LCCD/output'

rasterfiles = [os.path.join(root, file) for root, dirs, files in os.walk(root)
               for file in files if os.path.splitext(file)[1] == '.tif' ]

i = 0

for raster_path in rasterfiles:
    output_path =outputPath + '/' + str(i) + '_crop' + ".tif"

            with rasterio.open(raster_path) as tile:
                out_image, out_transform = rasterio.mask.mask(tile, shpfiles, crop=True)
                out_meta = tile.meta

                out_meta.update({"driver": "GTiff",
                                 "height": out_image.shape[1],
                                 "width": out_image.shape[2],
                                 "transform": out_transform,
                                 "nodata":0})
                i=i+1
                with rasterio.open(output_path, "w", **out_meta) as dest:
                    dest.write(out_image)


texts: Trying to understand Otsu method code in GEE

texts:
I hope this question is OK, I know it's not usual kind of questions that are being posted here.

I am trying to apply Otsu method on image in GEE. For that I have used the code that can be found here :
https://medium.com/google-earth/otsus-method-for-image-segmentation-f5c48f405e

The code is good and I believe it works, but my problem is that I don't fully understand the code and I would like to understand it better.

I have understood how the method works by reading this tutorial:
http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html,
but still' the steps are not clear to me.

I put here the code with my questions regard what different steps do when each row that I don't understand marked with // above the row.

///////////////////////Otsu

// Compute the histogram of the Difference one band image.  The mean and variance are only FYI.
var histogram = difference.select('VH').reduceRegion({
  reducer: ee.Reducer.histogram(255, 2)
      .combine('mean', null, true)
      .combine('variance', null, true),
  geometry: geometry,
  scale: 30,
  bestEffort: true
});
print(histogram);

// Chart the histogram
print(Chart.image.histogram(difference.select('VH'), geometry, 30));

///////////Otsu
  // Return the DN that maximizes interclass variance in VH (in the region).
var otsu = function(histogram) {
  //does count here take the histogram and convert it to array?
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
  //what is bucketmeans? is it the mean of each bar in the histogram?
  var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
  //I know in array size will show the dimensions. here we want to number of rows in the array from hostogram?
  var size = means.length().get([0]);
  //this step: not sure what happened here. we sum all the rows? what is get[0] in the end?
  var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
  var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
  var mean = sum.divide(total);
  //and also not sure what indices is
  var indices = ee.List.sequence(1, size);

  // Compute between sum of squares, where each mean partitions the data.
  var bss = indices.map(function(i) {
    //hard to understand what is i because I don't know what indices is
    var aCounts = counts.slice(0, 0, i);
    var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
    var aMeans = means.slice(0, 0, i);
    var aMean = aMeans.multiply(aCounts)
        .reduce(ee.Reducer.sum(), [0]).get([0])
        .divide(aCount);
    var bCount = total.subtract(aCount);
    var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
    return aCount.multiply(aMean.subtract(mean).pow(2)).add(
           bCount.multiply(bMean.subtract(mean).pow(2)));
  });
    print(ui.Chart.array.values(ee.Array(bss), 0, means));

  // Return the mean value corresponding to the maximum BSS.
  return means.sort(bss).get([-1]);
};

var threshold = otsu(histogram.get('VH_histogram'));
//var threshold=25;
print('threshold', threshold);

My end goal: To understand the Otsu method code, what happens in each step- specifically the ones I asked inside the code.


texts: Speed up colouring of contour lines by creating new layer?

texts:
Related to my previous question: Colouring contour lines differently when they are inside polygons of another layer I would like to create a new layer based upon the result of the geometry generator linestring/multilinestring expression:
intersection(collect($geometry),aggregate('NZ_Snow_Polygons__Topo__1_50k__72b39e76_c3e8_4878_a811_8e3dc5ade0d4','collect',$geometry,intersects($geometry,geometry(@parent))))

That as the expression for a linestring/multilinestring geometry generator works but is very very slow. Hence I want to use the result for a new static layer and render like the contour layer - or at least that is my goal.
I've tried "Vector Selection", "Select by expression" with the expectation I can save the selection as a layer, but the expression above runs for over 24 hrs at which point my patience is exhausted and I kill it.
What is the correct way to create a new layer that is the collection of contour lines sections where they intersect the snow polygons? i.e. Gives me the results I get with the geometry generator with none of the constant regeneration and slowness.


texts: Cascade filtering with ArcGIS API for JavaScript?

texts:
I want to populate my dropdown with attributes of shapefile for example in one dropdown on loading name of states given then after selecting the specific state next drop down loads its counties and also on every select map zoom to that level.

<!DOCTYPE html>

  <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
  <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no">
  <title>Cost Estimation</title>
  <link rel="stylesheet" href="https://js.arcgis.com/3.20/dijit/themes/claro/claro.css">
  <link rel="stylesheet" href="https://js.arcgis.com/3.20/esri/css/esri.css">

District INFORMATION PORTAL

             Select Area from the List

I have added dropbox in the left panel now I want to populate my drop box with the feature layer attributes. I have explained earlier about the action.


texts: Running geoserver behind nginx from custom location

texts:
I was required to make GeoServer to be accessible from /beta/geoserver however I always get redirected back to /geoserver which is 404 of course.

Here is the current location rule:

location /beta/geoserver/ {
   proxy_set_header Host               $host:$server_port;
   proxy_set_header X-Forwarded-Host   $host;
   proxy_set_header X-Forwarded-Server $host;
   proxy_set_header X-Forwarded-For    $proxy_add_x_forwarded_for;
   proxy_pass http://127.0.0.1:8081/geoserver/;
}

And my observations:

/beta/geoserver/ redirects to /geoserver/web/ which is 404

/beta/geoserver/index.html redirects to /beta/geoserver/web/ which is fine.

I can try login there but I'm immediately redirected to /geoserver/web;jsessionid=... which is 404 again and I'm stuck.

The only thing I've configured is the Proxy Base URL but that's used by the services but not by the GeoServer AFAIK.

I'm neither a Tomcat nor a GeoServer expert so I guess the solution is one configuration file away from me but I haven't found a single clue of how to solve it.


texts: Unable to find all points in extent area using ArcGIS API for JavaScript?

texts:
I am trying to create a map where I can draw a selection area, and get info of all graphic pointers under it, I tried to modify ArcGis example but the problem now is that I can only able to select only last pointer and its data.

Below is my code, You can check the fiddle as well.

var gl;
var highlightSymbol;
var defaultSymbol;
var jsn = [{"Name" : "Wessam AlDahin","Email" : "email@address.com","City" : "Dubai","Latitude" : 23.408847304,"Longitude" : 54.41760682,}, {"Name" : "Danish Adeel","Email" : "email_2@address.com","City" : "Sharjah","Latitude" : 23.325462417,"Longitude" : 54.390645617,}, {"Name" : "Kashif Makhdoom","Email" : "email_3@address.com","City" : "Dubai","Latitude" : 23.27369932,"Longitude" : 54.285749593,}];

require([
    "esri/map",
    "esri/geometry/Point",
    "esri/symbols/SimpleMarkerSymbol",
    "esri/graphic",
    "esri/layers/GraphicsLayer",
    "esri/toolbars/draw",
    "esri/symbols/PictureMarkerSymbol",
    "esri/InfoTemplate",
    "dojo/domReady!"], function (Map, Point, SimpleMarkerSymbol, Graphic, GraphicsLayer, Draw, PictureMarkerSymbol, InfoTemplate) {
    function initToolbar(map) {
        var tb = new Draw(map);
        tb.on("draw-end", findPointsInExtent);
        tb.activate(Draw.EXTENT);
    }

    function findPointsInExtent(result) {
        var extent = result.geometry;
        var results = [];
        require(["dojo/_base/array"], function (array) {
            array.forEach(gl.graphics, function (graphic) {
                if (extent.contains(graphic.geometry)) {
                    graphic.setSymbol(highlightSymbol);
                    results.push(graphic.getContent());
                }
                //else if point was previously highlighted, reset its symbology
                else if (graphic.symbol == highlightSymbol) {
                    graphic.setSymbol(defaultSymbol);
                }
            });
        });

        //display number of points in extent
        dojo.byId("inextent").innerHTML = results.length;

        //display list of points in extent
        dojo.byId("results").innerHTML = "

" + results.join("") + "

";
    }

    map = new Map("map", {
            basemap : "streets",
            center : [54.285749593, 23.27369932],
            zoom : 10
        });
    map.on("load", function () {
        $.each(jsn, function (key, val) {
            gl = new GraphicsLayer();
            name = this.Name;
            email = this.Email;
            content = '

' + email + '

' + this.City + '

' + status + '

';
            pointVal = new Point(this.Longitude, this.Latitude);
            symbolVal = new PictureMarkerSymbol({
                "url" : 'https://webdesign.danols.com/static/template/images/icons/light/pin_map_icon&48.png',
                "height" : 20,
                "width" : 20
            });

            g = new Graphic(pointVal, symbolVal, {
                'title' : name,
                'content' : content
            }, new esri.InfoTemplate('${title}', '${content}'));
            gl.add(g);

            map.addLayer(gl);

        });

        initToolbar(this);
    });
    highlightSymbol = new SimpleMarkerSymbol().setColor(new dojo.Color([255, 0, 0]));
    defaultSymbol = new PictureMarkerSymbol({"url" : 'https://webdesign.danols.com/static/template/images/icons/light/pin_map_icon&48.png',"height" : 20,"width" : 20});
});


texts:
I want to make style like:

This is my style for start point:
<PointSymbolizer>
    <Geometry>
        <ogc:Function name="startPoint">
            <ogc:PropertyName>geom</ogc:PropertyName>
        </ogc:Function>
    </Geometry>
    <Graphic>
        <Mark>
            <WellKnownName>shape://oarrow</WellKnownName>
            <Stroke>
                <CssParameter name="stroke">
                    <ogc:Literal>#FF0000</ogc:Literal>
                </CssParameter>
                <CssParameter name="stroke-linecap">
                    <ogc:Literal>butt</ogc:Literal>
                </CssParameter>
                <CssParameter name="stroke-linejoin">
                    <ogc:Literal>miter</ogc:Literal>
                </CssParameter>
                <CssParameter name="stroke-opacity">
                    <ogc:Literal>1.0</ogc:Literal>
                </CssParameter>
                <CssParameter name="stroke-width">
                    <ogc:Literal>2</ogc:Literal>
                </CssParameter>
            </Stroke>
        </Mark>
        <Size>30</Size>
        <Rotation>
            <ogc:Function name="startAngle">
                <ogc:PropertyName>geom</ogc:PropertyName>
            </ogc:Function>
        </Rotation>
    </Graphic>
</PointSymbolizer>

And result:

How can i flip start mark?


texts: Qgis Profile Tool plugin with wrong units

texts: I am using qgis 2.18. I installed the plugin Profile Tool and i am trying to do a terrain profile from a DEM. I am working in projection WGS 84. The profile i made with that plugin doesn't have the proper lenght units. The vertical scale is right, in meters, but the horizontal one is wrong. Insted of the numbers of meters, it appears..0,01...0,02....Does anyone knows how to fis this? When i use the qgis scale the units are ok. Only Profile Tool is wrong.


texts: Looking for a pure Python method to convert Lat/Long To MGRS

texts:
Does anyone know of any packages out there that implement a pure python method to convert lat/long coordinates to MGRS coordinates?

PyPi has a package, but it requires some c-extensions from another source, and I didn't want to use that.  I figured before I go an recreate a complicated process, I'm hoping the community can make my life a tad easier.

Thanks everyone!


texts:
I am just learning how to use QGIS (3.4.1) and am having some issues overlaying some data (CSV file) with a QuickMapServices base map. I feel like this is an issue with the CRS, but am not 100% sure.

I have set the base map from QuickMapServices outlined below:

Outlined in the bottom right pane, this is set to ESPG: 3857, however, the coordinate (-11315, 6710551) is way out of place when I hover over London UK etc.

I have then added my own data, outlined below:

Which is also in ESPG: 3857, but the coordinates seem to be the wrong way around?  I.e. the point -0.37, 51.36, should be 51.36, -0.37. But then even if it was it still wouldn't be lined up correctly with my base map.

I am pretty confident in my wich is formatted as below, but I'm just at a bit of a loss as I seem to have changed the projections for both layers to ESPG: 3857?


texts: Adding polyline in ArcMap using ArcObjects?

texts:
I am trying to add polyline in ArcMap by using tool for this I have written some code, but the polyline is not adding in the proper extent.

How can I change my code to draw a polyline?

public override void OnClick()
        {
            IMxDocument mxd = m_application.Document as IMxDocument;

            ESRI.ArcGIS.Carto.IActiveView activeView =  mxd.ActiveView;
            IScreenDisplay sd = activeView.ScreenDisplay;
            DrawPolyline(activeView);

        }
        private void DrawPolyline(ESRI.ArcGIS.Carto.IActiveView activeView)
        {

            if (activeView == null)
            {
                return;
            }
            ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;

            // Constant
            screenDisplay.StartDrawing(screenDisplay.hDC, (System.Int16)ESRI.ArcGIS.Display.esriScreenCache.esriNoScreenCache); // Explicit Cast
            ESRI.ArcGIS.Display.IRgbColor rgbColor = new ESRI.ArcGIS.Display.RgbColorClass();
            rgbColor.Red = 255;

            ESRI.ArcGIS.Display.IColor color = rgbColor; // Implicit Cast
            ESRI.ArcGIS.Display.ISimpleLineSymbol simpleLineSymbol = new ESRI.ArcGIS.Display.SimpleLineSymbolClass();
            simpleLineSymbol.Color = color;

            ESRI.ArcGIS.Display.ISymbol symbol = (ESRI.ArcGIS.Display.ISymbol)simpleLineSymbol; // Explicit Cast
            ESRI.ArcGIS.Display.IRubberBand rubberBand = new ESRI.ArcGIS.Display.RubberLineClass();
            ESRI.ArcGIS.Geometry.IGeometry geometry = rubberBand.TrackNew(screenDisplay, symbol);
            screenDisplay.SetSymbol(symbol);
            screenDisplay.DrawPolyline(geometry);
            screenDisplay.FinishDrawing();
        }


texts: Varying extent of each data driven page using ArcGIS Desktop?

texts:
I would like to export data driven pages, each with a different extent. I understand that within the Setup for the data driven pages, there is an option 'Extent' and one can choose 'Data Driven Scale'. I'd expect that if I select this and choose a column from the table of the data driven page layer, the extent would change according to the value. Apparently this doesn't happen.

How does this option work exactly?


texts: I would like to create a few subplots within my project, I don't really know how to do this in one layer I've tried creating polygons but that isn't really what I want I just want the boundary to outline the area I am focusing on. Is this possible?


texts: EVI coefficient values for Sentinel 2 image analysis

texts:
I am estimating EVI for Sentinel imagery and wanted to know the values of the coefficients L, C1, C2 etc.

All the literature I have found use values from MODIS.


texts: Patch over raster A with values from raster B where raster A values are nodata in GDAL

texts:
I'm trying to come up with a way to patch over raster A with data from raster B. Extent of raster A can be smaller than raster B.

More specifically, I'd like to patch over nodata values from raster A with raster B values as well as create pixels in raster A where they exist in raster B.

Simple gdal_merge doesn't work for me as this doesn't work:

gdal_merge -o out.tif rasterA.tif rasterB.tif

as it is replacing values in raster A with values of raster B and would like values in raster A to be kept where they are not nodata.

and this won't work

gdal_merge -o out.tif rasterB.tif rasterA.tif

as it will replace values in raster B with nodata from raster A.

gdal_calc won't work, as rasters have different dimesions.

I'm aware that this such conditional logic can be implemented in Arc, but I'm looking for a solution using GDAL. Any suggestions?


texts:
I am mapping the accessibility of job opportunities in Detroit for those who rely on public transportation.  To do this, I need to create a circular buffer surrounding each bus stop with a radius extending out a mile in every direction. One mile represents the average time it takes a person to walk 20 minutes, a value chosen based on empirical research that an average person walks 4.5 feet a second (or 3 miles per hour).

In ArcGIS 3.0, I only see a "Variable Distance Buffer" tool under Vector--> Geoprocessing Tools.  (See pic below).  But these fields are not intuitive to me, even after clicking the help button.  Where and how do I specify 1 mile for the radius measurement to extend out from the point in all directions?  P.S. I should add I do not yet know python, so looking for alternative solutions.


texts: Displaying data from csv (or other file format) onto GeoJSON popup?

texts:
I'm trying to display data from another file such as a csv into a GeoJSON popup window. So far I know I can match the data using a nested for loop with a common id between the two files, but I'm wondering if it's even possible to display this data in the GeoJSON popup window. This is all done with a leaflet map as well.

My reason for displaying data from another file is that it needs to be updated frequently which is easier to do on a csv.


texts: How to convert UTM projected map to slippy map tilenames

texts:
I work with a map of Paris in ETRS89 / UTM zone 31N projection (EPSG:25831) and need to write slippy map tiles in filenames like /z/x/y.png.

I made a custom python script to read the QGIS project file and write the tiles but I can't get the tile names right.

The extent of my map [minX, minY, maxX, maxY] is something like [230000, 6270000,  295000, 6225000]. So for example in zoom 10 they should be named something like:

/10/518/351.png
/10/519/351.png
/10/518/352.png
/10/519/352.png

How could I calculate the x and y tile names correctly? For now I try to calculate them starting with the transformation of the upper left corner from degrees to meters. The meter coordinates origin of slippy maps (85W 180N) in meters are:

originX = 471118.87
originY = 9446194.39

Now I calculate the tile size in meters using the resolution in meters (for example 76.43702828517624 for z=10) and tile size in pixels (for example 256):

tileSize = resolutionsInMeters[z] * tileSizeInPixels

With that numbers in mind I calculate the x and y tile names:

x = int( math.floor( ( minX-originX ) / tileSize ) )
y = int( math.floor( ( originY-maxY ) / tileSize ) )

The result definitely is wrong, any hints to get valid tile names?


texts: How to identify feature vertices that are part of a donut hole in ArcGIS 10?

texts:
I´m exporting the coordinates of polygon vertices to a CSV-File using the "Feature Vertices to Points" and then the "Export feature attributes to ASCII" tools.

Some of these polygons have holes (donuts). Can I do the analysis in a way, that in the resulting CSV-File it´s possible to identify if a point has been part of the outer border or the inner border (donut) of the polygon?


texts: Creating equirectangular projection from 3d spatial model?

texts:
Apologies if this is not the right group for this question, but I think there may be people here who can help. (I'm open to suggestions for better places to ask.)

I have a program that contains a model of a planet. The model contains the surface of the planet and objects that have height e.g. mountains. The program accepts the (3d) coordinates of a viewer and the direction the viewer is looking and computes what the viewer would see on a flat screen (using the rules of perspective).

What I would like it to do now is to compute the equirectangular 2:1 projection corresponding to the viewer's position - i.e. the full-sphere panorama at the viewpoint. The maths involved presents no difficulties to me but I want to avoid a pixel-by-pixel calculation of the projection when I am sure this has been done before.

How can I create an equirectangular projection from a 3d spatial model?

My search so far has been limited to what Google mostly turns up i.e. the panoramic photgraphers where the assumption is that your source data is captured on a camera that rotated about its null parallax point. I have considered simulating taking photos in this way (including the overlapping, sigh) and using a stitcher like Autopano, but this will not work because stitchers only find control points in images with "texture" which my images would not have. I suspect the people I need to talk to are game programers but I haven't found them yet.

My projections will be gigapixels in size so the implementation needs to be efficient.


texts: How to build tile server from Openstreetmap data for beginner?

texts: I am a beginner at GIS and I am developing a Java Application which needs a map. (The map would specifically be of just one city in the Philippines.) I have researched on different ways to do this and I have found that I will need a tile server to put the map because I need to be able to access this without the connection of the internet. I am a beginner and most of the tutorials I have seen have been complicated and have steps that I do not know how to do yet. Does anyone know of any tutorials or can help me start creating the tile server considering that I am a beginner at this?


texts: Solving Geoserver - No service: ( ows ) - exception

texts:
I have an application running on PHP which makes an AJAX request to a OWS service.

The OWS service is provided by GeoServer running on Tomcat and protected with a Nginx reverse proxy

The condition specified in nginx is:

location ~ /geoserver/MY_WORKSPACE(.*) {

    valid_referers server_name mydomain1 mydomain2;

    if ($invalid_referer) {
        return 403;
    }

    proxy_set_header "attribute" "value";
    proxy_pass http://localhost:81/geoserver/MY_WORKSPACE$1;
}

Based on the referrer (domains allowed) the Ngnix reverse proxy lets the request pass or not.

The request is a full ows request https://mydomain2.com/geoserver/MY_WORKSPACE/ows?service=WFS&version=1.0.0&request=getFeature&typeName=MY_LAYER&maxFeatures=3000&transparent=true&outputFormat=application%2Fjson

Without the condition in the reverse proxy the request gets back the json data flawesly.

With the condition instead of getting back my JSON file what I get from GeoServer is a XML document (status code 200) with the following information:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:ows="http://www.opengis.net/ows" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows https://mydomain2.com/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
  <ows:Exception exceptionCode="InvalidParameterValue" locator="service">
    <ows:ExceptionText>No service: ( ows )</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts: PackageMap script failed to execute in IDLE but runs in ArcMap Python window?

texts:
This runs in the ArcGIS Python window but not from the IDLE command prompt.
import arcpy
import os
from arcpy import env

env.overwriteOutput = True
arcpy.env.workspace = "D:/Production/MobileData/FireAssets.gdb"

mxdLocation = "D:/Production/MobileData/MobileAssets.mxd"
tempFile = "D:/Production/MobileData/FireMobile_Assets.mpk"

arcpy.PackageMap_management(in_map=mxdLocation, output_file=tempFile, convert_data="CONVERT", convert_arcsde_data="CONVERT_ARCSDE", extent="MAXOF", arcgisruntime="RUNTIME", version="10.2")

This works fine in my ArcMap 10.7 Python Window but fails in IDLE or Python.
This is the error I get:

enter code here
    Traceback (most recent call last):
  File "D:\Production\MobileData\TEST.py", line 12, in <module>
    arcpy.PackageMap_management(in_map=mxdLocation, output_file=tempFile, convert_data="CONVERT", convert_arcsde_data="CONVERT_ARCSDE", extent="MAXOF", arcgisruntime="RUNTIME", version="10.2")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.7\ArcPy\arcpy\management.py", line 8884, in PackageMap
    raise e
ExecuteError: ERROR 999999: Error executing function.
Failed to execute (PackageMap).

This map packages layers from Enterprise geodatabase, FGDB, and a shapefile and exports to a 10.2 runtime map package for our mobile map.  The basic script ran fine for years on ArcMap 10.3 - then one day it did not, and we are not sure what update jacked it up.
I ended up upgrading the machine and upgrading to ArcMap 10.7. The script runs in ArcMap python window, but not standalone.
I ran analyzer and corrected all errors for descriptors, tags, etc.  I made sure all the layers were using the same connection reference using the find and replace sde connection command.  I rebuilt the FGDB for 10.7.  I ran check geometry and repaired null values on one layer. All common solutions for PackageMap general error 999999.
I checked the PATH variables C:\Python27\ArcGISx6410.7.
I added a set PATH for C:\Python27\ArcGIS10.7 and removed the 64bit path but that was not it.  IDLE still launches in 64bit.
I read somewhere that PackageMap will not work with 64-bit background processing and that you need to call the 32-bit version of Python.  I tried changing the PATH to C:\Python27\ArcGIS10.7 but IDLE is still running in 64-bit.  I tried to find the executable for 32 bit IDLE but had no success.  Still researching the correct executable for IDLE 32.
I have learned how to call the script from a DOS CMD prompt and that basically worked by changing my directory to location of 32bit Python and then running myscript.py.
I still want to use 32bit IDLE to test and debug my scripts.

Can IDLE run in 32-bit after installing a 64bit version?


texts: Fixing broken gdal-tools plugin in QGIS?

texts: I have QGIS installed on my 2011 macbook pro and I recently upgraded my OS to 10.9. After the upgrade, I noticed that QGIS gives me this error (see attachment 1) saying that the plugin is not installed (but it is; the second screenshot shows that it is "local only"). I suppose I can re-install the plugin but I'm hoping that it isn't necessary. Can someone provide some guidance as to how to get the plugin working again? Is it going to require that I completely remove the plugin and re-install?


texts: v.generalize does not affect some boundaries

texts:
All in all, I would describe myself as a QGIS (3.12.2) beginner and have now contact with GRASS (7.8.2) for the first time.

I want to simplify the boundaries of my polygons considering the topology. I came across v.generalize and did it with the following settings:

layer=1 type=area type=boundary method=douglas threshold=4000 look_ahead=7 reduction=50 slide=0.5 angle_thresh=3 degree_thresh=0 closeness_thresh=0 betweeness_thresh=0 alpha=1.0 beta=1.0 iterations=1 output=generalized

It actually looks exactly the way I want it to, except that some borders have been left out of the algorithm. I get the message: "8 boundaries were not changed because changes would damage the topology"

It doesn't seem to me that fundamental topology problems of my polygons are responsible for the error, because depending on the threshold value, more or less (and sometimes also other) boundaries are affected. But i could be wrong.

Could someone give me a hint where the problem could be and how I can solve it?

My ShapeFile via Dropbox: ShapeFile


texts: Error whilst creating spectral graph on composite image in Google Earth Engine

texts:
I'm trying to create a spectral graph (https://developers.google.com/earth-engine/charts_image_regions)based off a composite image. When i use this code, multiple errors come up, that all end with pattern pixel_qa not not match any bands.

Does someone know how to overcome this issue?

Here is the link to the map
https://code.earthengine.google.com/e4e45f74b30028629cdbb6c95f0da9de


texts:
Suppose I have a number of .gpx files (these contain GPX waypoint data from a Garmin eTrex).  I want to load them into R with different names and manipulate them.

I can read one file thus:

library(maptools)
gpx.raw <- readGPS(i = "gpx", f = "file1_w_12_f_ddf.gpx", type="w")

Suppose I want to read a number of them into memory.  I could try a for loop:

files <- list.files(".",pattern = "*.gpx")
for(x in files){

    #Create new file name
    temp <- strsplit(x,"_",fixed=TRUE)
    visit.id <- sapply(temp,FUN=function(x){paste(x[1],x[4],substr(x[5],1,3),sep="_")})

    #read file with new filename
    assign(visit.id, readGPS(i = "gpx", f = x, type="w"))
}

Running above program yields following error:

  Error in read.table(con <- textConnection(gpsdata), fill = TRUE, ...) :
    no lines available in input
  In addition: Warning message:
  running command 'C:\PROGRA~2\GPSBabel\gpsbabel.exe -w -i gpx -f file1_w_12_f_ddf.gpx -o tabsep -F -' had status 1

Note that I was able to read this file on its own, so it would seem it has nothing to do with the file itself but with running readGPS in a loop.

In general I still find it very confusing how R treats variables like x above.  I am not sure how to modify the argument to readGPS from the stand alone instance f = "file1_w_12_f_ddf.gpx":  Should it be x, or f = x, or f = "x", or what?  Or maybe the problem is in the call to GPSBabel...

I include a sample file below so you can copy it to text editor, and save as .gpx. twice with different names and try yourself.

<?xml version="1.0" encoding="UTF-8"?>
<gpx
 version="1.0"
 creator="GPSBabel - http://www.gpsbabel.org"
 xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
 xmlns="http://www.topografix.com/GPX/1/0"
 xsi:schemaLocation="http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd">
 <time>2010-09-14T18:35:43Z</time>
 <bounds minlat="18.149888897" minlon="-96.747799935" maxlat="50.982883293" maxlon="121.640266674"/>
<wpt lat="38.855549991" lon="-94.799016668">
<ele>325.049072</ele>
 <name>GARMIN</name>
 <cmt>GARMIN</cmt>
 <desc>GARMIN</desc>
 <sym>Flag</sym>
 </wpt>
 <wpt lat="50.982883293" lon="-1.463899976">
 <ele>35.934692</ele>
 <name>GRMEUR</name>
 <cmt>GRMEUR</cmt>
 <desc>GRMEUR</desc>
 <sym>Flag</sym>
 </wpt>
 <wpt lat="25.061783362" lon="121.640266674">
 <ele>38.097656</ele>
 <name>GRMTWN</name>
 <cmt>GRMTWN</cmt>
 <desc>GRMTWN</desc>
 <sym>Flag</sym>
 </wpt>
 </gpx>


texts: Google Earth Issue Adding Overlay - Unknown error on image ingest

texts:
I'm trying to add a .tif into Google Earth Pro as a superoverlay. The file is a large island survey I georeferenced in ArcGIS, and I'm trying to preserve its high quality and coordinates.

When I try to open the file, I get the error "Unknown error on image ingest" and the file shows up under Places as a camera icon image, but even then there's no image to be seen in the view window.

A couple people on the internet have had this happen before, but they didn't receive answers. Do you have any advice?


texts: KMZ with local image reference not working in Google Earth web viewer

texts:
I'm developing an iOS geospatial application that can export to KMZ.  Geotagged photos are bundled along side the KML into a zipped KMZ.  Here is a snipped of the placemarks description tag:

<![CDATA[
    <img src="myphoto.png" />
]]>

When loading the KMZ into the online Google Earth it shows as a failed image but it does work on the Google Earth Pro desktop viewer.

I've been through the docs and tried different flavours of local references but nothing has worked.  Frustratingly, the documentation for KMZ has this sample file which also doesn't work in the web viewer.


texts:
Although this question has been asked before, my case seems to be different. If I have a raster in QGIS that has it's CRS set to EPSG:26915 (NAD83 / UTM zone 15N) and I reclassify it using r.reclass, the output raster has a CRS of EPSG:3159 (NAD83(CSRS) / UTM zone 15N).

This happens regardless of what my project CRS is set to or what the default CRS is set to in "Settings/Options/CRS". If, however, I warp my original raster to another CRS (e.g., EPSG:32615 WGS84 / UTM zone 15 N), r.reclass of the warped layer yields a reclassified output that is in the same CRS as the warped layer.

I understand that the difference between NAD83 and NAD83(CSRS) is small; however, for some programs (e.g., R) that require layers to have the same projection before analysis, this behavior can be troublesome. Is there any way to set QGIS or GRASS up to behave differently?

I am running QGIS 3.12.0 and GRASS 7.6.0 on a Debian 10 machine.


texts: Constraining Voronoi polygons to other layer using FME?

texts:
Using FME 2017.0.1.1 (Build 17291) is there any way to constrain the extent of created Voronoi polygons to polygonal features of another layer? Instead of clipping the Voronoi polygons after creating them.

Example below, I need Voronoi polygons to the extent of the green background polygons. Not simply clipped to the extent of the entirety of the background shapefile (Layer B).

If clipped afterwards, some of the data from points within the green polygons may overlap with the red area. If the red area had its own points (and Voronoi polygons) some of that area would be lost to points in the green.

These will be Voronoi polygons inside already created Voronoi polygons, so the background data set is continuous and all polygons will most likely contain one or more points.

~~This is a section of my previous workbench that created the background Voronoi diagram. As you can see, it is clipped at the end, which is fine as I only wanted that layer to be clipped to an outer boundary. Now I would like this next layer to be constrained by many features.~~

UPDATE: After using a PK from the polygons (Layer B) as a foreign key on the points (Layer A) and using the "Group By" parameter on the VoronoiDiagrammer, I get this error message in the Translation Log.

  VoronoiDiagrammer(SurfaceModelFactory): Fail to set up the Surface
  Model. A minimum of 3 points are required for the Surface Model

Many of my 'Layer B' polygons only contain 2 points from 'Layer A'.

If I clone all the points, can I randomly offset the clones before the VoronoiDiagrammer?

UPDATE 2:

SpatialFilter the polygons (Layer A) to remove any which do not
contain any points (Layer B) before giving points the foreign key.
Create 2 clones of all the points and offset the 2 sets of clones to
different constants. (Creating very small triangles). This removes the error of:

  Surface Model cannot be constructed. Either all points are linear or initial set of points did not consist of 3 non NAN elevations

If all clones were offset using the same shift in coordinates.

To get the workbench to run all the way through, I have to set the offset to a minute 1mm from the original x,y ensuring that all 3 created points fall within the polygon of 'Layer A'. The last problem is that if a polygon from 'Layer A' only contains one original point form 'Layer B', the Voronoi process only outputs the tiny triangles and does not fill to the extent of the polygon in 'Layer A'.


texts: Alternatives to ArcGIS Extensions using ArcObjects?

texts:
I was wondering if there are any projects based around using ArcObjects and external libraries to create open source alternatives to extensions in ArcGIS.

For example 3D Analyst, Spatial Analyst etc or even stuff like implementing some ArcInfo functionality in ArcView?


texts: Adding ArcGIS Desktop map in WordPress web page?

texts:
I have designed a map in ArcGIS Desktop but I do not know how can I add it to a WordPress web site.

How can I do this?


texts: Mosaicing thousands of raster files in ArcGIS?

texts:
I'm trying to mosaic a very large number of rasters (from 1000 to a maximum of 30000, each one having a very small size, less than 100 kb) into a unique final raster.

I have tried to use the Mosaic To New Raster tool available at:

ArcToolbox > Raster > Raster Dataset

I encounter some issues when I choose the input files from the input folder:

because I need to wait a huge amount of time before viewing the files loaded into the main dialog (I know this is certainly related to my PC performances). Even if, in some cases, I'm able to load all the input files, the duration of the analysis is very long due to the number of the files.

From some comments on this question, it seems that it should be avoided the writing of raster data to a personal geodatabase. I set the Default.gdb for the Output Location parameter, but I don't understand very well its role.

Since I'm unfamiliar with ArcGIS, is there a more efficient way for the mosaicing of a large number of input rasters? Any workaround or method for avoiding the overloading of the PC resources?


texts: Error 001359 when publishing map service to ArcGIS server?

texts:
I encountered this error when trying to publish a map service to an ArcGIS Server 10.1 (installed on linux).

  --------------------------- Service Publishing Result
  --------------------------- Packaging succeeded, but publishing failed.

  ERROR 001359: Failed to connect to the server.

  Please consult the GP Results Window for enhanced error description

.

Other info:
The geodatabase is registered with the server. And it did work in the past.
We use a PostgreSQL gdb server.


texts:
I'm trying to mosaic some georeferenced images using gdalwarp. But if I change the resampling method the result is always the same, why?
I need to build a mosaic where the value of each pixel is the average of all overlapped pixels at that point.

For example, my images have a resolution of 0.05m and I want a mosaic with 1m  resolution so I try:

gdalwarp -tr 1 1 *.tif merged.tif

and:

gdalwarp -r average *.tif merged_average.tif

but the results are always the same 1m resolution raster. I also tried with min and max resampling methods, but even in this case the results are the same. How can this be possible? Maybe I did something wrong in the code?


texts: Create a horizontal legend in a split map

texts:
Do you know how to get a horizontal legend in a split map in GEE? I already have the legend, but in vertical. I would like to know if there is any way to transform the legend to a horizontal one.

Here is my code: https://code.earthengine.google.com/d0fb69c5abd1fdedff4b9c38d1926dae


texts: Positioning size labels on lines that delimit the polygon in OpenLayers

texts:
I want to position the label that defines the size of the side of the polygon exactly parallel to the line and not fall on the line above it ...
Here is a picture as it is now:

Now it is exactly on the line and after the orientation of the line in the middle also in the middle I just want it not to have contact with the line but to be inside or outside the polygon parallel to the line.
Here is the code I have managed to develop so far:
var line_stile = new ol.style.Style({
    fill: new ol.style.Fill({
        color: 'black'
    }),
    stroke: new ol.style.Stroke({
        color: 'blue',
        width: 10
    }),
    text: new ol.style.Text({
        font: 'bold 11px "Open Sans", "Arial Unicode MS", "sans-serif"',
        placement: 'center',
        fill: new ol.style.Fill({
            color: 'black',
            width: 42
        })
    })
});

var layerLines = new ol.layer.Vector({
    source: new ol.source.Vector({
        features: [new ol.Feature({
            geometry: new ol.geom.LineString(coordinates),
            name: 'line'
        })]
    }),
});

dist = dist.toFixed(2);
var data = dist + " m";

line_stile.getText().setText(data);
layerLines.setStyle(line_stile);

map.addLayer(layerLines);


texts: How to implement value-by-alpha map in GIS?

texts:
I've been doing some background reading about cartograms recently and stumbled upon interesting alternative - value-by-alpha maps.

Concept is described in more details here and full paper is available [not openly tho] [here]2.

There are some details about the implementation on indiemaps blog:

  ... all of the above graphics were produced with ActionScript 3, using just a text editor and the latest free Flex SDK...

How could you implement something like that in GIS package (ArcGIS? QuantumGIS?)? Or R and ggplot2? Or Python?


texts: Creating depressionless DEMs with rgrass7

texts:
I'm working with the rgrass7 package to use GRASS 7.6.0 functions within R. What I'm trying to do is the following: Import a raster DEM with a custom CRS, fill the sinks and export the depressionless raster - all via RStudio.
Some sample data can be found at: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT

I have never worked with GRASS before so I might just be missing a single command, but all my outputs are just empty rasters. Maybe this is a projection problem?

My code so far is this:

library(rgrass7)
initGRASS(gisBase = "/usr/.../grass76/",
          home = tempdir(),
          mapset = "PERMANENT",
          override = TRUE)

# modify current mapset with custom projection
execGRASS("g.proj", flags = "c", proj4 = "+proj=aea +lat_1=25 +lat_2=50 +lat_0=37 +lon_0=87 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")
execGRASS("g.mapset", flags = "c", mapset = "new_mapset")

# load DEM
execGRASS("r.in.gdal", input = "/path/to/.tif", output = "GRASS_raster", flags = c("overwrite"))

# set region (following the comment of @dmci = fixes the issue)
execGRASS("g.region", raster = raster("/path/to/.tif"))

# fill sinks
execGRASS("r.fill.dir", input = "GRASS_raster", output = "GRASS_filled_DEM", direction = "flow_test") # fill sinks
# export depressionless DEM
execGRASS("r.out.gdal", input = "GRASS_filled_DEM", output = "path/to/filled_DEM.tif")) # export filled DEM

This process does create a GRASS_raster and a filled_DEM.tif but they are just empty 1x1 raster files when I load them into RStudio.

I have tried to do this process with the GUI of GRASS. The imported DEM can be seen in the map display window, but the results of r.fill.dir are only empty rasters like this:


texts: Multicolor inside borders for polygons in QGIS 3

texts:
I was reading a report and I came across a figure that I quite liked the look of. I found that the author uses a black border for all of their polygons coupled with inside colored lines which just hug the outside borders. I have been trying to approximate this effect in QGIS 3.16 but have thus far not been successful.
Does anyone know if this effect can be achieved via QGIS?


texts: Opening Lowrance .at5 files in Mac/Linux/Windows

texts: I have some bathymetry data in Lowrance .at5 format. I would like to read those data to my computer and process it (convert it into point cloud etc). Is there any software available to do so? Or is there any Python package to read those data? Or any program?


texts: Defining spatial reference using ArcGIS Pro?

texts:
I want to define with ArcGis Pro the right spatial reference for a layer, which has no coordinate system asigned, but coordinates - 300.040,40940 Top and 424.407,907900 Right.
I am not really certain, how can I recognise with the given coordinates the real projection or CS.

Do you have any idea, or does my question make sense?


texts: How to Add Field Between Two Existing Field in Shapefile

texts:
Can you please let me know if it is possible to add a new Field between two existing field using ArcPy?
For Example if I have "Name", "Area", fiels already in my shapefile can I add a new field called "Scope" between these two? so the table looks like "Name","Scope", "Area"

Thanks


texts: Python OGR Execute SQL - Get console output

texts:
I am currently trying to scan around 1500~ MapInfo TAB files for invalid geometry using Python and OGR.

So far, I've managed to run the following command to filter features that contain errors:

driver = ogr.GetDriverByName("MapInfo File")
datasource = driver.Open(file, 0)
result = datasource.ExecuteSQL("select * from '" + file + "' WHERE not ST_IsValid(geometry)", dialect = "SQLITE")

Interestingly, when the SQL query is running I get outputs to the terminal (despite not using a print command) as follows:

Warning 1: Self-intersection at or near point 419991.27111455891 422299.0516533262
Warning 1: Self-intersection at or near point 416196.84999999998 411811.95000000001
Warning 1: Ring Self-intersection at or near point 419726.59999999998 421968.40000000002

Is there anyway to get this console output to be exported somewhere so that I can see exactly where/what the errors are?

Ideally, I would like this information to be added to the filtered layer as a new field.


texts: Intersection of Shapefiles, clean borders

texts:
I'm trying to create a shapefile of the state of zip codes (in this example, I pick Virginia). However, the zip codes do not line up exactly to the state borders. So I run the following query to get a shapefile of the zip codes for the given state.

SELECT uszips.geoid10, uszips.geom FROM uszips,va_shp
WHERE ST_Intersection(uszips.geom, va_shp.geom)
GROUP BY uszips.geoid10, uszips.geom

However, when I checked the shapefile, it looks like this:

When the state of Virginia looks like this:

Is there any way I can return the shapefile with the zip code geom truncated to match the exact borders of the state?


texts: Is there a way to run a Def Query on an AOI?

texts: I want to spatially only display features within an AOI boundary.  Is there a way to do this like with the def query tab in a feature classes attributes?


texts: R: Efficient addition of random noise to a raster dataset

texts:
Question

What is the fastest way to assign random noise (from a truncated normal distribution) to an existing DEM in R

What I have now

I have a DEM to which I want to add random noise for a normal distribution. This process is to be iterated numerous times as part of a monte carlo-esque type scheme. I start by defining a function to give me such a constrained random value as described:

#' Get random value from a truncated normal distribution
#' a = min value
#' b = max value
#' m = mean
#' sd = standard deviation
library(truncnorm)

rand_noise<-function(a=0,b=10,m=5,sd=2){
    val=rtruncnorm(n=1,a=a, b=b, mean=m, sd=sd)
    return(val[1])
}

For this example let's say we have a DEM as follows:

library(raster)

# get data
dem <- getData('SRTM', lon=16, lat=48)

# reproject to utm zone 33N (the data is for Austria)
tsrs <- "+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
dem <- projectRaster(dem, crs = tsrs)

# aggregate it
dem=aggregate(dem, fact=5)

I'd like to create a second (matching) raster field with random values using my pre-defined function as fast as possible. Currently I'm converting to a dataframe and using lapply to assign the function to each xy pair:

# convert raster to df to assign noise
df=as.data.frame(dem, xy=TRUE)
df$noise=apply(df, 1, function(x) rand_noise())

keeps=c('x','y','noise')
df=df[keeps]
noise=rasterFromXYZ(df)
crs(noise) = crs(dem)
plot(noise)

I can then add the noise to the original DEM:

noisy=dem+noise

Is there a faster way to create and add the noise (as described above) without having to convert to a dataframe? Something akin to lapply? The outlined approach above will be slow for multiple iterations.


texts: I need to label the simple-features/shapes in a Well-known text file. Is there some best practice on how to achieve this? Should I just keep the labels as separate entities linked to the WKT or is it possible to add this to the WKT itself? Thanks.


texts: Loading geometry point tables in SQL Server Spatial

texts:
I need to load point data in SQL Server Spatial.
Similar to this post but I really don't know how to get XY data to be accessible by SQL
I want to do what was done in the initial part of this post (Identify points that plot inside the polygon):
How do I find points that are NOT within a list of polygons using MS SQL?
SELECT points.name, polygons.name
    FROM points CROSS JOIN
    polygons
    WHERE  (polygons.shape.STContains(points.shape) = 1)

How do you get the points loaded into your table and SQL to recognize the points with out having to define them like the below?
insert into @points values
(geometry::Point(1,1,0)),
(geometry::Point(2,3,0)),
(geometry::Point(5,4,0)),
(geometry::Point(4,2,0))

I am having a hard time defining points without doing it manually like above.
Ideally I could build a table like this by selecting the xy data I have in a table already.


texts: Generating random points inside polygons that capture values from multiple attribute fields

texts:
I have some census block data where one attribute field is a total count (let's say number of homes), and then additional attribute fields disaggregate counts in different categories (e.g. home value might have 5 bins: 50k-100k, 100k-150k, etc...). I'd like to create a shp file in QGIS that generates random points for each home in each block polygon but also retains the attribute name of the disaggregated counts as a value in a new attribute field. EG, each point would have a new categorical attribute called "home_value" with a value like 50-100k or 100k-150k.

I think I could do this by generating a different layer for each attribute using "random points inside polygons (variable)" and then merging the layers, but I'm hoping there is some easier way.


texts: gdal_calc and then gdal2tiles - scale colour and interpret nodata correctly

texts:
I have a GeoTiff file with 4 layers (red, green, blue, NIR).
The image has been clipped to a polygonal area, all data are set to 0 outside that area. It has Float32 as data type.

I want to calculate NDVI from this file and I want the pixels outside that particular area to be transparent.

As of now, I make a GeoTiff with NDVI:

gdal_calc.py --calc="(B - A) / (B + A) * (A > 0)" --outfile=ndvi.tif \
             -A source.tif --A_band 1 -B source.tif --B_band 4 \
             --NoDataValue=0

Then I run gdal_tiles.py on ndvi.tif to display the image on top of Google Maps, but the resulting image is not properly clipped (it's just grey outside of the field's borders, while I expect it to be transparent).

gdalwarp has a -dstalpha option that does the trick of making the nodata areas transparent - does gdal_calc has something similar? If it doesn't, what is the easiest way around it?


texts: How to make style in Geoserver different for last X hours?

texts:
When I want to set style for points (for example) representing events I would like to achieve that points created in last hour have different color than others. Is there any easy way to achieve it?

I know that there is a possible solution to run cron every minute and set status (example: last_hour: true / false) regarding time and than use this as an attribute to style data. But this solution is so ugly.

I have been looking to do something like this (using ECQL):

CREATED_AT >= TIMESTAMP - '1 hour'

Is it possible?


texts: Why are my vector layers with different CRS not aligning?

texts:
I'm using QGIS 3.4.1 on Windows. I am trying to use a point vector source from a public authority with OpenStreetMap data, but the points do not align with their expected locations.

Try the following:

Create a new project.
From the browser, go to the XYZ Tileserver category and add the
OpenStreetMap layer to the project.
Zoom in to Berlin, Germany.
In the browser, right-click on WFS and choose "New connection..."
Enter
https://fbinter.stadt-berlin.de/fb/wfs/data/senstadt/s_messpunkte in
the URL field, name the layer "emissions" (these are air quality
testing stations) and click "OK".
Add the layer from the provider to the project (Entwicklung
Luftqualität - Immissionen (Umweltatlas))

The default CRS is EPSG:4326, the OSM "CRS" is EPSG:3857 and the CRS of the point layer with the air quality testing stations is EPSG:25833. (After adding the OSM data, the project's projection changes to EPSG:3857, as expected.)

The emissions layer seems to be displaced. Consider the point (1493533,6893411) in the emissions layer (the top point in the image below). Taking the OSM image as a reference, this point should actually be at (1493539,6893156) (roughly where the bottom point is).

Why are these layers not aligned? The potential causes I have considered are:

EPSG:3857 is incompatible with the CRS (EPSG:25833 - ETRS 89 UTM Zone 33N) of the point layer. But I pulled a set of street data from the same provider (which uses the EPSG:25833) and it is reprojected so that it accurately matches the streets in the OSM data.
I have a poor combination of CRS and/or am using the wrong CRS
for the project.
The data source ("emissions") is somehow broken.
There is a bug in QGIS :P

This installation is using proj 5.2.0 and GDAL 2.3.2.

EDIT: The cause turned out to be an error in the data source.


texts: GeoServer Integrated GeoWebCache and Layer Group Mode

texts:
I've been told that in previous versions of GeoServer it was possible to cache the individual layers which make up a layer group and that requests for the layer group would then return the relevant tiles if they are seeded.

In the latest versions of GeoServer you can integrate the GWC with WMS and there are 4 options for layer group types. As far as I can tell the options work as follows:
 - single = this appears to generate a seperate tile cache for the layer group
 - named tree = same result as single
 - container tree = error, doesn't produce any tiles
 - Earth Observation tree = error, doesn't produce any tiles

Is it possible to cache the individual layers rather than the layer group in the latest versions of GeoServer or do I have to double up the file space and cache the tiles twice?

Do any of the layer group types impact the caching process or is that option unrelated to tile caching?


texts:
I have absolutely no experience with GIS, it just so happens that one of the hobby projects I am doing at the local hackerspace includes tracking a model boat with a GPS... I hope this is the right place to ask! [If not, I would appreciate suggestions on where to turn for help].

Here's my problem: today we went testing the new model boat, but when I processed the log, I noticed that for some strange reason, some of the GPS readings are totally off (see picture below: we never actually sailed on the left of the trees).

I was wandering if there is a standard way to filter a GPS log to spot and eliminate such noisy signals. The only obvious way I could think of was to tell the filter the a maximum speed for the boat that was tracked, and simply exclude any signal that would imply the boat to have travelled fastest than that maximum speed. Is there a better way?
Since our boat is a model one and most of the times it travels at speeds below the knot, I was wondering if there is a clever way to "rectify" points that are probably wrong. An example: suppose the boat is travelling east to west and the track looks something like:

...⋅...

it is quite possible that the central dot is a glitch in the GPS, but rather than dropping the log signal altogether, I would like to "correct" it [each log record contains information about the electronics too, and this data is always correct]. Again: is there a standard way to "smooth" a GPS track? I thought some sort of solution could pass through the standard deviation, but I am clueless on how to proceed.

Additional information: we are using python for processing the log, but while code snippets are always appreciated, I'd be more than pleased with a general explanation of the algorithm that I would implement myself.

Many thanks in advance for your time and expertise.


texts: What the difference: mapnik and mapnik2?

texts:
I use Ubuntu 12.04 LTC, postgres 9.1, postgis 2.0.

I install mapnik.

$ sudo add-apt-repository ppa:mapnik/nightly-trunk
$ sudo apt-get update
$ sudo apt-get install libmapnik mapnik-utils python-mapnik

Which version of mapnik I install? https://github.com/mapnik/mapnik/wiki/Mapnik2
in this documentation write that mapnik may be 0.7 or 2.1

>>> import mapnik
>>> mapnik.mapnik_version()
300000

In synaptic: 3.0.0+dev20130906.git.f3b8645-1~precise1


texts: Enabling/disabling NMEA sentences on u-Blox gps receiver?

texts:
I want to disable the $GPGSV sentences on u-Blox NEO-7M because I have serial buffer overflow on Arduino.
I found in the Receiver Discription p.69 and tried sending this:

$PUBX,40,GSV,0,0,0,0,0,0*59

But the receiver keeps sending the $GPGSV.


texts:
I have added vector layer(features) on the map.
when i click on the feature , in the back end it has to take the latlon values and retrieve the data from database(using php) and display the values in the popup.
i am able to retrieve the values from the database, but the values are displaying in another page but not in the popup my code is

function createPopup(feature) {

    var temp=temponclick(feature.geometry.x ,feature.geometry.y,name);
            popup1 = new GeoExt.Popup({
        title: 'My Popup',
        location: feature,
        width:200,

        html: temp,
        maximizable: true,
        collapsible: true,
        anchorPosition: "auto"
    });

    popup1.on({
        close: function() {
            if(OpenLayers.Util.indexOf(vectorLayer.selectedFeatures,
                                       this.feature) > -1) {
                selectCtrl.unselect(this.feature);
            }
        }
    });
    popup1.show();
}

vectorLayer.events.on({
    featureselected: function(e) {

        createPopup(e.feature);

    }
});

  mapPanel.map.addControl(selectCtrl);
  selectCtrl.activate();

 function temponclick(x, y,name) {

     var url = 'temperature.php?';
 var query = 'x=' + x + '&y=' + y+ '&name=' + name;

 window.location = url + query;
}

temperature.php

while( $row = mysql_fetch_assoc($result) ) {
echo '<input type="label" id="temp"  value="' . $row['temperature'] . '"  />';

}

?>


texts: Stacking overlapping rasters, keeping cell values on top?

texts: I have a land cover classification that consists of three files. One file is the LCC for the western half of my study area, another covers the eastern half (these two overlap in the middle), and the third represents a singe vegetation type which is present in both of the other two images (overlaps both). I would like to stack them western, eastern, single veg type, from bottom to top respectively, and end up with a raster that keeps only the top most cell values. Basically, I want a single image with the same appearance as the three when stacked. I'm working with ArcMap 10.2 and Erdas Imagine 2013 - so solutions in either software would be helpful.