texts: How to estimate the width/height size of output image in a target CRS without gdalwarp?

texts:
I have a hard problem to find the width/height of a 2D GeoTiff file
which is reprojected from a source CRS (e.g: EPSG:4326) to target CRS
(e.g: EPSG:3857) without running gdalwarp. Because, the file could be
large (GBs) and I don't need the projected output file, just need the
width/height of it.

What could be done with gdalwarp is:

gdalwarp  -s_srs EPSG:4326 -t_srs EPSG:3857 full.tif full_3857.tif

gdalinfo full_3857.tif

Size is 879, 811 (width/height).

I've searched a lot and what seems to be good without doing the real
projection is this gdal function: GDALSuggestedWarpOutput2
https://github.com/OSGeo/gdal/blob/106c8288e7a05f4efc1a588c5a3b2da7ec52d915/gdal/alg/gdaltransformer.cpp#L354.
However, it doesn't help because my application developed in Java and it
uses gdal-java (GDAL JNI) https://search.maven.org/#artifactdetails%7Corg.gdal%7Cgdal%7C1.11.1%7Cpom
as library.

Unfortunately, I cannot find this GDALSuggestedWarpOutput2()
from gdal-java http://gdal.org/java/overview-summary.html then cannot
invoke this C++ function from my application to test.

Can anyone please give me a hint to solve this problem?


texts: Calculating distances between two geometry columns using R?

texts:
I'd like to calculate the pairwise distance between two geometry columns.  That is, I want to calculate the distance between several pairs of geometries without calculating the full distance matrix between every possible pair.  I've come up with a way to do it with purrr::map2_dbl, but it feels messy.  Is there a better approach?

library(sf)
library(purrr)
library(dplyr)

# Load sample data and project to EPSG:5070.
nc <- st_transform(crs = "+init=epsg:5070",
  st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE))

# Here's the full distance matrix. I don't want that.
nc_dist_mat <- st_distance(nc, nc)

# Define a function to calculate distances pairwise:
pairwise_distance <- purrr::partial(purrr::map2_dbl, .f = sf::st_distance)

# distance to self is always zero:
pairwise_distance(st_geometry(nc), st_geometry(nc))
# distance to random county isn't always zero:
pairwise_distance(st_geometry(nc), st_geometry(dplyr::sample_frac(nc)))


texts: How to create a table which showing attributes from one layer and another?

texts:
I have two layers.

The A12_zone system which is basically the zone system of UK

and then another layer with some developments.

Each layer has its own attribute table.

I would like to create a table where each development will have a corresponding number of zone depending on in which zone it falls.
So for instance the development zone mm falls in the 151 zone and that should be indicated in the final table!

Can you help me doing this in QGIS?


texts: Split line into the different lengths from attribute table data

texts:
Is it possible in QGIS to split one line in one direction into multiple segments of varying length which I have in the attribute table?

For example: overall line is 10 m long, in table I have data: 1.5m, 2.8m, 3m,.. In this order will the line be split to specific follow-up lengths in table.


texts: Lat/Long co-ordinates are not plotting onto EPSG:3857 OSM basemap correctly

texts:
Excuse me if this is a stupid question, but I've honestly tried to find the answer in a few threads here on SE, and I cannot seem to follow the suggested steps correctly.

Project is in EPSG:3857 (WGS 84 / Pseudo-Mercator)
BaseMap (OpenStreetMap) layer is set to EPSG:3857
Vector data has been imported from CSV, and has Lat/Long coordinates, set to EPSG:3857.

The layer appears right in the middle of the map, around 0,0 obviously.

Why is this, what am I doing wrong?


texts: Removing No Data/Black Space from Raster (*.tif) using ArcGIS Desktop?

texts:
I am running ArcGIS Desktop 10.4.1 with a basic license and I am working on a community plan for a first nations community. We contracted a company to do some aerial photography of the area. The images they provided to us are in multiple formats, all displaying the same information (*.IMD, *.TIF, *.TIL, *.JPG). I have been working with the *.TIF versions of the files because they have the greatest detail. Our study area required them to provide us with 3 separate images to span the entirety of the area.

The issue that I am having is the 'black space' or what I assume to be no data due to the clipping that this company performed on the data is creating visual slivers when displaying due to the overlapping nature of the three images.

To confirm that the black space is in fact 'no data' I need an attrubute table. From what I have read online, it appears that I need to flatten the raster (currently 3 bands) to a single band to be used in the build raster attribute table tool. From ESRI support (Article ID: 000009462) it states to use the copy raster tool followed by the build raster attribute table. I was able to do this, the tool ran successfully but was unsuccessful in flattening the raster to a single band, it remains in three bands (RGB), not allowing me to build an attribute table.

Will flattening and building an attribute table for the raster allow for me to remove these 'no data' values so I can display these images without slivers as they are to be used in meetings and final reports?

If so, what would be the steps/tools needed to make this happen?

If not, I would like a definite answer if this is at all possible, I have scoured the interwebs for answers and many people appear to have similar issues with no definite answers.

It is important to note that I have reviewed steps to make this happen using the greyscale function in the image analysis window (which worked in flattening the raster). I need to keep colours in tact for these images, greyscale is not an option.


texts: Generation of pixel statistics in multiband raster

texts:
I want do generate for each pixel the average and standard deviation or better the variance coefficient for raster with 57 bands. The average should be average of the value for all bands. The result should be one raster, text file or something similar.

I tried a lot of tools. so far I have not been able to solve my problem.

I have no experience with Python, which is maybe my main problem.
Have you some ideas? If you recommend me solutions with Python could you do it foolproof?

The similar post from 2014 couldn't help me.


texts: How to access and send coordinates from third party application to the current running instance of ArcMap?

texts:
I'm quiet new to ArcObjects and I need some help. My main issue is that I have an external WPF application (a video player exactly), and from that application I want to  perform the following actions in the current instance of ArcMap:

From my application get the currently running instance of ArcGIS desktop or ArcMap
Send coordinates and display a bitmap at that location

Frankly I don't know from where I should start, so here are my questions:

Is it possible to access the current running session of ArcGIS from a third party program?
Is it possible to send coordinates to ArcGIS this way?
Should I convert my application to an AddIn? How can i perform that, if you have any useful indications it will be great.


texts: ESRI Javascript API Geocoder VS Locator module

texts:
What's the difference between ESRI's Geocoder and Locator module in the javascript API? I know the Geocoder defaults to ESRI's geocoding service, or both modules can be set to a custom Address Locator service on ArcGIS Server, which is my use scenario. Is the Geocoder intended to be used more like a "widget" because it adds a search box? Is there a reason you'd use one module over the other, or a case where'd you'd want to use both?
Thanks for any info or opinions offered.


texts: Failing to have visible features in ArcMap 10.5

texts: How do I increase visibility of my map? Feature labels are too tiny and features are too blurred. I am using ArcMap 10.5. I tried to change the resolution to 3000 dpi but still didn't work out. I also tried to increase the layout view.


texts: How to join tables to create a new shapefile in ArcMap?

texts:
I'm unsure how to join a table and reference shapefile in order to create a new shapefile. I'm using ArcMap. I made a simplified version example below.

I have a shapefile of points (locations where an indvidual can reside). The points have a unique ID. The points don't have a specified lat/long. I also have a table of individuals. Each individual has a has residence information, in the same form as the point ID. More than one individual can be found at the same point.

I want to create a point shapefile that shows every individual's residence. I know that some of the points will overlap.

I thought I would have to use the "Add Join" tool, but when I used it, it didn't create a shapefile. Is there some special join tool that will create a new shapefile/table that will hold the point data and the individual data?

I don't know much about geocoding, but would it be useful in this situation? Any other thoughts would be a huge help?


texts: Adding new field in ArcMap from data of existing field?

texts:
I want to create a new field with years data.

There is an existing field with the timestamp YYYYMMDDHHMM but I just want the years from that field to the new field.

That way I can display just the years when I create a graph.

Is there anyone that can explain the code to do this?

I am not a programmer.


texts: Multi-sheet printing of large-scale maps

texts:
Mapping a specific region at a larger scale, this region does not fit onto a single sheet of paper. For example mapping Austria at 1: 1 mio needs three sheets of din A4 paper.
Question: Can I do multi-sheet printing of large-scale maps in QGIS?


texts: Applying ImageCollection reduction to just one band in Google Earth Engine

texts:
I have an image collection that is a binary layer that indicates whether there is a burn or not in each pixel, for the specified time period by which I filtered the satellite images. This is calculated from the differences in the Normalized Burn ratio (dNBR), calculated for each image. A pixel has value 1 (which means it is burned) if the dNBR value is above the threshold value for a pixel, both when looking at the difference between NBR of day1 and day2, and looking at the difference between NBR of day1 and day3. Now I want to reduce the image collection to an image such that each pixel that burnt at any point of time (within the given time period) has a positive value and attach a timestamp as a band to the image, such that each burnt value pixel has a time of burn. The way  I am implementing it is by adding a date band to the original filtered image collection, then finding the dNBR and then reducing the image collection by max() reducer. But the problem is the reducer is maximising even the time band and I am getting only one value as the time stamp. Is there a way to apply reduce to only one value and retain the value of the other band which the maximum value was part of? I have been breaking my head over it for quite some time. Please help

var vis = {"opacity":1,"bands":["B4","B3","B2"],"gamma":3.464},
sent = ee.ImageCollection("COPERNICUS/S2")

var land = sent2.filterBounds(geometry).filterDate("2019-01-01","2019-06-01")
.filterMetadata("MGRS_TILE","equals",'44RLT')
.filter(ee.Filter.lt("CLOUDY_PIXEL_PERCENTAGE",5));
var landsort = land.sort("system:time_start");
print(landsort);

var addDate = function(image){
  var doy = image.date().getRelative('day', 'year');
  var doyBand = ee.Image.constant(doy).uint16().rename('doy')
  doyBand = doyBand.updateMask(image.select('B8').mask())

  return image.addBands(doyBand);
};

var a = ee.ImageCollection(landsort.toList(landsort.size()));
var a = a.map(addDate)
var a = a.toList(a.size())
var b = function(image){
 var nbr3 = ee.Image(image).normalizedDifference(["B8","B12"]);
 var nbr2 = ee.Image(a.get(a.indexOf(image).subtract(1))).normalizedDifference(["B8","B12"]);
 var nbr1 = ee.Image(a.get(a.indexOf(image).subtract(2))).normalizedDifference(["B8","B12"]);

 var e = nbr1.subtract(nbr2).gt(0.3);
 var f = nbr1.subtract(nbr3).gt(0.3);
 var g = e.multiply(f);
 var h = g.addBands(ee.Image(a.get(a.indexOf(image).subtract(2)))
.metadata('system:time_start'));
 return h;
};
var dnrcollect = ee.ImageCollection(a.map(b)).max();
Map.addLayer(dnrcollect)


texts: Configurable Map Viewer Advanced Search Widget

texts: I add Advanced Search Widget on my configurable map viewer application, I have also configure proxy file and paste Attributes Tables, the problem is that search by attributes, by adresse works but search by spatial don't. where is the problem ? I am using WFS layers.


texts: How to modify a LYR file using ArcObjects?

texts:
I'm creating an Add-in for ArcCatalog using C# that modifies the datasources of layers in a MXD and .LYR files by themselves.  I have gotten the MXD code to work but I'm lost with how to modify a .LYR file without the use of IMap.  Does anyone have any code suggestions?

To load an MXD and loop and start modifying layers in the TOC was simply something like:

UID pUID = new UIDClass();

            pDoc.Open("C:\\Test\\File1.mxd");

            for (int i = 0; i < pDoc.MapCount; i++)
            {
                pMap = pDoc.get_Map(i);

                UID uid = new UIDClass();
                uid.Value = "{40A9E885-5533-11d0-98BE-00805F7CED21}"; //IFeatureLayer
                IEnumLayer layers = pMap.get_Layers(uid, true);       //Returns all layers

But I'm unsure how to get this process going for standalone .LYR files..

ILayerFile l = new LayerFileClass();
l.Open(fileName);


texts: Opening InMemory Table to execute some SQL gives Runtime Error 445 - Object doesn't support this action?

texts:
I have a InMemory that gets created, what I then want to do is run INSERT, UPDATE and DELETE queries on the table. The code i have at the moment is below, I create a new InMemory table, add it to the map (which works).

    Dim pInMemTable As ITable
    'create Temp DP table
    Set pWSF = New InMemoryWorkspaceFactory
    Dim PName As esriSystem.IName
    Set PName = pWSF.Create("fred", "InMemory", Nothing, 0)
    ' Debug.Print ">>" + PName.NameString + "<<"
    Set pFWS = PName.Open

    Dim pFlds As IFieldsEdit
    Set pFlds = New esriGeoDatabase.Fields
    pFlds.AddField MakeField("OBJECTID", esriFieldTypeOID)
    pFlds.AddField MakeField("DP_REF", esriFieldTypeString, 10)
    pFlds.AddField MakeField("WANTED", esriFieldTypeString, 1)

    Dim pUID As New UID
    pUID.Value = "esriGeodatabase.Object"

    Dim ptable As ITable
    Set ptable = pFWS.CreateTable("DP_Temp", pFlds, pUID, Nothing, "")
   ' Debug.Print ptable.

    Dim pSAT As IStandaloneTable
    Set pSAT = New StandaloneTable
    Set pSAT.Table = ptable
    pSAT.Name = "DP_Temp999"

    Dim pMxDoc As IMxDocument
    Set pMxDoc = ThisDocument

    Dim pSATC As IStandaloneTableCollection
    Set pSATC = pMxDoc.FocusMap

    pSATC.AddStandaloneTable pSAT
    pMxDoc.UpdateContents

    Set pInMemTable = pSAT

    Dim path As String
    path = "InMemory/DP_Temp999"
    If pWSF Is Nothing Then
        MsgBox "pwsf is nothing"
    Else
       ' MsgBox "pwsf is something"
        'MsgBox pWSF.WorkspaceDescription(False)

        'pwsf.Open(
        Set pWS = pWSF.OpenFromFile(path, 0) ' **********PWS doesnt get set***********
        If pWS Is Nothing Then
            MsgBox "PWS IS NOTHING"
        Else
            MsgBox ("Yes!!")
        End If

    End If

but the code fails where PWS is getting set as it always returns as Nothing. I have tried various ways of setting the path value e.g. "InMemory/DP_Temp999", "In_Memory/DP_Temp" and trying the "\" as opposed to "/" but it still comes back as nothing.
Does anyone have any ideas please? or is there a better way of executing SQL on a table that will allow SELECT, DELETE and INSERTS.

****UPDATE**
Right I have been working on this for a few days have have been able to set pWS, by casting by ptable to a idataset object and then setting pWS to the pdataset propertyset properties which works. However I have then tried the following code, to see if SQL can be executed

    If pWS Is Nothing Then
        MsgBox "PWS IS NOTHING"
    Else
        'MsgBox ("Yes!!")
        Dim pWorkspaceProperties As IWorkspaceProperties
        Set pWorkspaceProperties = pWS

        Dim pWSProperty As IWorkspaceProperty
        Set pWSProperty = pWorkspaceProperties.Property(esriWorkspacePropertyGroup, esriWorkspacePropCanExecuteSQL)
        If pWSProperty.IsSupported = True Then
            MsgBox "SQL is supported"
            Dim sql As String
            'sql = "select * from 'DP_TEMP';"
            sql = "UPDATE dp_temp set dp_temp.[DP_REF] = 'Test';"
            pWS.ExecuteSQL sql
        Else
            MsgBox "SQL is not supported"
        End If
    End If

When the code is run, the pWSProperty.IsSupported returns as True, however when the .ExecuteSQL method is run, I get the following error Runtime Error 445 - Object doesn't support this action. So I am a bit confused as to whether SQL can be run or not.


texts: Combining Data from a CSV file to a Shapefile based on matching FID

texts:
I am relatively new to GIS.

I have an existing shapefile of New Zealand islands that details island name, Lat/Long, island group, area, etc. Included in this shapefile is a unique FID for each island. This information is quite extensive (polygons are provided down to the rock stack level: <0.01 ha in size) but is missing some information. Islands that are quite small, for example, often do not have a name. For my analysis I wanted to add in some of this classification information.

I have imported a CSV file from excel that attempts to fill in some of this missing information: I have manually added in names and archipelago groups that were missing from many islands. The FID of these islands remains the same as the shapefile, I just used excel to "fill in the holes" because excel is so easy to manipulate.

I want to append the information that I have added in the CSV file to the existing NZ island shapefile. From my (limited) experience using a JOIN in GIS adds in a new column. Instead, I want to fill in the existing holes of the shapefile with information from my CSV file, based on the unique FID.

Does this make sense?


texts: How to center the map to a KML placemark with OpenLayers 3?

texts:
I'm using OpenLayers 3 to display my KML layer on the map.
My KML contains two placemarks, one contains a linestring of waypoints and the second one is a point, like this:

<kml xmlns="http://www.opengis.net/kml/2.2"
xmlns:gx="http://www.google.com/kml/ext/2.2"
xmlns:xal="urn:oasis:names:tc:ciq:xsdschema:xAL:2.0"
xmlns:atom="http://www.w3.org/2005/Atom">
<Folder>
 <Placemark>
   <name>Asset 5</name>
    <description>
      Path
    </description>
   <Style>
   <LineStyle>
    <color>41FA7800</color>
    <width>2.5</width>
   </LineStyle>
  </Style>
 <LineString>
   <extrude>0</extrude>
   <tessellate>1</tessellate>
   <altitudeMode>relativeToGround</altitudeMode>
  <coordinates>
    ... lots of coordinates ...
  </coordinates>
 </LineString>
 </Placemark>
<Placemark>
  <name>Current location</name>
  <description>Current location </description>
  <Point> <coordinates>50.1384717,11.2500542</coordinates> </Point>
</Placemark>
</Folder>
</kml>

I would like to center the map at the position of the second Point placemark.

How can I do this?


texts: ArcPy Export SQL Server to Shapefile gives RuntimeError: Object: Error in executing tool?

texts:
I have executed a stored procedure and want to export this layer to a shapefile.  I feel like this should be an easy solution.  I tried using arcpy.CopyFeatures_management and it threw an error.  I am not sure if CopyFeature is meant for sql server tables.

I am running ArcGIS 10.1.  I am using SQL server 2012 with geography data type.  The table is not registered with the geodatabase.  The error is below:

RuntimeError: Object: Error in executing tool
<<< Process finished. (Exit code 1)

Below is my code:

import pyodbc
import arcpy

sqlserver = pyodbc.connect('DRIVER={SQL Server};SERVER=sqlserver;DATABASE=AndrewSandbox;Integrated Security=True')
cursor = sqlserver.cursor()
table = cursor.execute("SELECT * FROM TEMPLINE")

arcpy.CopyFeatures_management(table, "D:\Python\WaterSage\Viability\TestExport.shp")


texts: Georeferencing screen shot of Bing imagery in ArcMap?

texts: I've been using the world imagery layer in ArcMap for hand-digitizing possible restoration sites. This has generally been adequate except for a few instances in which the imagery is a bit old and significant changes (such as the addition of an irrigation pond) have occured. In general the imagery available from Bing is adequate but I don't have a key to access it. I was thinking about just taking screen shots of the specific sites where I've notices discrepancies and georeferencing- does anyone have experience with this, or can you recommend another source for up-to-date imagery of the US pacific northwest?


texts: Using "Join Attributes by Location" to perform spatial join of layer with joined data on dissolved layer

texts:
I'm relatively new to QGIS and have been using it for a project related to gerrymandering and electoral redistricting. My apologies if the question may seem a bit novice.

I have two layers that I've imported into QGIS: One with the boundaries of all zip codes in Virginia and one with all electoral precincts in Virginia. The precinct file also included a column for cluster which notes that electoral district a precinct is in. I used GDAL Dissolve polygons to dissolve the precincts up into the associated cluster/district as shown below:

After importing the zip code layer, I used Joins to attach data to the zip code attribute table, matching based on the zip code identified (labelled ZCTA5CE00). The joined table looks like this:

One of the columns I attached to the open attribute table for the zipcode layer is idealpnum2.

I want to do a spatial join using "Join Attributes by Location" to map the zipcodes onto the cluster/district layer and get an average idealpnum2 for each cluster/district. The layers do not match perfectly and district/cluster boundaries cut through zip codes in several places. Is this possible?

When I've tried to do this using the dissolved cluster/district layer as the target vector layer and the zipcode layer as the join vector layer, it returned a joined layer with an attribute table that was null for everything. . In the pictures below, I had first went into the zip code layer and changed the CRS to EPSG:4326 so that it had the same CRS as the dissolved cluster/district layer.

I also tried to run it without changing the CRSs and instead "forcing" it through. Did not work there either.

Does anyone know what I should be doing differently here? Is it even possible?


texts:
I'm using QGIS-2.4.0 to developing an application.

I created my own custom MapCanvasItem for showing SVG. It is great to load a SVG from symbol library that QGIS provides(such as d:\QGIS Chugiak\apps\qgis\svg\symbol\poi_cave.svg), or a SVG file that QGIS Composer exports.

But when I created a SVG using Adobe Illustrator or Inkscape on my own, the application broke down to loading a SVG. I tested the SVG file created by Illustrator/Inkscape, it is well to display in QT project when I create a class Inheritanced from QgraphicsItem.

I don't know why, Can anyone help me?
Here is my code to rewrite the function from QgsMapCanvasItem:

void TFMapCanvasSvgItem::paint(QPainter* painter)
{
  if (!painter)
  {
    return;
  }
  QImage svgImage(50, 40, QImage::Format_ARGB32);
  QPainter imagePainter(&svgImage);
  mSvgRenderer.render(&imagePainter);
  painter->drawImage(QRect(0, 0, 100, 80), svgImage, QRect(0, 0, 100, 80));
}

EDIT more:

The SVG symbol created by Illustrator/Inkscape couldn't be used in QGIS API application, but could be used in QGIS release software(such as Symbol selector, SVG annotation), as showed in 1.PNG, 2.PNG.

1.PNG

2.PNG


texts: Generating random points stratified by category within sampling hexagons using ArcGIS for Desktop?

texts:
I am trying to use two different layers together to generate a layer of random points that represent sample sites. The two layers are:

A grid layer of hexagon-shaped polygons created using the [Create Hexagon Tessellation][1] tool.
A polygon layer with four categories of percent impervious surface (%IS), created by reclassifying the NLCD %IS raster layer and converting to a vector polygon layer. Let's say the categories are 0-25%, 25-50%, 50-75%, and 75-100%.

I would like to generate 25 random points for each category of %IS (100 points total), but only one point per hexagon. So 25 hexagons that have a point placed within a %IS polygon representing 0-25%, 25 hexagons that have a point placed within a %IS polygon representing 25-50%, etc.

I have tried using the built-in Create Random Points tool, the [GME Generate Stratified Random Points tool][1], and the Sampling Design Tool. However, I have only been able to create either

a points layer with one random point inside each hexagon but not associated with the %IS information (Create Random Points), or
a points layer where the points are stratified correctly by %IS but not constrained by the hexagon layer (GME, Sampling Design Tool).

Is there a way to do what I'm attempting?

Here's an example of what a small portion looks like. The colors of the %IS polygon layer represent the four different classes (increasing %IS from blue to red).


texts: When we type an address in search bar, in some web mapping applications like google maps, a list of words automatically displayed from the database.

texts:
What is actually happening behind the code and how should be the code?
when we search for a specifc word, will it look in the database (here it is postGIS) and lists the similar words.


texts: How to create Concentric Circles in QGIS?

texts: I'm trying to build a planning tool for center pivot irrigation system.  It involves selecting a field (polygon), a center point, and then an outer circle to set the outer limits of the system.  From there you select the inner rings based on fixed length sections. to optimize the coverage area.  Anyone have an idea where I can start?


texts:
I am using a NTv2-grid for transformation of geodata.
The grid can be found here:

gsb-file:
https://www.lgl-bw.de/lgl-internet/web/sites/default/de/05_Geoinformation/Galerien/Dokumente/BWTA2017.zip

gsa-file:
https://www.lgl-bw.de/lgl-internet/web/sites/default/de/05_Geoinformation/Galerien/Dokumente/BWTA2017_gsa.zip

I set up a directory watch with an automatic transformation workflow so workmates can just drop a shapefile in my input-folder and get the transformed geodata saved in the output-folder.

In order to prevent transformation of data which lies outside the NTv2-grid-extent I would like to get the extent of the NTv-2-grid so I can do a Within-Check before transforming the data.

Does anyone know a way to extract the points of the grid to calculate a convex hull? The given extent in the header is way too imprecise and using the polygon of the federal state for which the grid was set up is also not accurate enough for a within check.


texts: Distance from within boundary in Shapely

texts:
I'm new, although I am reasonably experienced in Python.

My question concerns Shapely, a library for spatial data. I am trying to assign some points as being within some distance d of the boundary of the shape. The main problem is that these shapes are in some cases very, very nasty indeed, with self-overlap and intersection galore - some examples illustrate what I am dealing with:

My initial approach was to try constructing buffered objects that were +/-d from the original shape using the buffer() method. But this doesn't appear to work for me - the shapes are either overlapping, completely all over the place, or way too far in to be sensible.

I figured that another approach was to simply get the distance() to the shape. But this does not work for those points within the shape - distance() merely returns 0.

I am wondering if there is any way of getting the distance of a point to the boundary of a shape from within that shape, or more generally a way to do what I am trying to do. Any suggestions?

An example Python snippet of the approach that seems to work best is as follows:

# A function that reads in a likely invalid shape and tries to return
# a 'valid' one.
def construct_shape(poly, distance, delta, update):
    difference = 0
    if type(poly.buffer(distance)) is MultiPolygon:
        difference = update(difference, delta)
        while type(poly.buffer(distance + difference)) is MultiPolygon:
            difference = update(difference, delta)
        new_poly = poly.buffer(distance + difference)
        return new_poly.buffer(-difference)
    else:
        return poly.buffer(distance)

for (file_i, polygon_file) in enumerate(sorted(glob.glob("./line/*.annotations"))):
    # Do some processing with the files, which contain the points that define the shapes
    # tags is a list of points extracted from the current file
    # ...

    poly = Polygon(tags)

    # Iterative increasing approach

    # dilated is +500 from the shape
    dilated = construct_shape(poly, 500, 1, lambda x, y: x + y)
    # contracted is -500 from the shape, or in this case -1000 from dilated.
    contracted = construct_shape(dilated, -1000, 1, lambda x, y: x - y)

    if type(contracted) is MultiPolygon:
        print("Uh oh. Contracted for " + polygon_file + " invalid.")
        int_x = [i[0] for po in contracted for i in po.exterior.coords]
        int_y = [i[1] for po in contracted for i in po.exterior.coords]
    else:
        int_x = [i[0] for i in contracted.exterior.coords]
        int_y = [i[1] for i in contracted.exterior.coords]

    if type(dilated) is MultiPolygon:
        print("Uh oh. Dilated for " + polygon_file + " invalid.")
        ext_x = [i[0] for po in dilated for i in po.exterior.coords]
        ext_y = [i[1] for po in dilated for i in po.exterior.coords]
    else:
        ext_x = [i[0] for i in dilated.exterior.coords]
        ext_y = [i[1] for i in dilated.exterior.coords]

    # Plot the results
    plt.plot(ext_x, ext_y, color="r")
    plt.plot(int_x, int_y, color="b")

    patch = PolygonPatch(poly, fc="red", alpha=0.5, zorder=2)
    ax.add_patch(patch)

    plt.show()

The above is my current approach. What I had in mind for it would be something like this:

inner_margin = contracted
outer_margin = dilated

point = Point(...)

if point.within(outer_margin) and not point.within(inner_margin):
    # do some stuff with point
else:
    # ignore

For the other approach regarding points being some distance from the margin:

point = Point(...)

# Assume we just got boundary from Polygon(tags)

if point.distance(boundary) <= 500:
    # do some stuff with point
else:
    # ignore


texts: How do I combine bands in a tiff then convert to grayscale?

texts:
I have a large mosaic of tiffs in which I need to do the following:
-Combine bands 4,3,2, for an infrared image then convert the image to grayscale without producing another dataset
-This is easy on QGIS but I need to execute this on ArcGIS 10.1

Is there a way to render the image like on QGIS?


texts: How to use Python console under QGIS 2 with gdal module?

texts:
I am trying to move my plugin from QGIS 1.8 to QGIS2.0. However, the gdal module cannot use under QGIS2.0.
I try to use Python console, typing:

    import gdal  #it's OK
    gdal.GetDriverByName( 'GTiff' ) # dosen't work, got an AttributeError

But the 'smart spelling' give 'GetDriverByName' attribute.(and has 'two'?! one for method and another one for attribute?? )
I also typing:

    dir(gdal)

But No these methods or attributes found

What can I do??  Please help me.


texts: r.watershed analysis error from raster map being in different zone and unable to be opened

texts:
I tried running my DEM (geotiff) file through the r.watershed program by the GRASS GIS plugin of QGIS.  I believe the error message results from a raster map in the different zone than the current region but I set my GRASS GIS mapset in WGS84-Zone 16N coordinate system where it ought to be.  Thank you in advance.

Here is the error message I get:

r.watershed elevation=dem2@Tejeras threshold=30 accumulation=dem2_accum drainage=dem2_drain stream=dem2_stream basin=dem2_basin

SECTION 1a (of 5): Initiating Memory.

￼Raster map is in different zone (16) than current region (0)

￼Unable to open raster map

￼unable to open elevation map layer

￼Subprocess failed with exit code 1

￼category information for [dem2_accum] in [Tejeras] missing or invalid

￼category information for [dem2_drain] in [Tejeras] missing or invalid

￼category information for [dem2_basin] in [Tejeras] missing or invalid

￼category information for [dem2_stream] in [Tejeras] missing or invalid
Finished with error


texts: Mustache functions possible in CartoDB tooltip templates?

texts:
CartoDB.js uses the mustache.js library to render tooltip windows. One of the features of mustache.js is allows functions to be defined, which are called when the template is rendered. An example is given below:

{
  "beatles": [
    { "firstName": "John", "lastName": "Lennon" },
    { "firstName": "Paul", "lastName": "McCartney" },
    { "firstName": "George", "lastName": "Harrison" },
    { "firstName": "Ringo", "lastName": "Starr" }
  ],
  "name": function () {
    return this.firstName + " " + this.lastName;
  }
}

Is it possible to pass such a function to a CartoDB tooltip template? I understand how basic variable substitution works, I am talking about the case where a function is called to generate the substitution. In the above example that is:

 "name": function () {
    return this.firstName + " " + this.lastName;
  }

which would be referenced as {{name}} in the template.


texts:
The current behavior: I migrated from @types/openlayers": "4.6.17" to "@types/ol": "5.3.6" and I don't find the function getVectorContext(event) in ol/render so I added it manually to render.d.ts but the next npm install I will lose this modification

The wanted behavior: I want to add a file *.d.ts that contain a custom types to my project but I got the error :module not found can t resolve '/path/customtypes'


texts: Split millions of points in shp with 100'000 points each

texts:
I havie over 8 million points in a shapefile. Now I would like to split these 8 mio points into parts of 100'000 points each, in order to facilitate further calculations with these points.

In the end I need a table containing the r,g,b,nir values of every pixel of an image, at the location of these points AND within a specific radius around these points (buffer).

I intended to use the "sample" tool to get this table.
But unfortunately it seems that sample can't handle millions of points (gave me an empty table after 20 hours of calculating when I tried to sample 2 million of my 8 million points (worked perfectly fine for 10'000 points)).
This is why I'm looking for a method to automatically iterate my 8 mio points and split them into smaller parts, over which I then can iterate in further steps to do the buffering, rasterizing and sampling.

The smaller shapefiles should still contain the attributes: x,y,Radii and NewID.

As I'm new to Python, I have no idea how this could be done, but I seek a Python solution if possible.


texts: RS & GIS (QGIS Plugin) not producing Land Surface Temperature Layer

texts:
Recently I downloaded this plugin as it has demonstrated capability to produce Land Surface Temperature Calculated layer. I tested this plugin using the Landsat 8 contents straight from USGS Earth Explorer. I opted for using compressed zip file as the base data for the plugin to run the calculation. I ran this plugin to calculate FCC (False Colour Composite), NDVI and LST. The plugin was able to produce FCC and NDVI without a single issue (I have tested few times). However, I am encountering issue where the LST layer is not able to be formed.  Here is the error I am facing:
For Land Surface Temperature:
Error:
Traceback (most recent call last):

  File "C:/Users/Acer/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\NITK_RS-GIS_17\RSGIS_M.py", line 2147, in run

    lst_cal_w_del(num)

  File "C:/Users/Acer/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\NITK_RS-GIS_17\RSGIS_M.py", line 1850, in lst_cal_w_del

    temp_cal(num)

  File "C:/Users/Acer/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\NITK_RS-GIS_17\RSGIS_M.py", line 1804, in temp_cal

    (metadata_required[num][3][0][1] / output_ra[num][t_info[1][1]]) + 1)) - 273.15

TypeError: unsupported operand type(s) for /: 'float' and 'list'

As I am not very familiar with Python, may I ask kind assistance for me to explore a workaround with the dataset or in the plugin to rectify the issue?


texts: Access / modify QGIS copyright decoration through python

texts:
Is there a possibility, i.e. API, for PyQGIS to access / modify the copyright decoration (View -> Decorations > Copyright Label) in QGIS?

I've looked in the Python QGIS API (http://geoapis.sourcepole.com/qgispyapi/qgsproject), however, the read/write methods are only for plugin settings, as far as I can tell from the examples I found, e.g. at https://docs.qgis.org/2.14/de/docs/pyqgis_developer_cookbook/settings.html.

Unfortunately, at this point I don't know what to search for and hope somebody can point me in the right direction. Hints are also welcome at: https://github.com/geometalab/Vector-Tiles-Reader-QGIS-Plugin/issues/158


texts: Return value from Django view to JS function

texts:
I am building a webmap using GeoDjango, Leaflet and PostGIS. I am trying to write a function that returns a raster value from a database when the user clicks on the map (i.e., via a JS alert, or other means). I have read a number of posts on SE but can't seem to get the code to work myself.

The relevant code is below... first, the JS functions:

map.on('click', function(e) {
    var lng = e.latlng.lng;
    var lat = e.latlng.lat;
    console.log("I clicked!");
    resp = updateLatLng(lat, lng);
    console.log('The response is: ', resp);
});

function updateLatLng(lat, lng) {
    var data = {'lat': lat, 'lng': lng, csrfmiddlewaretoken: '{{ csrf_token }}'};
    $.post('{% url 'legend' %}', data);
    var response = '{{result}}';
    return response
};

The legend view in views.py is as follows:

def legend_view(request):

  lat = request.POST['lat']
  lon = request.POST['lng']

  query = "SELECT ST_Value(raster, ST_TRANSFORM(ST_SetSRID(ST_MakePoint(%s, %s), 4326), 4326)) FROM my_raster_layer;" % (lon, lat)

  conn = psycopg2.connect("dbname='mydb' user='postgres'")
  cur = conn.cursor()
  cur.execute(query)
  result = cur.fetchall()[0][0]

  print("\n\n++++++++++++\nInside the legend view\n")
  print("Lat: ", lat, "   Lon: ", lon)
  print("Value: ", result)
  print("++++++++++++\n\n")

  return HttpResponse(result)

When I click on the map, it successfully sends the lat and long values to views.py, as well as successfully queries the database. I can see this because it prints in the terminal I'm running the webmap from. However, it does not seem to be returning the value to the JS, as result is printing nothing.

Some other relevant SE items that I've followed, but can't get to work for my particular case:

https://stackoverflow.com/questions/1445989/passing-python-data-to-javascript-via-django
https://stackoverflow.com/questions/298772/django-template-variables-and-javascript?noredirect=1&lq=1


texts: Setting Numeric scale bar in map composer QGIS

texts:
I'd like to display the scale bar at a specific cartographic scale.  I'm able to figure out how to display it as a Numeric value but would like to assign it a different value than the one that it's set at.

For example it loads as 1:16260 but I'd like to set it to 1:10000.


texts: How to filter time-related images in Google Earth Engine and conduct regression analysis to verify the consistency of different sensors？

texts:
Landsat data harmony is shown on the Google developer website(https://developers.google.com/earth-engine/tutorials/community/landsat-etm-to-oli-harmonization), but I want to verify this harmony result by selecting sample points to merge. What should I do?
Below is a link to the original code(https://code.earthengine.google.com/798ae7a268f8e2c8022433b9562785c0)


texts: How to increase maximum memory usage cap for running Arc 10.1?

texts:
Not sure if this is strictly GIS related or perhaps belongs on StackOverflow...it is regarding ESRI ArcGIS 10.1, so I'm sticking it here.

I'm looking for options to increase processing capacity, and I've read unhelpful bits and pieces about increasing the cap on memory usage for ArcGIS.

I'm working on ESRI ArcGIS 10.1 (all licenses) on a 32-bit Windows 7, 3.33GHz Intel DuoCore, 4GB RAM computer. I'll eventually be upgrading to a more powerful, 64-bit computer with 8GB Ram, and would like to know what to do in the meantime, or when the new one gets in.

Thanks!


texts: QGIS Rule-Based Labeling - Toggle child rules visibilty

texts:
So it works to add labels by dragging it to another label, like it does in styles. But if I try to hide this upper filter, the contained labels won't hide.

Is this a bug, or is it not wanted, to set upper filters? Is there a workaround?

I've got to set 10 filters for 4 different labels, so I will get 40 filters after that. I've got no problem with that, but there will be many people, that will have to work with that styles/labels, and I don't want to confuse them.

My QGIS Version ist 2.18.


texts: Intersecting multiple feature classes with just one feature class?

texts:
I have multiple feature classes (fc) and would like to intersect them with one other fc, is this possible without merging/joining the multiple fc's?
Using ArcGIS Pro I've attempted to intersect numerous point layers and a line layer but I actually want to intersect the numerous points to just the line. I've just joined the points to do this but I wanted to avoid this, if there is a way around this?


texts: Where are the button icons stored in QGIS?

texts: I'm trying to create a custom icon theme to make my install of QGIS a bit nicer, but I can't find where the icons are stored. All I've found is the QGIS logo in /usr/share/qgis/images/icons.


texts: Clipping from a tiled LAS dataset with LASTools (for ArcGIS)

texts:
I have huge amounts of lidar data stored as tiled datasets per project. The "tiles" are separate LAS files based on our national map sheet setup (800 by 600 meter). To keep an overview of the data I have set up LAS dataset definitions per project in ArcGIS 10.1.

I've been testing out the LASTools for ArcGIS toolbox to see if I could use this for clipping, as that's a function I'm really missing in ArcGIS. What I often want to do is to define a clipping area, and copy data within the area to a new LAS file. This area will normally contain parts of multiple LAS files:

However, it seems the lasclip tool only will clip from one LAS file at a time. What I'm looking for is a way to say - "give me the data from this LAS dataset in this area, and save it to a new LAS file". Am I missing something, or is that not possible here?

Edit: I've also tried posting this to the Lastools GoogleGroups, but without any response.


texts: R get_map returns map correctly, but filename argument does not work?

texts:
In R library ggmap there is a function get_map which saves a file (or at least have saved files automatically two or three years ago). Now get_map is not saving anymore, although the parameter filename gives no warning, no error.

m <- get_map(c(lon0,lat0,lon9,lat9),maptype="satellite",filename='map1.png',messag=T)

Is get_map still expected to create a PNG file?


texts: Error running Arc Hydro tools function in ArcMap?

texts:
I am trying to run a function from Arc Hydro tools in ArcMap 10.1; Sink evaluation. But before it can finish it returns an error message, that, frankly, I don't understand at all.
It goes System.runtime.interopservices.comexception (0x80004005): error hresult e_fail .

Can anyone help me solve this? I have tried all I can think of. Run the function from the toolbox instead of dropdown menus, shorten paths, restart ArcMap, reinstall AHT, open new map documents, restart my computer, wait till next day and see if it automagically have fixed itself...

After searching this forum I found a possible duplicate, although it's got no answers. Now I'm all out of ideas, can anyone please help me?!

I'm using ArcMap 10.1 SP1 basic license and Windows 7.

I think I have narrowed it down a bit, still not clear what's wrong though. Sink evaluation uses a series of functions in Arcmap, of which one is Zonal statistics as table. That's the one it gets stuck on. The input raster somehow does not exist, path unknown. I guess this is some intermediary raster.

After consulting the Arc Hydro team, I got the advice not to use grids stored in a gdb (which I didn't). I have also tried to change target locations, thanks to @G.Oxsen. Unfortunately neither of those solved the issue. According to the Arc Hydro team, GPL0 that I asked about earlier is an in-memory grid, so probably no problem there.

I am currently pursuing that it is related to the installation of SP1. Anyone out there using AHT with SP1 and got Sink evaluation working?

Uninstalling SP1 did not change anything, unfortunately. Still working on this, but with less intensity. Seems like I'll have to find workarounds...

By reducing the number of sinks in need of processing I managed to get Sink evaluation running. Having the files not stored in a gdb seem to help as well, as mentioned earlier, and using integer grids.
Changing the "Maximum number of unique values to render" in ArcMap Options (Customize menu), tab Raster --> Raster dataset is another tip that I got (from the Arc Hydro Team). Restarting ArcMap after running every one or two functions also helps keeping things stable.


texts: Algorithm to optimally split a map comprising many separate areas into panels for printing

texts:
I have a map of set of land parcels (polygons) that need to inserted into a document. Some of these polygons are contiguous, but some can be scattered randomly. The area is typically too large for a one page map. So the map must be divided into several tiles/panels. Each tile/panel is of fixed scale and size, but the corner coordinate doesn't have to be on a fixed grid.

A naive approach is to use a regular grid and printing any tile that touches one of the polygons. However this leads to really sub-optimal results. For example if a polygon contains a grid corner you have 4 tiles when one tile could do the job.

The following image illustrates the problem: 4 tiles/panels should be sufficient to print this map, but the naive approach doubles that to 8. This will be fun algorithm figure out, but I'm sure this problem has been solved before. Any pointers (or even keywords to use in a search) would be helpful.

I'm writing custom non-interactive/batch oriented software (in Java), so algorithms, papers or libraries is what I'm most interested in. However pointers to software products (particularly if open source) might be helpful also.


texts: I want create a custom script but the Toolset menu appear disable.


texts: Adding video to ArcGIS API for JavaScript popup template?

texts:
I was wondering how do I add a video to be played in the popup template?

I've looked at the documentation here: https://developers.arcgis.com/javascript/jssamples/widget_extendInfowindow.html

and can't seem to find any information on it?


texts: PostgreSQL: How to use more than one CPU?

texts:
I'm creating a real time traffic application. I read the AVL car data from a web service and insert the records in a PostgreSQL table.

Then a trigger calculates the closest link and if there is a previous position for the same car, I calculate the route to get the distance travelled in that time frame, then I get the average speed for that route.

The thing is, I get around 3000 points each 5 minutes, so I need to process ~600 records per minute (10 each second). I'm very close to those numbers because I filter the road network before calling the routing function and the routes are very small.

I understand PostgreSQL doesn't support multiple threads in one connection but supports multiple connections, so could use the other CPU cores.

I can make a different table for each State (25), but how can I partition the problem so each trigger creates a new connection and usse more CPU cores instead of only one?


texts: Problem with the synchronization between QField and QGIS

texts:
I have a serious issue with synchronizing QField and QGIS. I'm working with .qgs project that I export to QField via QFieldSync, then I added a new feature and every thing works perfectly even when I synchronize with QGIS I can see my new feature, but when I want to edit the feature that I've recently created (change the value of a field) the synchronization with QGIS doesn't seem to be working, it preserves the old values.

If you have any idea about what is happening tell me please.


texts:
I've recently noticed that the UTM zone displayed in the Identify box and the "Goto XY" box display a different UTM zone than the one I'm in.  I'm working in WGS84 NAD83, Zone 13N but when I use the Identify tool and look at the location of the point it shows the location as '13S X,Y'.  Also, if I type in a set of coordinates manually using the "Go To XY" tool, I type it in as 13N X Y, and hit enter, after it flashes the location the zone I had input is changed to 13S.

Can someone explain why the zone is being displayed as 13S and not 13N?

I'm assuming this is one of those weird geometric things that I never picked up in school.  Hopefully someone can enlighten me!


texts: I'm trying to install TimeManager for QGIS. I can't find it under Plugins/Fetch Python Plugins and I have verified that the official repository is listed in the Plugins menu, so I'm wondering if there's some other way that I'm supposed to install it. Underdark, am I missing something obvious here?


texts: How to speed up simultaneous working with PostGIS layers in QGIS?

texts:
we are working on one big layer at the same time (10-20 people). Some of us use QGIS 2.8.1 and the other 2.6.
The user with 2.8 have a lot of crashes (QGIS). The 2.6 users don't have this crashes but to add a feature is incredible slow.
The Postgres 9.4 server is running on a Virtual Machine (Hyper V) with Win7.
Maybe it's the network and not QGIS? Maybe someone have an idea. Would be great :)


texts:
What I'm trying to do is to merge to geotif files with gdal_merge, but the result is not what I expected.

Here are the files I'm trying to merge.

gdal_merge.py -n 0 -a_nodata 0 -o result.tif file1.tif file2.tif

The two tif files

The result


texts: How to connect QtDesigner with python console on QGIS

texts:
I've created a script with python console on QGIS and an interface for my plugin with Qt Designer with QGIS 2.14.0 custom widgets.

I would like to connect my sript with the interface but I can't find the way to do it, should I use python console in QGIS or jus the command in my computer?

And what should I write to connect these two things?


texts: GeoServer: Legend title with single feature layer

texts:
Been working with GeoServer for some time now. While working with the publishing tab of the layer, I noticed a problem.  Every time that I set the default style with a single featured CSS style, the defined title doesn't show. For some clarification I've added some screenshots:

The first image shows the desired visualization. Notice the two titles in the legend
The second image shows the desired style but not the title, this always happens if a layer is only styled by one feature.

I would like to figure out a way to make my titles visual in all my layers, can anyone help me out?


texts:
Learning the js api i tried to make a querytask that will query a layer and update graphics on change.

I think the problem is in my success("succ") function.

I made this codepen:

https://codepen.io/segev-salman/pen/joBwyW

in the 3.28 version of the api I saw a code example building an array and a for loop inserting features into a graphic from the querytask reault, but the 4.11 doc says that querytask returns a graphic so does that mean i can throw the result straight into a new graphic layers ?


texts: How to stop highlighting using QgsHighlight?

texts:
I am highlighting certain features of a layer using QgsHighlight.
On a certain event, I would like to stop the highlighting and highlight something else or nothing at all. I have tried QgsHighlight.dropEvent(), QgsHighlight.updateCanvas() or just to refresh the whole canvas but that won't help and I don't see anything else in the documentation.


texts: Creating graphical marker/point on Geotools map in Java?

texts:
I'm writing a program that tracks the International Space Station, and plots its current position on a geotools map.

I currently have a geotools map that pops up, and I have the current location of the ISS that updates every second in another class:

public void doIt() throws IOException, URISyntaxException {

    URL url = this.getClass().getResource("/coastlines_z8.shp");
    System.out.println(url);
    File file;
    file = new File(url.toURI());

    FileDataStore store = FileDataStoreFinder.getDataStore(file);
    SimpleFeatureSource featureSource = store.getFeatureSource();

    MapContent map = new MapContent();
    map.setTitle("map");

    Style style = SLD.createSimpleStyle(featureSource.getSchema());
    Layer layer = new FeatureLayer(featureSource, style);
    map.addLayer(layer);

    createPoint();
    JMapFrame.showMap(map);
}

This seems like it should be simple, but I've been scouring the internet for the past day and I can't find a definitive answer as to how to go about doing this. I've looked into SimpleFeatureTypeBuilder, but I'm unsure how to implement it.

Additionally, once I get the point drawn, how do I ensure that it is updated every second?


texts:
I am using ArcGIS 10.2.2 for Desktop.

I have a model made in ModelBuilder which currently does:

input XY coordinates excel table
data display XY event
converting a format event shape
add field "name1"
add field "name2"

So far so good but my problems are:

the model only works so edition
I cannot add tool to make a selection of my new shapefile


texts: Qgis -> Settings -> Customization : how to import ini file with python?

texts:
With QGIs, from "Settings"->"Customization" i can hide/show the components of the GUI.
I can export a INI file with my customization.

There is a way to import my INI file with a python command in the console, without open the customization dialog?


texts:
I load dbf files in R and then I manipulate it because I only want some features and after that  I export it.

But when I will open the shapefile with the same name  in a GIS software features are not located in the correct place. Can it be the result of the change of the primary key of the dbf file?

Code:

recfe<-read.dbf("RECFE.dbf", as.is = FALSE)
recin<-recfe[which(recfe$MUNICIPIO==86),]
write.dbf(recin, "RECFE.dbf", factor2char = TRUE, max_nchar = 254)

Any suggestion?


texts: I need to be able to sample many raster files (imported from GeoTIFFs) using a previously defined set of spatial points, using R.  What's the best way to do this?


texts: How to get feature on vector layer in Openlayers 3?

texts:
I'm try to get feature of a vector layer, but I'll have 3 or more Vector layer. When I click the map, I didn't get the feature. My vector layer definition;

var vectorAltAlan = new ol.layer.Vector({
        source: new ol.source.GeoJSON({
            projection: 'EPSG:2322',
            url: 'data/geojson/guneykaya_alt_alan.geojson'
        }),
        style: function (feature, resolution) {
            ...
    });

So how can I detect the layer and getting featue by click ?


texts: Is it possible to export dynamic labels from ArcMap to KMZ? (Or alternate approach)

texts:
I am preparing increasingly more maps and layers for export to KML/KMZ as it seems to be the most convenient way to quickly analyze spatial data by my superiors.
In order to export labels I seem to have to convert them to annotations. Even this approach is not without limitations; the annotations seem to be rasterized and almost always illegible unless I make them very large and adjust the pixels and DPI values of the output image.

Is there a way to export labels directly, or as vector? If not is there another simple way to export labels to KML/KMZ?

COMMENT:
Thank you all for your suggestions. It seems that the existing workflow seems to be the best option. I exported my annotations to polygons and this workaround seems to work pretty well.

Arc2Earth produces erroneous KMZ files that are riddled with errors.

The recommended ArcGIS script is only available for legacy ArcGIS 9.x versions.

Dynamic Labels DO NOT export at all into KMZ using the default tool.


texts: How can i turn off the minus to zoom control in OpenLayers while using OpenLayers.Control.KeyboardDefaults?

texts:
I have a map using OpenLayers.Control.KeyboardDefaults but I also have a form next to the map that acepts coordinates. As is stands pressing '-' for a negative coordinate does not type but zooms the map out.

Is there a parameter I can add between the brackets of new OpenLayers.Control.KeyboardDefaults() that will allow me to turn off the listener for '-' while retaining all the other controls?

Thanks


texts: Extracting Geometry Coordinate from a feature

texts:
I am using Google Earth Engine, I have a FeatureCollection in which I extract a single Linestring into a "Feature" item.
(JSON output is included below)
I can extract any items from the "properties" objects using the ee.Feature.get() function.
However, I am unable to grab the coordinates from the "geometry" section.
How can I grab the coordinates of this line and export them to separate variables?

-edit-
I was able to find a workaround by using GIS, making table entries for begin and end X and Y values, the read the data in Earth Engine with a .getInfo() call from the new properties entries.


texts: Loss of edge pixels when displaying clipped rasters

texts:
I am looking for tips to clip properly and cleanly 1 raster in 4 rasters without loss of pixels at the edge.

I am using QGIS Valmeria. I have a raster named 'A' representing the Percent of trees cover (PTC), from 1000m elevation to 4000m elevation. From it, I want to get 4 rasters representing the PTC at 1000m, 2000m, 3000m and 4000m elevation.

So using polygons of the four elevations, i clipped the raster 'A' (command raster>extraction>clip). for the rasters obtained, i write '0' for the additional null value (in order to remove the frame around).

But when i display together the four rasters, there is a kind of stream of 1 or 2 pixels large between the rasters displayed. It looks like the pixels at the edge of the clipped rasters are missing. The pixels lose represent PTC data from 0 to 100%.


texts: SAGA GIS, how to generate groundwater contour lines from points?

texts:
I'm trying to generate contour lines/isolines from groundwater points (groundwater depth) with SAGA (I tried before with QGIS, contour plugin, but the result was not good looking, too angular lines).

I have uploaded my groundwater points into SAGA, but I can't find the right tool or the steps to draw the contour lines. Can somebody tell me the steps to do it?

Do I need to generate a grid first ?


texts: Not all KML attributes import into attribute table

texts:
I have a .KML file that has placemarks that have four main pieces of data that I want to use in QGIS.  When I do 'Add Vector Layer' and select my KML file it imports all the placemarks as feature points as expected. But when I right click on the layer and select 'Open attribute table only two of the attributes are there.

I need these four attributes in QGIS.

<name>9255 LECALDERWOOD LN</name>
<description>9255</description>
<address>9255 LECALDERWOOD LN, San Fransi, CA</address>
<ExtendedData>
<Data name="Terminal Address"><value>r 9252 meridian av</value></Data>
</ExtendedData>

Example Placemark:

<Placemark>
<name>9255 LECALDERWOOD LN</name>
<description>9255</description>
<address>9255 LECALDERWOOD LN, San Fransi, CA</address>
<ExtendedData>
<Data name="Terminal Address"><value>r 5252 meridian av</value></Data>
</ExtendedData>
<Point>
<coordinates>-121.893896,37.249197,0</coordinates>
</Point>
<Style>
<IconStyle>
<Icon>
<href>https://www.qcdata.s3.amazonaws.com/img/earthit/house2.png</href>
</Icon>
</IconStyle>
</Style>
</Placemark>

Attribute table:

Does anyone have any ideas how to get it to import all attributes?


texts: Dimension Difference preventing plot of series

texts:
I have performed linear regression between the 'ro','pr' series (runoff,precipitation) and printed the chart. I want to print the linear fit on top of it and get the R**2 value, but dimensionality prevents it. yValues is list of 1 element with size 98 and xValues.map(lin) (x values mapped in regression function) is a list of 98 elements. I tried converting both to ee.List, ee.Array but does not work. Any solution?

var sro = ee.ImageCollection('IDAHO_EPSCOR/TERRACLIMATE') // selection based on filters
.filterDate("2017-12-30","2019-12-30")
.select("ro","aet",'pr')
.filterBounds(geometry02);

var linearFit2 = imgRegress_pre.reduceRegion({
  reducer: ee.Reducer.linearFit(),
  geometry:small_region,
  scale: 30,
});

//------------------------------------------Regression---------------------------------------
print("---------------Pre-fire statistics-----------------------------------")
print('OLS estimates:', linearFit2);
print('y-intercept:', linearFit2.get('offset'));
print('Slope:', linearFit2.get('scale'));

var s1=sro_pre.mean();

// reduce region
var avPixelVals_pre = s1.reduceRegion(
    {reducer: ee.Reducer.toList(), geometry: small_region, scale: 1500});

var yValues = avPixelVals_pre.toArray(['ro']);

// Get the red band value list; to be plotted along the x-axis.
var xValues = ee.List(avPixelVals_pre.get('pr'));

var lin = function(x){
  var y =ee.Number(x).multiply(linearFit3.get('scale')).add(linearFit3.get('offset'));
  return y;
};

var Yarr=ee.Array.cat([yValues,xValues.map(lin)]);

// Define the chart and print it to the console.
var chart8 = ui.Chart.array.values({array: Yarr, axis: 1, xLabels: xValues})
                .setSeriesNames(['Runoff','Reg'])
                .setOptions({
                  title: 'Relationship between Runoff and Precipitation',
                  colors: ['red','blue'],
                  pointSize: 4,
                  dataOpacity: 0.7,
                  hAxis: {
                    'title': 'Precipitation (mm)',
                    titleTextStyle: {italic: false, bold: true}
                  },
                  vAxis: {
                    'title': 'Runoff(mm)',
                    titleTextStyle: {italic: false, bold: true}
                  }
                });
print(chart8);


texts: QGIS 1.7 + MacOS 10.6 : georeferencing from map canvas bugs

texts:
I am experiencing issues with the gdal plugin to georeference a raster. It works fine when I enter the coordinates manually but it does not work when I try to link the plugin with the map canvas.

Does anyone had this issue and was able to solve it?

Thanks

Guillaume

QGIS 1.7 MacOS 10.6.7


texts: How to edit a tiff aerial image in ArcMap 10.3

texts:
(beginner with arcmap)

I have aerial images in tiff format and I need to edit them in order to draw polygons on them.

When trying to edit them I get a message that says that there are no layers editable.

I guess I need to convert the tiff layers I have in order to be able to edit them; but how do I convert them? to what do I convert them?


texts: Mapserver not displaying a raster properly

texts:
I have a 32bit floating point raster containing only five values (-2, -1, 0, 1, 2 + -999000000 for noData don't ask me why, it was delivered to me), no compression, projected (EPSG:3857). I am using mapserver in combination with OpenLayers to display it in a web application as a tile layer. Everything seems ok and no errors are reported however this is what I get:

It looks like some tiles are not loaded (straight border of yellow colour in the upper middle) moreover there shouldn't be any transparent pixels over the land, water is noData. As I zoom in an d out as well as pan pixels and tiles seem to randomly appear and disappear.

Here is my basic mapfile I also tried different combinations of PROCESSING configurations but none of them work:

MAP
  NAME "qqqq"
  CONFIG "PROJ_LIB" "../../share/proj/"
  STATUS ON
  UNITS METERS
  WEB
    IMAGEPATH "/ms4w/tmp/ms_tmp/"
    IMAGEURL "/ms_tmp/"
    METADATA
      "wms_title"     "WMS qqqq"
      "wms_onlineresource" "http://xxx.xxx.xxx.xx/cgi-bin/mapserv.exe?"   ##required
      "wms_srs" "EPSG:4326 EPSG:3857"
      "wms_enable_request" "*"
      "WMS_FEATURE_INFO_MIME_TYPE" "application/json"
    END
  END

  PROJECTION
    "init=epsg:3857"
  END

  #
  # Start of layer definitions
  #

  LAYER
    NAME "products"
    PROCESSING "NODATA=-999000000"
    TEMPLATE "templates/NDVITemplateJSON.js"
    TOLERANCE 5
    TOLERANCEUNITS PIXELS
    METADATA
        "wms_title"    "WMS qqqq"
    END
    VALIDATION
        "provider" "[a-z]+"
    END
    TYPE RASTER
    STATUS ON
    DATA '/xxx/xxx/xxx/%provider%/%spaceresolution%/%timeresolution%/%year%/%month%/%day%/%name%.tif'
    PROJECTION
      "init=epsg:3857"
    END
    PROCESSING "RESAMPLE=NEAREST"
    CLASSITEM "[pixel]"
    CLASS
      EXPRESSION ([pixel] = -2)
      STYLE
        COLOR 227 26 28
      END
    END
    CLASS
      EXPRESSION ([pixel] = -1)
      STYLE
        COLOR 231 245 54
      END
    END
    CLASS
      EXPRESSION ([pixel] = 0)
      STYLE
        COLOR 221 221 221
      END
    END
    CLASS
      EXPRESSION ([pixel] = 1)
      STYLE
        COLOR 0 245 0
      END
    END
    CLASS
      EXPRESSION ([pixel] = 2)
      STYLE
        COLOR 26 150 65
      END
    END
  END
END

I also tried to define colours based on ranges:

EXPRESSION ([pixel] >= 1 AND [pixel] < 2)
…

It didn't work.

Any ideas what might be wrong?

EDIT:

I found out the tif was produced by GrADS software

gdalinfo T2m-60_perc_2018_231.tif

…
Metadata:
  TIFFTAG_SOFTWARE=GrADS version 2.0.a7.1
  TIFFTAG_XRESOLUTION=0.018242123
  TIFFTAG_YRESOLUTION=0.02681388
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  AREA_OR_POINT=Area
…

And this might be the cause as I also found this question, unfortunately left unanswered.


texts: Cannot get CMV / Disclaimer Widget up and running

texts:
I am a Noob to javascript so forgive my ignorance.

I have moved all folders and associated files from the download to the CMV Dijit folder with exception of the config file which I placed in the config folder.

Do I need to modify any existing script to get this thing running?
Does it run "out of the box" like the actual CMV app?

I am running Windows OS.


texts: Field calculator add sequential number to text field to populate a unique ID

texts:
I am using arcmap 10.1 and working in a personal geodatabase. I am trying to add a sequential number to an [encounterID] field to populate a [sightingID] field. What code would i use?

For example I want to sequentially number my encounter ID's and populate the sighting id field as it appear in the sighting ID below and for multiple records.

Column 1

[encounterID]
20170327.0301
20170327.0301
20170327.0301

Column 2

[sightingID]
20170327.0301_1
20170327.0301_2
20170327.0301_3


texts: How to show attributes of a new feature in Attribute table area when using IEditor3 and IEditTemplate, ArcGIS 10

texts:
ArcObjects SDK for .NET and ArcGIS 10 Desktop

Quick snippet to set the stage:

IEditor3 editor.StartOperation();
IFeature newFeature;
// fill out newFeature.Shape, etc.
newFeature.Store();
editor.StopOperation("Adding new feature");

At this point I'd like to automatically "show" the attributes for this new feature in the Attribute display/table, which could be visible alongside the Create Features template list.

Equivalent process through ArcMap itself:

Start the Editor via the Editor Toolbar.
Click a feature template from the Create Features list.
Click the Attributes button to open the Attributes display.
Mouse-down in ArcMap.  See new/default attributes for the new feature be displayed in the Attributes table area.

Is the possible with ArcObjects for .NET?


texts: What's the proper ArcPy usage of 'TemplateExtent' in CreateFishnet_3d?

texts:
While running the standalone tool in ArcMap, the optional parameter TemplateExtent is a wonderful shortcut, as it autofills 10 other parameters. However, I need to run it in a standalone script, and have been unsuccessful in using the parameter as it's used in the tool: leaving these fields to be autofilled using the TemplateExtent layer's info. Instead, as seen in Creating fishnet from template feature class using ArcPy? (courtesy of radouxju), the best solution seems to involve using Describe to find and enter each parameter explicitly:

desc = arcpy.Describe(fc)
arcpy.CreateFishnet_management(fc[:-4]+"_c200.shp",str(desc.extent.lowerLeft),str(desc.extent.XMin) + " " + str(desc.extent.YMax + 10),"200","200","0","0",str(desc.extent.upperRight),"NO_LABELS","#","POLYGON")

That approach squares with the fact that the first of these parameters are required. If I'm understanding it correctly, when run in ArcPy it has no utility since all other fields are required.

Question: So should I always ignore this paramter in the ArcPy environment, just using '#' as a placeholder; or is there some other utility?

Following PolyGeo's and Midavalo's suggestions below, here's the Python snippet that results from running CreateFishnet as a tool:

arcpy.CreateFishnet_management(out_feature_class="[path]/test_fishnet", origin_coord="447677.5 3111557.5", y_axis_coord="447677.5 3111567.5", cell_width="50", cell_height="50", number_rows="0", number_columns="0", corner_coord="482222.5 3142082.5", labels="LABELS", template="mid_5m", geometry_type="POLYLINE")


texts: QGIS: Comparing attribute values from different layers with Expression string builder

texts:
I want to compare attribute values from different attributes from multiple layers in the Expression string builder in QGIS (3.16.1) in of these layers, but I don't know how to do that.
I've tried out aggregate to get values from other layers, but I can only count, summarize etc. values. But I want to compare the values.
For example there is the attribute "OBJID"(from the first layer). And I want to say OBJID (from the first layer) = OBJID (from another layer). So that all objects in the first layer are displayed where this equation is true.


texts: Color ramp in multiple attributes symbology (Arcmap)

texts:
I am trying to use the "multiple attributes" symbology tool in Arcmap 10.

I would like to divide my data in 2 categories. For each categories, I would like to display a color ramp.

I read ESRI's documentation about that, but I can't figure how to obtain the proper color ramp:

Screenshot

I want to have a "green" color ramp for the first category (value '0'), and a "red" color ramp for the second one (value '1').

What am I doing wrong?

Thanks for your help


texts:
I have a vector layer containing watercourse location which I have rasterized. I also have a DEM containing elevation data for the region in which the watercourses are located. I would like to reclassify my rasterized watercourse layer such that pixels which contain a watercourse are equal to 1, and pixels which do not contain watercourses are equal to 0. I would then like to combine this new watercourse layer with my existing DEM to generate a layer containing watercourses showing the elevation at each pixel along their length.

I have attached a picture describing what I would like to achieve.

Would anyone please be able to provide a step-by-step guide for reclassifying my rasterized watercourse layer and for generating my desired end product? I am using QGIS 3.4 Madeira.


texts: Debugging network analyst solve error 030024

texts:
I'm trying to build service areas for a low-voltage utility network, but am having trouble solving the network analysis. Here is the code to build the service area:

Treatment_SA = arcpy.na.MakeServiceAreaLayer (network, district+"_Treatment_SA", "Length", "TRAVEL_FROM", "200" , "DETAILED_POLYS", "NO_MERGE" , "DISKS", "NO_LINES" , "NON_OVERLAP", "NO_SPLIT", Feeders, "", "", "", "TRIM_POLYS", distance_from_line)
Treatment_SA = Treatment_SA.getOutput(0)

#Load Elements
arcpy.na.AddLocations (Treatment_SA, "Facilities", Injection_TX, "", "25 Meters")
arcpy.na.AddLocations (Treatment_SA, "Polygon Barriers", Existing_TX_Polygons, "", "25 Meters")

#Solve
arcpy.na.Solve(Treatment_SA,"SKIP")

I'm not too worried about errors in the code since it has run successfully with other data, but I'm including in this post in case it is helpful. When I run this code, I receive error 030024:

  Solve returned a failure. No edges were traversed. Failed to execute (Solve).

How can I find out more about what is causing this error?


texts: Pathfinder office (ssf) to ArcGIS (shp) projection issue in ArcMap

texts:
I'm honestly not sure if that's even the problem but here's what's going on:

Collected points with Trimble Geo7X set to NAD83 UTM zone 15N. Post-processed the points in Pathfinder Office, transformed to a shapefile and opened in ArcMap. There was some issue there where I had to open in ArcCatalog and set the projection again.

The points seemed to display correctly in ArcMap at first, but when I added a TIF as background imagery it was clear that they were in the wrong spot. I used Calculate Geometry to see the coordinates and they were correct but in decimal degrees.
Ex. (-93.652034, 42.031155). If I use Identify, the coordinates are different and in dms (97 29'22.498", 0 0'1.365") for the same point. Seems like they displayed in the wrong projection and coordinate system but I'm not sure how.
Adding the correction summaries:

Differential Correction Summary:
    1 file processed.  In this file:
        1164 (100.00%) of 1164 selected positions were code corrected by post-processing
        1163 (99.91%) of 1164 selected positions were carrier corrected by post-processing

    Estimated accuracies for 1164 corrected positions are as follows:
            Range   Percentage
            ----------  ----------
              0-5cm 99.66%
             5-15cm  0.26%
            15-30cm      -
            30-50cm      -
             0.5-1m      -
               1-2m  0.09%
               2-5m      -
                >5m      -

Differential correction complete.

Differential Correction Summary:
    1 file processed.  In this file:
        3585 (99.97%) of 3586 selected positions were code corrected by post-processing
        1422 (39.65%) of 3586 selected positions were carrier corrected by post-processing
        1252 (88.05%) of code positions chosen over carrier, as they were of higher quality

    Estimated accuracies for 3585 corrected positions are as follows:
            Range   Percentage
            ----------  ----------
              0-5cm      -
             5-15cm  0.59%
            15-30cm  0.31%
            30-50cm  0.11%
             0.5-1m 20.22%
               1-2m 63.26%
               2-5m 14.28%
                >5m  1.23%

Differential correction complete.


texts: QGIS 3.X CAD tools not working with snapping on PostgreSQL layers?

texts:
If I create a scratch layer in EPSG 32720 CAD Tools are working fine with snapping options. The vertex are highlighted etc. etc.

If I use a PostgreSQL layer, same EPSG, the snaping is working, but not highlighted and angles are not recommended etc. etc.

Does anyone else have the same problem?


texts: How to select features with same values?

texts:
How can I select all those features that has exact same values in their attribute field called "X"?
I think that the solution could be something like the Countif function in Excel, but I haven't found such function in QGIS. I know that there is the Group Stats plugin, but I would like to select values and later work with them.


texts: Where should the data update date be entered into the FGDC metadata form (esri)?

texts: if you've edited data recently (a gdb feature class), either spatial or tabular, and you want to capture the "update date" and possibly a note about what was updated, without overwriting the "publication date" or "edition" feilds in the metadata, where should it go?  Perhaps i'm interpreting the meaning of the publication/edition dates incorrectly and that's what they're for...? any known/definitive answer to this?


texts: What is the relationship between ITRF and WGS84?

texts:
Why are there people saying that "I am using WGS84(ITRF96)"?

There are newer version of ITRF. Will WGS84 follow it?

If so, what will be the influence to us?


texts: OTB Pansharpening Error: Adapter for adaptPansharpening-bayes not found

texts: When I try to use the OTB pansharpening tool inside the processing tool box I always get this message: "Adapter for adaptPansharpening-bayes not found". I guess it's an issue withe the configuration. Has anyone a solution for my problem?


texts: LiDAR GUI tool for selecting and saving groups of points?

texts:
I'm working in a small group that is processing LiDAR data with some custom code. We need to see how well our stuff is working so we need samples of "truth" data that we can compare our output against.

Ideally we could load a .las file into a Graphical User Interface (GUI) then select a series or group of points that can be saved off. Then, we could programmatically compare the output to the truth.

We will have thousands of point groups to manually save so I'm trying to find the most efficient way to do so.


texts: Interpolating OS terrain data - national grid reference NU

texts:
When using the interpolation plugin over the OS terrain data for national grid square NU, the pixels of the resultant raster are not square, but rectangular. Furthermore, I'm not able to produce viewsheds from the interpolated terrain data. I've conducted multiple tests against other national grid squares which work perfectly well.

I've shared the contours for you to have a go at interpolating. Settings are shown below.

If you are able to produce terrain data from my shared shapefile with square pixels, please could you share the result and mention any ideas as to what I might be doing wrong?


texts: Getting number of new features to be stored in layers using PyQGIS?

texts: In PyQGIS, how can I get the number of features that have been created in a given layer and that are pending to be saved in the layer when "Toggle editing" or "Saving layer edits"?


texts: Quickly setting value of selected features in QGIS?

texts:
I am manually selecting vector features (points from a Spatialite db) in a QGIS project, and then setting values manually for all the selected features using the Field Calculator button. This takes a lot of clicks, entering a value in the expression box and selecting the correct field (out of dozens) from the 'update existing field' option.

2016-09-26 edit: The point layer is > 100,000 rows, so using the attribute table takes a considerable amount of time.

Is there a quicker way to do this, perhaps a one-line Python command from the Python console?


texts: python calling def __init__(self, uri=None): method

texts:
I am working on python plugins for QGIS. I have my QueryBuilder.py and postgis.py files.
I connected it to database using Connection class and under def __init__(self, uri=None)  in postgis.py file.
I want to connect to database through QueryBuilder.py file.I imported postgis.py into QueryBuilder.py But how to call def __init__(self, uri=None) of postgis.py from  QueryBuilder.py? I need to create instance of 'Connection' class in order to call
def __init__(self, uri=None) method of postgis.py in QueryBuilder.py How do I go about it?

postgis.py file:

class Connection(DbConn.Connection):

@classmethod
def getTypeName(self):
    return 'postgis'

@classmethod
def getTypeNameString(self):
    return 'PostgreSQL'

@classmethod
def getProviderName(self):
    return 'postgres'

@classmethod
def getSettingsKey(self):
    return 'PostgreSQL'

@classmethod
def icon(self):
    return QIcon(":/icons/postgis_elephant.png")

@classmethod

def __init__(self, uri=None):

    DbConn.Connection.__init__(self, uri)

    #self.host = uri.host()
    #self.port = uri.port()
    #self.dbname = uri.database()
    #self.user = uri.username()
    #self.passwd = uri.password()

    self.host ='localhost'
    self.port =5432
    self.dbname ='ssss'
    self.user = 'j'
    self.passwd = 'jkj'

    if self.dbname == '' or self.dbname is None:
        self.dbname = self.user

    try:
        self.con = psycopg2.connect(self.con_info())
    except psycopg2.OperationalError, e:
        raise DbError(e)

    self.has_spatial = self.check_spatial()

    self.check_geometry_columns_table()

    # a counter to ensure that the cursor will be unique
    self.last_cursor_id = 0
               '''
               '''
if __name__ == "__main__":
window = Connection()

and QueryBuilder.py is:

import sys
import pdb

from PyQt4 import QtGui, QtCore
import psycopg2
import pprint

#from qgis.core import *
#from qgis.gui import *

import postgis

from DatabaseModel import TableItem, SchemaItem, DatabaseModel
import DbConnection

from ui.DlgQueryBuilder_ui import Ui_Dialog
from DlgDbError import DlgDbError

from QueryParser import QueryParser
from QueryManager import *

class DlgQueryBuilder(QDialog, Ui_Dialog):
    so on.....


texts: Google Earth Engine: How to obtain coordinates of centroids of polygons of a FeatureCollection

texts:
I need to obtain coordinates (longitude, latitude) of the centroids of a FeatureCollection of polygons in Google Earth Engine. I would like to download the coordinates if possible.
I am new to GEE and to JavaScript (normally I use R).
Here is what I have so far:
//Mapping GAUL
var gaul1 = ee.FeatureCollection("FAO/GAUL/2015/level1");

Map.setCenter(12.876, 42.682, 5);

Map.addLayer(gaul1, {}, 'GAUL');

//Subset Italy
var italy = gaul1.filter("ADM0_NAME == 'Italy'");

print('italy', italy);

// Compute centroids of each polygon
var getCentroid = function(feature) {
  // Keep this list of properties.
  var keepProperties = ['ADM0_CODE', 'ADM0_NAME', 'ADM1_CODE', 'ADM1_NAME', 'ADM2_CODE', 'ADM2_NAME',];
  // Get the centroid of the feature's geometry.
  var centroid = feature.geometry().centroid();
  // Return a new Feature, copying properties from the old Feature.
  return ee.Feature(centroid).copyProperties(feature, keepProperties);
};

// Map the centroid getting function over the features.
var centroids = italy.map(getCentroid);

// Display the results.
Map.addLayer(centroids, {color: 'FF0000'}, 'centroids');

// Print centroids
print('centroids', centroids);

How does one abstract the coordinates from the centroids?


texts: Use if statements in Calculate Field Code Block to give features in a field different values

texts:
I've been self learning Python. Now I want to play a little with the Python code block so I can be more comfortable building my own expressions later on. Right now I am just playing with this water resource layer. I created the new field called "play", basically what I want to do is write a simple expression so that if the "WETLAND_TY" field holds the value of "Freshwater Forested/Shrub Wetland". the new "play" field will return a value of 1, otherwise it will return value 5.

Below is the screenshot of my current window

I don't know what goes wrong here.

I know this issue can be dealt simply with select by attribute. But I am learning to build my expressions for some more complicated projects, so I'd like to finish it with Python expression.

In the process of playing, I made the "play' field like this:

and changed my code block to:

def NewValue(Area):
'''I was playing with a different field: Area'''
if Area >= 0 and Area < 30000:
    return 1
elif Area >= 30000 and Area < 800000:
    return 2
else:
    return 3

And this time it worked. I just don't understand what's wrong with what I wrote the first time.

Could anyone help me with this?


texts: Using the feature class name to populate a field in the features attribute table

texts:
I'm trying to duplicate the model shown in Using Feature Class "Name" from Iterate Feature Classes to Calculate Field with ModelBuilder?

I have over 700 feature classes in a gdb and I need to use the name of the feature class to populate a new field within the feature class. I tried to duplicate the model above and I can't figure it out.

What is The Feature Type Parameter?


texts: Error when trying to run GRASS 7 via processing toolbox in QGIS 2.18 on Mac

texts:
I have gone through the following steps:

Installed the following packages that came with QGIS 2.18:

1 GDAL Complete.pkg
2 NumPy.pkg
3 matplotlib.pkg
4 Install QGIS.pkg

Installed GRASS-7.4.1 from GRASS GIS for Mac by drag-and-droping it into the /Applications folder.

When I now try to run v.net.nodes in the processing toolbox, I get the following errors:

The log reads the following:

2   Uncaught error while executing algorithm
Traceback (most recent call last):
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/core/GeoAlgorithm.py", line 203, in execute
self.processAlgorithm(progress)
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 297, in processAlgorithm
Grass7Utils.executeGrass7(self.commands, progress, self.outputCommands)
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 298, in executeGrass7
command, grassenv = Grass7Utils.prepareGrass7Execution(commands)
File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 289, in prepareGrass7Execution
command = Grass7Utils.command + ' ' + os.path.join(Grass7Utils.grassMapsetFolder(), 'PERMANENT')
TypeError: unsupported operand type(s) for +: 'NoneType' and 'str'
1   There were errors executing the algorithm.

I am guessing that Grass7Utils.command is Null, which tells me, that Grass7Utils.py does not know where to search for GRASS 7.

Do I need to configure QGIS for it to know where to find GRASS 7? If yes, how is this done?


texts: Buffer Csv gets repeated when more than one feature is selected in mapcanvas in PYQGIS

texts:
I have a function which will allow the user to select a feature or more features in map canvas and performs buffer for the selected feature(s) and writes the output to a csv file. I dont have any problem if a single feature is selected in map canvas. But when many features are selected, the name of the cities in the output csv gets repeated which is shown below and my code is also attached.

def mapselection(self):
    for layer in QgsMapLayerRegistry.instance().mapLayers().values():
        if layer.name() == "us_cities":
            lsNearCities = []
            filePath = self.dlg.fileNameEdit.text()
            if len(str(filePath)) ==0:
                QMessageBox.warning(self.dlg, "File Name","Give it a file name and try again")
                return
            if self.dlg.DistancelineEdit.text():
                dist = float(self.dlg.DistancelineEdit.text())*(1609.344)
                if dist != 0:
                    selfeature = layer.selectedFeatures()
                    for feature in selfeature:
                        mbuf = feature.geometry().buffer(dist, 2)
                        iterFeat = layer.getFeatures()
                        for f in iterFeat:
                            geom2 = f.geometry()
                            valTest = QgsGeometry.within(geom2, mbuf)
                            if valTest:
                                idx = layer.fieldNameIndex('name')  #name of the attribute column with the cities in it
                                idx1 = layer.fieldNameIndex('state')
                                nearCity = f.attributes()[idx]
                                state = f.attributes()[idx1]
                                state_and_cities = str(nearCity) + "," + str(state)
                                lsNearCities.append(state_and_cities)
                elif dist == 0:
                    selfeature = layer.selectedFeatures()
                    for feature in selfeature:
                        idx = layer.fieldNameIndex('name')  #name of the attribute column with the cities in it
                        idx1 = layer.fieldNameIndex('state')
                        nearCity = feature.attributes()[idx]
                        state = feature.attributes()[idx1]
                        state_and_cities = str(nearCity) + "," + str(state)
                        lsNearCities.append(state_and_cities)

            else:
                selfeature = layer.selectedFeatures()
                for feature in selfeature:
                    idx = layer.fieldNameIndex('name')  #name of the attribute column with the cities in it
                    idx1 = layer.fieldNameIndex('state')
                    nearCity = feature.attributes()[idx]
                    state = feature.attributes()[idx1]
                    state_and_cities = str(nearCity) + "," + str(state)
                    lsNearCities.append(state_and_cities)

    #call the function to write them into a file
    self.writeToFile(lsNearCities, filePath)


texts: Retrieving Google Map data into Matlab to plot points

texts:
First of all, I hope this question isn't too basic, obvious or inappropriate for this forum. I'm very new to using map data. I'm working on a project on more accurate map routing using Matlab.

My question is:

I'm getting data from Google Map but how can I retrieve it to Matlab and plotting it?


texts: Setting mask (halo) for labels in Python script?

texts:
I currently am writing a python script and using a lot of the arcpy mapping module functionality. I have been using the LabelClass of the arcpy mapping module to automate the display of labels with my ptyhon script.

Is there a way to set a mask (halo) for labels with ArcPy?


texts: r - turn off automatic ordering of linestrings when applying sf::st_intersection

texts:
This question is about how while using sf::st_intersection on a LINESTRING and a set of POLYGONs to get several LINESTRINGs in the CORRECT order. sf::st_intersection automatically orders the newly created LINESTRINGs based on the names of the POLYGONs you have. To demonstrate this behavior is the following example:

# sample data for linestring
df <- data.frame(
  id = c(2,2),
  lon = c(4,9),
  lat = c(9,5)
)

df <- df %>% st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  group_by(id) %>%
  summarise(geometry = st_union(geometry), do_union = FALSE) %>%
  st_cast("LINESTRING")

# function to create set of polygons
create_grid <- function(lonmin, lonmax, latmin, latmax, no_squares) {
  # a matrix defines the bounding box as required by st_make_grid()
  globe_bb <- matrix(c(lonmin,  latmax,
                   lonmax,  latmax,
                   lonmax, latmin,
                   lonmin, latmin,
                   lonmin,  latmax), byrow = TRUE, ncol = 2) %>%
  list() %>%
  st_polygon() %>%
  st_sfc(., crs = 4236)

grid <- st_make_grid(globe_bb, n = c(no_squares, no_squares),
                   crs = 4326, what = 'polygons')  %>%
  st_sf('geometry' = ., data.frame('ID' = 1:length(.)))
}
grid <- create_grid(0, 10, 1, 10, 4) %>% st_transform(4326)

# applying `sf::st_intersection` to linestring (df) and set of polygons (grid)
df_in <- df %>%
  group_by(id) %>% st_intersection(grid)

the current and expected output are explained in the following image:

Based on where the original LINESTRING starts and ends the order of the newly created LINESTRINGs is not correct since it is ordered by the names of the POLYGONs.

Is there a way to set an argument in sf::st_intersection to turn off ordering by names of polygons by default, or is there maybe another approach to this issue ?


texts: ECW support for QGIS 2.8 (Wien) and Ubuntu 14.10 (Utopic Unicorn)

texts:
Getting QGIS 2.8 (Ubuntu 14.10) ready for ECW files is really a pain. I've no terminal experience and it took me several hours to install the ECW JPEG 2000 SDK 3.3.

What I did so far:

#download ECW JPEG 2000 SDK 3.3 to /home/username
wget -nd https://api.opensuse.org/public/source/home:jluce2:GEO/libecwj/libecwj2-3.3.tar.bz2

#unpack archive
tar -xjf libecwj2-3.3.tar.bz2

#change directory
cd /home/username/libecwj2-3.3

./configure

make

sudo make install

#unstable releases of Ubuntu GIS packages
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable

sudo apt-get update

sudo apt-get install libgdal-ecw-src

Unfortunately the last command isn't working :/  It seems that there is still no libgdal-ecw-src package for newer Ubuntu versions than 13.10. Does anybody know how to fix this issue? I won't spend hours and days just for loading orthophotos.


texts: How to assign time interval to take regular cache using geowebche. For instance every afternoon at 13:00 hrs geowebcache should be started automatically and it should start creating tiles. ?


texts: Adding path to file containing photo on QT form using QGIS/PostGIS/PostgreSQL

texts:
I have just designed my first QT form.

I have a spatial table in PostgreSQL, loaded into QGIS, I have set the "photo_link" column to be "photo" in the fields/widget section of the layer's properties.  When I use the "drag and drop designer " as the layer's attribute editor, I get the expected file picker dialog and button and it all works, I see the photo above it. So I create my form in QT designer and add a "WebView" and associate it with the "photo_link" field.  I get the frame where the photo will appear and when there is already a path to the photo, it appears on the form as expected.

How do I add the path to the photo for a record where the column is null?

There is no picker mechanism on the form.  I can drag and drop a file onto the form and it displays the photo, but pressing OK and saving the edits doesn't save the path to the photo.  So how do you do it, I'm sure its something simple I've missed.


texts: Labeling at specific place using Maplex?

texts:
I have a bridge(.shp) point data which use symbol representation  and  I need to label it using field types and the label must be place under the second "ear"  like image below:

I already try using placement properties in maplex but the result is not really like i want,

Is there any way or any vb script that can i used in order to label it like i want?


texts: OpenLayers 6.x: Incorrect Static Image Location On Map -- React NextJS

texts:
VvI have been searching across the interwebs to find an answer to my problem that was answered a decade ago for a much older version of OpenLayers. I am trying to overlay a NWS Radar mosaic on top of an OSM map.
There have been plenty of complications upon the way because of the platforms I am using (mainly NextJS). I am pretty sure I resolved all the issues, but the alignment of my image still does not work.
Below is my constructor from the component that will manage the running of the OpenLayer app:
    constructor(props) {
        super(props);
        this.state = { center: [-10720000, 4550000], zoom: 4 };
        //IMPORTANT: Next.js doesn't work well with openLayer
        //SO: https://stackoverflow.com/questions/56508322/openlayers-compilation-error-with-react-next
        const KML = require("ol/format/KML").default;
        const Map = require("ol/Map").default;
        const OSM = require("ol/source/OSM").default;
        const TileLayer = require("ol/layer/Tile").default;
        //const VectorSource = require("ol/source/Vector").default;
        //const VectorLayer = require("ol/layer/Vector").default;
        const ImageLayer = require("ol/layer/Image").default;
        const Static = require("ol/source/ImageStatic").default;
        const View = require("ol/View").default;
        const Projection = require("ol/proj/Projection").default;
        const TransformExtent = require("ol/proj");
        //Establish the vector layer for the KML
        //SO: https://stackoverflow.com/questions/32135596/open-layers-2-12-and-3-kml-groundoverlay-and-screenoverlay-support
        //SO: https://gis.stackexchange.com/questions/8981/openlayers-all-image-layers-are-out-of-place-on-google-maps-v3-and-osm-base-lay/9014
        //SO: https://gis.stackexchange.com/questions/13249/reprojecting-base-layers-in-openlayers
        //SO: https://gis.stackexchange.com/questions/126979/i-cant-get-my-openlayers-radar-overlay-to-display-correctly
        //var extent = [-14206635.22, 2469859.05, -7404742.97, 6516996.45]; // [minx,miny,maxx,maxy]
        //const extent = [-14206593.415,2205071.940,-7394741.764,6517585.728];
        const extentBounds = [-14206635.218164166, 2469859.049484826, -7404742.9730659835, 6516996.454856468]; // [minx,miny,maxx,maxy]
        const extent = new TransformExtent.transformExtent(extentBounds, 'EPSG:900913', 'EPSG:3857');
        console.log(extent);
        const projection = new Projection({code: 'EPSG:3857', extent:extent});
        const radarLayer = new ImageLayer({
            source: new Static({
              url: 'https://cdn.wxyard.com/splashradar.png',
              projection: projection,
              imageExtent: extent,
              imageSize: [3400,1600]
            }),
          });

        console.log(radarLayer);
        this.olmap = new Map({
            target: null, //Needs to be null for app to initially load
            layers: [new TileLayer({source: new OSM()}), radarLayer],
            view: new View({
                center: this.state.center,
                zoom: this.state.zoom
            })
        });
        this.componentDidMount = this.componentDidMount.bind(this);
        this.updateMap = this.updateMap.bind(this);
        this.shouldComponentUpdate = this.shouldComponentUpdate.bind(this);
    }

I went and tried to use the bounds from this exchange link: Source Of Transformed Bounds without using transformExtent and sticking to EPSG:3857 and the image still didn't line up.
From what I can tell by looking at the code, transformExtent should work the same as the old bounding processes that I have commented out in the snippet.
Below is a reference to what the image looks like on the map vs an actual radar image of where the reflectivity should be.
Can you provide any ideas or guidance?
OSM Reference

Reality