texts: Is there an alternative to the Integrate tool that will not change one of the inputs

texts:
I have two polygon layers, and would like to snap the outlines of one of them (colored poly's) to the other (grey) without changing the grey in anyway (so Integrate seems out). It should only snap if the edge of the colored poly's is within a distance of the grey poly edges already, something like 200ft. There are several thousand of them that need snapped so it can't be a by hand method. Am I forgetting some tool that could do that since Integrate changes both inputs? I have ArcInfo 10.1 most if not all extensions and Xtools Pro.

Thanks.


texts: In what ways has QGIS 2.x API changed making Plugin Builder Tutorial invalid?

texts: I am working on the QGIS workshop Tutorial for QGIS Plugin Builder and notice that it is for QGIS 1.1. As this seems like a popular tutorial, in general, what are some areas of this tutorial that need to be changed since the QGIS 2.0 release?


texts: Using Toggle PopUp Functionality in ArcGIS API for JavaScript?

texts:
I've look through every example on Stack GIS.Stackexchange which might deal with my problem. Of course there was good input, but unfortunately I'm not able to fix it.

Admittedly my JavaScript skills are not the best, so "easy" thing are also very hard to overcome.

Until now I've stuck at the following level: http://jsfiddle.net/84vvR/

I'm using nliu's TOC http://gmaps-utility-gis.googlecode.com/svn/tags/agsjs/latest/examples/toc.html and the editor and popup functionalities overed by esri's API.

As you can see, theoretically all work. The mapapp contains an editor an the option for popup.

However, if the user wants to add new features by clicking on the map, the popup pops out and gives information about the underlying area, which has the defined popup. The same problem, if the user want to use the measurement tool. On the first click, the popup pops out and so it's not possible to add a starting point for measurement.

I've tried to implement a workaround with dojo.connect and a clickHandler or onClick like (button dojoType=”dijit.form.Button” onClick=”....(),...();”>Disable Popup Function

I assume that it's not too hard for an advanced developer to solve this problem to bind the Popup Functionality to a checkbox or Button.


texts:
In ArcGIS 10.2 ArcMap am attempting to label a grid of varying size polygons overlaid by a grid of evenly spaced (1km) points with an associated value X by the the value of the point. If a polygon is intersected by more than 1 then I want it to be labelled with the maximum X value. If a polygon doesn't intersect with a point I want it to be labelled with the nearest point.

I used spatial join with one to one intersect and maximum value for X to be selected. This allowed me to get the maximum value for multiple points that intersected. However it doesn't allow me to label the polygons where a point does not intersect it. I put a search radius on  it and then it can label polygons with points where X is greater but are not contained in the actual polygon.


texts: How can I geocode CSV by postcode, using existing point data, in QGIS?

texts:
I have a very simple MapInfo workflow for geocoding a CSV file which contains UK postcodes, based on an SQL select query referring to point data for postcodes.

This is simple. Using Mapinfo open the CSV, open the postcode layer, run a single select query which combines according to the matching column and creates a temporary table (combining CSV data with the Postcode data AND geometry), save the table.

I normally work in QGIS. I'd like to be able to omit Mapinfo from this workflow, but I cannot find any sensible way to achieve this, and certainly nothing so simple (despite already using table joins, relationships, filters, and all sorts of other techniques elsewhere). I'm sure that there's a process which involves a Spatialite database, but I'm out of my depth at this point, and it's also going to be a process with multiple (relatively involved) steps too.

Am I missing something?


texts: QGIS Layout: Locking the scale of loaded map even though the location changes

texts:
Curently I have made some templates of layouts where there just needs to be set the map theme and the location that can be panned to a "favorite" location. The problem annoyance is of course when you pan to that location, we have reset the scale back to what it was intended. Now my Question is there somewhere an option to "lock" the scale so that even though we pan to this favorite location, the scale does not change?
This favorite location is more meant to "center" the map on a project area, making the project area be in the center of the map.


texts: Request a Ajax with GET to remove one or more geometries in OpenLayers

texts:
I have a html selection of the form:

          Afiseaza:

          <select id="layer" onchange="search()" class="form-control">
               <option value="Green Area">Teren</option>
               <option value="Building">Constructii</option>

          </select>

Here we have a function in js which on an onchange event calls for the following:

function search(){

    var layer = $('#layer').val();

    console.log(layer);

      proj4.defs('EPSG:31700','+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000 +ellps=krass +towgs84=28,-121,-77,0,0,0,0 +units=m +no_defs');
      ol.proj.proj4.register(proj4);

//var geojson = $.getJSON('json/geojson.json', function(data) { return data });
//var geojson1 = {"type":"FeatureCollection","features":[{"type":"Feature","id":"2","geometry":{"type":"LineString","coordinates":[[698282.236,377375.758],[698289.591,377372.136],[698183.88,377182.298],[698182.72,377197.045],[698282.236,377375.758]]},"crs":{"type":"EPSG","properties":{"code":"31700"}},"geometry_name":"wkb_geometry","properties":{"ogc_fid":"2","area_id":"2","building_id":"2","type":"test","area":"test","owner":"test","description":"test","administration":"test","address":"test","property":"test","layer":"Green Area","extendedentity":null,"linetype":"test","entityhandle":null,"text":"test","wkb_geometry":"01020000000500000000FED478544F2541002731087F081741408D972E634F25418195438B70081741295C8FC28F4E254181E92631790517410AD7A3708D4E2541027B142EB405174100FED478544F2541002731087F081741"}}]}

//in cazul in care dorim direct din GeoJson
//$.getJSON('json/geojson.json', processGeoJSON );

  //console.log(geojson);

  if (ol.proj.proj4 && ol.proj.proj4.register) { ol.proj.proj4.register(proj4); }

    var features;

        $.ajax({
        type: "GET",
        url: "get_geojson.php?layer="+layer,
        dataType: 'json',
        async: false,
        success: function (response) {
           //console.log(response);
           //console.log(b);
           features = new ol.format.GeoJSON().readFeatures(response);
        }
    });

  var extent = ol.extent.createEmpty();
  features.forEach(function(feature){
    feature.setId(undefined);
    feature.getGeometry().transform('EPSG:31700', 'EPSG:3857');
    ol.extent.extend(extent, feature.getGeometry().getExtent());
  });

  var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM(),
        }),
        new ol.layer.Vector({
            source: new ol.source.Vector({
                features: features
            })
        })
    ],
    target: 'map',
    controls: ol.control.defaults({
            attributionOptions:
                ({
                    collapsible: false
                })
        }).extend([
            new ol.control.OverviewMap(),
            new ol.control.FullScreen(),
            new ol.control.ScaleLine()

        ]),
        interactions: [
            new ol.interaction.DragRotateAndZoom(),
            new ol.interaction.MouseWheelZoom(),
            new ol.interaction.DragPan(),
            new ol.interaction.DoubleClickZoom()
        ],
    view: new ol.View()
  });

  map.getView().fit(extent);

}

Which would normally have to refresh the map with new geometries ...

This does not happen does anyone know why this phenomenon is happening?

In the console there is no error, just when the page load is executing the function but when I give the onchange do not ...


texts:
Is there a direct way to load a CSV as an sf? I've looked through the sf documentation and nothing jumps out at me, currently using:

my_csv <- read.csv('my_data.csv', header = TRUE)

my_sp <- sp::SpatialPointsDataFrame(coords = data.frame(my_csv$y,
my_csv$x), data = data.frame(my_csv$data), proj4string = CRS(BNG))

my_sf <- sf::st_as_sf(my_sp)


texts:
I have a polygon shapefile layer. I want to add a rectangle polygon (just like in Arcmap), and another one using the previous polygon edge. I need it to be accurate, my hands do not meet this condition.

Maybe this draw will explain better:


texts:
I would like to reproduce this map which was originally printed in mare, a German magazine.

I have done some research and think it is a projection by Athelstan Spilhaus. I would like to be able to produce it using QGIS version 3. This is an interesting article about the projection: https://le-cartographe.net/blog/archives/342-athelstan-spilhaus   And this is another interesting page from Jason Davies: https://www.jasondavies.com/maps/spilhaus/11/

I have tried with an antarctic polar stereographic projection (EPSG 3031) where I changed the standard parallel to various different lines of latitude but this does not show the north pole correctly, nor does it have the coastline around the edge. I have tried using a lambert azimuthal equal area projection with the latitude or origin at 0, but this doesn't get close to it either.
It seems to be a mix between a Hammer and Petermann projection, according to the french article above.

How would I make this in QGIS? Is it possible? (I have a feeling that the answer is no, but I wanted to ask others).


texts: Inserting new features to PG\PostGIS table with generated column using QGIS?

texts:
I'm using PG 12.2 with PostGIS 3.0.0, and I've created a table with a generated column:

CREATE TABLE public.test1
(
    id integer NOT NULL DEFAULT nextval('test1_id_seq1'::regclass),
    geom geometry(Polygon,3857),
    area_sqrm double precision GENERATED ALWAYS AS (st_area(geom)) STORED,
    CONSTRAINT test1_pkey PRIMARY KEY (id)
)

I want to create new features in this table using QGIS 3.12.2. Therefore I did the following steps:

Import table into QGIS project
Toggled editing
Sketched new valid geometry polygon entity
Approved QGIS default form

Saved layer edits

This last step resulted with the following error:

Could not commit changes to layer test1

Errors: ERROR: 1 feature(s) not added.

  Provider errors:
      PostGIS error while adding features: ERROR:  column "geom" does not exist
    LINE 1: ...") VALUES (st_geomfromwkb($1::bytea,3857),st_area(geom)) RET...
                                                                 ^
    HINT:  There is a column named "geom" in table "test1", but it cannot be referenced from this part of the query.

After omitting the default value QGIS applied for the area_sqrm column, I tried to save again. This resulted with the following error:

Could not commit changes to layer test1

Errors: ERROR: 1 feature(s) not added.

  Provider errors:
      PostGIS error while adding features: ERROR:  cannot insert into column "area_sqrm"
    DETAIL:  Column "area_sqrm" is a generated column.

How can I create new features using QGIS, and let their generated columns to be populated DB-wise?

In this case, it can be thought as working with ArcGIS' FGDB feature class which has two dedicated computed geometry fields: shape_area; shape_length. In its environment, the two fields are being computed after creation and altered after the entity's geometry is edited.


texts: How to export dxf of window view from postGis MasterMap layer

texts:
I have managed to get qgis to show my Mastermap data which is stored in a postGis database. It's national coverage so seems to take ages to load but that's another question.

I basically want to export the map that I can see (or at least one layer of it) as a DXF.

When I try to do this (layer save as) it starts to generate a temp version of the dxf file but then just hangs and I have to use task manager to close it and start all over again. It's as if it's trying to read the whole database.

The postGis database was created using Edonica. Some of the field names are different to those you get if you load a gml file directly in to qgis so this may be the problem. Or do I have to   make a selection of the area I want?

Thanks for any help

Graham


texts: ArcGIS Online WebMap authentication timeout

texts:
I have an ArcGIS Online public account and add WebMap to my website.

My ArcGIS Online WebMap looks like this ESRI's sample: LINK

And I am trying to add my WebMap to my website like this ESRI's reference page. You will see there is a map in the center of page: LINK

My WebMap is displayed on my webpage well. When I access my webpage, my WebMap asks my ID and Password. If I entered it, then it shows my map.

However, my question is, if I moved to different page and then come back to map page, it asks again. Is it possible to set a timeout so I don't have to sign in everytime I access the page?


texts: Opening table from file geodatabase in ArcPy?

texts:
I have tables in a file geodatabase. I want to open and work with them in my script. The first step I appear to have to do is scripted underneath. My result of this script is all names of all tables in my Geodatabase. I would like to also extract all data within these tables. So that I have a matrix or a list of lists in python.

Is this possible?

import arcpy
from arcpy import env

env.workspace = r"D:\data\sre20821\Documents\Voorspelmodel\Geodatabases\2_Component.gdb"

datasetList = arcpy.ListTables("*")

for dataset in datasetList:
     print dataset

http://resources.arcgis.com/en/help/main/10.2/index.html#//002z00000011000000


texts: How to sort layers in TOC in custom order?

texts:
Recently I saw a script that moves layers according to their names in TOC.
I wanted to make a complete script for my work in which all needed layers are listed and probably should be sorted in a way I want.

Here it is:

import arcpy
mxd=arcpy.mapping.MapDocument("CURRENT")
df= arcpy.mapping.ListDataFrames(mxd,"")[0]
from arcpy import mapping
for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "ramka":
        moveLayer = lyr
    if lyr.name == "address":
        refLayer = lyr
mapping.MoveLayer(df,refLayer,moveLayer,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "address":
        moveLayer1 = lyr
    if lyr.name == "river":
        refLayer1 = lyr
mapping.MoveLayer(df,refLayer1,moveLayer1,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "river":
        moveLayer2 = lyr
    if lyr.name == "allhouse":
        refLayer2 = lyr
mapping.MoveLayer(df,refLayer2,moveLayer2,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "allhouse":
        moveLayer3 = lyr
    if lyr.name == "topo":
        refLayer3 = lyr
mapping.MoveLayer(df,refLayer3,moveLayer3,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "topo":
        moveLayer4 = lyr
    if lyr.name == "village":
        refLayer4 = lyr
mapping.MoveLayer(df,refLayer4,moveLayer4,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "village":
        moveLayer5 = lyr
    if lyr.name == "quarter":
        refLayer5 = lyr
mapping.MoveLayer(df,refLayer5,moveLayer5,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "quarter":
        moveLayer6 = lyr
    if lyr.name == "town":
        refLayer6 = lyr
mapping.MoveLayer(df,refLayer6,moveLayer6,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "town":
        moveLayer7 = lyr
    if lyr.name == "district_t":
        refLayer7 = lyr
mapping.MoveLayer(df,refLayer7,moveLayer7,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "district_t":
        moveLayer8 = lyr
    if lyr.name == "district_a1":
        refLayer8 = lyr
mapping.MoveLayer(df,refLayer8,moveLayer8,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "district_a1":
        moveLayer9 = lyr
    if lyr.name == "district_a":
        refLayer9 = lyr
mapping.MoveLayer(df,refLayer9,moveLayer9,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "district_a":
        moveLayer10 = lyr
    if lyr.name == "admin2":
        refLayer10 = lyr
mapping.MoveLayer(df,refLayer10,moveLayer10,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "admin2":
        moveLayer11 = lyr
    if lyr.name == "admin1":
        refLayer11 = lyr
mapping.MoveLayer(df,refLayer11,moveLayer11,"BEFORE")

for lyr in mapping.ListLayers(mxd,"",df):
    if lyr.name == "admin1":
        moveLayer12 = lyr
    if lyr.name == "admin":
        refLayer12 = lyr
mapping.MoveLayer(df,refLayer12,moveLayer12,"BEFORE")

arcpy.RefreshTOC()
arcpy.RefreshActiveView()

Unfortunately, testing of this script in other databases throws an error like

Runtime error
Traceback (most recent call last):
  File "<string>", line 39, in <module>
NameError: name 'refLayer4' is not defined

It was because of one of the layers did not exist in TOC.
So what would you recommend if it is needed to have an order in TOC like in script despite the absence of some layers.
Or is there any other script sample that sorts layers?


texts: ArcGIS10.0 (32 bit) waits for splwow64.exe indefinitely and stops responding

texts: I am running ArcMap 10.0 (32 bit) on Windows 7 (64 bit), while using some extensions the application waits indefinitely for splwow64.exe and then either crashes or stops responding, shall I try running application on a 32 bit system ?


texts: How to reduce area boundaries resolution?

texts:
I managed to plot zip codes boundaries for an entire US state using ArcGIS data and python tools such as arcpy and Bokeh. The result:

The dataset has high resolution area boundaries so, I have like 2000 data points for each polygon and there are 1700 polygons. Thus, rendering 3.4 million data points approx. is very slow and the interface lags.

The questions is, do you know any methodology to reduce the area boundary resolution so the number of data points to be plotted is reduced? The resulting polygons should not overlap at all. Any other alternatives?


texts: Python script fails on stand alone but works in ArcMap as a toolbox

texts:
I have a fresh install of ArcGIS Server and ArcGIS for Desktop both 10.2, latest patches applied. The data is served in this application from a PostgresSQL enterprise geodatabase.
The problem is that I can connect to the database using ArcMap, on the ArcGIS Server (the maps that have their data source set in the database) and I can run python scripts that interact with the database in the python console in ArcMap. Yet, the same scripts, or even the smallest interactions (e.g. listing the feature classes with arcpy.ListFeatureClasses()) does not work stand alone in python command line. I have tried both the 32 bit and 64 bit Python installations that came with the server and desktop and none run any smallest interactions with the geodatabase.
Any thoughts?

For reference, a small snippet of code that fails:

import arcpy
arcpy.env.workspace = 'path/connection.sde'
arcpy.ListFeatureClasses()

returns:

[]


texts: Changing value for features overlapping another feature in the same layer

texts:
In a single layer, I want to use the feature with the value '33500' in the field "Seeding_Ra" to override the value of any feature that overlaps it, to that value, but just within the overlaping geometry.
Is there a way in the vector tools to do this?


texts:
I need to use this file on a map but I am not able to match the projection with the others I have, primarily data downloaded from OSM (city of Boston). I read many articles, I have repeatedly tried to change the layer CRS, convert to Shapefile but to no avail. Can anybody help?

This is the link to the dxf file:

http://www.bostonredevelopmentauthority.org/documents/planning/udtg/2d/city-blocks


texts: Packaging geospatial libraries for cross-platform application, or using docker container?

texts:
I am developing some python tools to help some less technical users process geospatial imagery. They generally need to do simple stuff like compress images, reproject images, etc. The users will be on Windows, Mac, and Linux. I use libraries like rasterio, shapely and geopandas in the background.
So I was able to develop the package in conda-forge and setup CI testing for all platforms, and I thought that was pretty good. But then I started to find little things, like one of my alpha-testers on Mac was able to execute compression using rasterio, but the resulting file was stripping out the CRS for the file. This problem only occurred on Mac, and not on the other platforms. There were other things like the webp compression standard was not available for Windows, but it was for linux and mac.
My sense is that these little differences between platform is due to some dependency mismatches. Now even though I use conda-forge package when installing both my own package and the dependencies, this still might not be enough control.
So my alternative idea is to create a "rootless" docker container that holds the tool, and then create a simple command line application that just forwards all processing to the docker container. Then I could just make a very simple pip installable application and have the user simple download the docker software.
So in deciding whether to stick with the original conda-forge approach, or to go with creating a docker container, I was hoping someone could answer a simple question. Should I expect there to be a lot more of these small cross-platform nit-picky things that I will have to spend a lot of time resolving? If that is the case then having a docker container makes sense since I can then focus on adding tools instead of debugging. But if these types of cross-platform issues are relatively easy to handle, then I can just add more test cases and try to resolve on a case by case basis.
I was hoping that someone with more experience maintaining some cross-platform libraries might be able to answer this question. Thanks.


texts: Changing value field when applying symbology from layer using Python Add-in?

texts:
I have a python addin button that creates a layer file from a polygon feature class in ArcMap. I want to take that layer file and, with user input from a combobox, apply it back to the feature class using a different value field. For example, I have a feature class called "names" with the following attributes:

And I want to create a layer file based on the unique values for Name and then use that layer file to display Name_Change1 and Name_Change2 so that my colors remain consistent for each unique value (ex: Finnley is always purple, Dylan is always green).
This is possible in ArcMap by importing the symbology from the layer file and changing the value field. However, I need this functionality in a python toolbar.

Right now, I am creating separate layers for each field and using those. This works but it requires a lot of set-up in order to use the addin:

if combobox_value == "Name_Change1":
    if arcpy.Exists("Name_Change1.lyr"):
        arcpy.ApplySymbologyFromLayer_management(names, "Name_Change1.lyr")

But I want it to be something like:

if combobox_value == "Name_Change1"
    if arcpy.Exists("Name.lyr"):
        #assign/change value field of Name.lyr
        arcpy.ApplySymbologyFromLayer_management(names, "Name.lyr")

Where symbology is always pulled from "Name.lyr".

I guess I really wish that the ApplySymbologyFromLayer tool had a parameter for the value field, but that doesn't seem to be the case. And I haven't found another tool/method that will allow this. Does anyone have any ideas?


texts: What is the purpose of PostGIS on PostgreSQL?

texts:
PostgreSQL already supports spatial data types, operators and indexing.

What does PostGIS provide exactly that made it necessary to exist as an extension to PostgreSQL?

Why don't we all just use PostgreSQL's spatial functionality?


texts: Getting message if shapefile does not exist in folder when executing Merge process?

texts:
I´m following this link Merge shapefiles in folders and subfolders with arcpy very very close because it really helps me to solve my problem. However, I rebuilt the script to accomplish my job but when I run the script to merge some shapefile, I get the next error:

this my rebuilt script:

import arcpy, os, sys, string, fnmatch
import arcpy.mapping
from arcpy import env

rootpath = r'F:\My_Shps\My_zone'
shp = 'mojonera50_p.shp'
counter = 0
matches = []
for root, dirs, files in os.walk(rootpath):
    for shps in files:
        if shps == 'mojonera50_p.shp':
           match =(os.path.join(root, shps))
           matches.append (match)
           counter = counter +1

arcpy.Merge_management(matches,r'F:\My_Shps\My_zone\Result')
print "Se unio la capa {0}".format(shp)

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid. ERROR 000735: Input Datasets:
Value is required Failed to execute (Merge).

I think when I running the script it doesn't find the shapefile, but how Can I modify the script to get a simple message warning that shapefile is not in the rootpath


texts: Framing argument for insert cursor when writing geometries using ArcPy?

texts:
I am writing polygon geometries and I have points holding the geometries that represent vertices of the polygon.

When I run the code below I get:

  TypeError: sequence size must match size of the row

cur = arcpy.da.InsertCursor(fc_path, ["SHAPE@"])

# Stuff
    X_Min = DF_extent.XMin
    X_Max = DF_extent.XMax
    Y_Min = DF_extent.YMin
    Y_Max = DF_extent.YMax

    point_0 = arcpy.Point(X_Min, Y_Min)
    point_1 = arcpy.Point(X_Min, Y_Max)
    point_2 = arcpy.Point(X_Max, Y_Max)
    point_3 = arcpy.Point(X_Max, Y_Min)

   # this is the part that is causing errors:
    cur.insertRow((point_0, point_1, point_2, point_3, point_1))

When I have the vertices that make up a polygon, what is the correct way to pass that information to the insert cursor so it knows that all of those points are, in fact, one argument to be taken at once?


texts:
We are using OpenLayers for our GIS Application and we have base files from Ordinance Survey for UK maps.

For creating WMS layers, we are generating tiles of 256 x 256 through Sharpmap Library but facing issues when the base layer is Google or Bing.

We know Google uses EPSG:900913 projection and we are able to transform the projection from EPSG:900913 to OSGB36 successfully. For this we are using GeoTools.js methods to achieve the transformation.

But the issue is that the tiles being created are not projected properly. There is a slight misalignment at the connecting point of each tile with other tile.


texts:
I have a folder of a few thousand GeoTIFF files of several years of satellite data that I wish to upload into my assets folder on the Google Earth Engine.
The only way I know how to upload the data is to manually press 'upload' each time and add the time data to the file.
Obviously, this is not ideal for my current circumstances.
Is there any way to batch upload the files while adding the 'system: start time' data?


texts:
I created a raster with the function rasterFromXYZ() in R using a dataframe with coordinates and the value of a variable, the result was this:

plot(raster1)

which is exactly what I need. Now I just want to convert/export it exactly like that as a shapefile.
I guess it should be easy but I haven't found the answer... Somebody?


texts: Metadata server and spatial server - a must?

texts:
I have been using ArcMap (only) for the past two years (read: newbie) in the oil and gas (upstream) industry.
I currently have about 200 GB of files (shp, tif, grids, a few fgdb, all sitting in a windows file server) in total, spread over 10 GIS projects. each project folder have their own files and mxd as they are of different regions.

i only learned abt this site a couple of days back and i kept reading about metadata server and spatial server.

I would like to hear from experienced GIS users/administrators here, what are the capabilities /functionalities of metadata servers and spatial server that ArcMap/ArcCatalog are not capable of? Is adding data to mxd / data loading using metadata one of them?

from what i read here, i can search for files much easier with a metadata server than AcrCatalog, and why is that so?

also, can anyone point me to a site that is like a "GIS infrastructure for dummies"?

thanks.


texts: ArcGIS Polygon to Raster cell size issue

texts:
I am trying to convert a polygon to raster and I noticed that the cell size of the resulting raster appears different for different coordinate systems, despite me setting the cell size the same. Both coordinate systems are in meters. Here's what I mean:

To begin, I set the polygon to raster cell size as 28 meters.

I run the tool with projected coordinate system of "WGS 1984 Web Mercator (auxiliary sphere)" and the resulting cell size is about 20.8 meters:

What's even weirder is that when I check the properties of that layer, it still says that the cell size is 28 and that the units are meters:

Now if I run the polygon to raster tool with the projected coordinate system of "NAD 1983 UTM Zone 16N", the cell size is as I expect, 28 meters.

Can anyone explain why the resulting cell sizes are different?


texts: How to convert Mapinfo “.Map” file (without .TAB) into SHP?

texts:
I've got several MapInfo files but i don't have the .TAB.

I only have:

M044.DAT
M044.ID
M044.IND
M044.MAP

I need to convert these files to SHP.

I tried to read those files with MapInfo Proviewer 11 and QGIS 2.2 with no success...

Regards,


texts: Possible to get Lat/Lon coordinates of triangle vertex point?

texts:
I have two known points (A and B) on latitude / longitude coordinates.
Certainly I know the distance between A and B in meters.

I have third point C with unknown coordinates, but I know the distance between AC and BC in meters.

Is it possible to get the lat/lon coordinates of point C with math formula?

Here is an example:

The all known distances between coordinate points are:

AB = 120.9234 meters
AC = 226.7206 meters
BC = 213.1376 meters

The known coordinates of A and B point:

A lat = 47.49784;
A lon = 19.08939;

B lat = 47.49893;
B lon = 19.08941;

and I calculated all angles of triangle in deg:
B: 80.341932306607
A: 31.722055823679
C: 67.936011869714

and I measured the result what I looking for in map:
Lat: 47.498575374822
Lon: 19.09129858017

but I like to use a math formula to get these two coordinates of point C using the upper input values.


texts: How to open the Electronic nautical chart S-57 data?

texts: Anyone know how to open the Electronic nautical chart S-57 data by using GIS software e.g. Arcview, QGIS? and converted it to general GIS / CAD formats.


texts: ArcGIS Mobile projects access from Field device

texts:
I am trying to figure out how to allow field users to connect to ArcGIS Server based mobile services within an enterprise environment. As part of this, I am asking network security to allow this access from Trimble GPS units. We are going to be doing this for the first time and data sensitivity is an issue. I am looking for some feedback from other users with similar experience within enterprise environment.

I know that I need to allow access to mobile projects but I am curious if there is any other access that we need to provide. I want to have as much detail as possible for network security team.
*http://ArcGIS/Mobile/Content/WindowsMobile/Projects

Appreciate any input.

Thanks


texts: Calculating beta-diversity w raster data [renamed]

texts:
Using ArcGIS 10, it is not possible to generate a raster attribute table for 32-bit float values. Instead, I have tried the RasterToNumPyArray method, which crashes since the file input raster is so large (I'm talking >7GB).

I've also tried converting the raster into an ASCII text file, and from there reading the file, line by line, into numpy arrays:

e.g.

import numpy as np
file = open(r"c:/yourfile.txt", 'r')
array_list = []
for line in file:
    array = np.loadfromstring(line.strip(), dtype = np.float32, sep = " ")
    array_list.append(array)

My main issue is that this is still pretty slow since the file has 71058 columns and 27820 rows.

Any advice? The troubling thing is that I can generate a list of unique values in ArcMap from symbolization, but there's no easy real way to move them into a table (besides copy and paste of each value).

Many thanks!


texts:
I am working on a GeoServer API, which on request provides a WMS image in response.

My requirement is for the monitoring of the API hits, the number of objects displayed on the WMS response will be considered as the number of transaction or API hits.

But in the development of the above scenario, I am facing an issue as to how will I monitor the transaction hits of my clients who will be accessing my API?


texts: Specific div elem for click event within view.when/view.on

texts:
view.when(function(){
    view.on("click", function () { doSomething = !doSomething;});
});

The above works as a global click event. HOWEVER, I am trying to specify a specific div elem not a global elem. i.e. The logic I'd like:

view.when(function(){
    view.on("click", #myDivElm function () { doSomething = !doSomething;});
});


texts: How to disable the Escape key for a dialog?

texts:
I may have found the solution to my own question but I have another problem. It turns out that the user can just hit "Esc" and close the form. How can I disable the escape key? I'm looking for some python code, but I'm at a loss.

Edit

Below is my code, if it helps.

# -*- coding: iso-8859-15 -*-
from PyQt4 import QtCore, QtGui
import psycopg2
import hashlib

def authentication():
    conn_string = "host='host' dbname='dbname' user='username' password=pass"
    conn = psycopg2.connect(conn_string)
    global userText
    username = str(userText.text())
    global passText
    hashpassword = hashlib.md5(str(passText.text())).hexdigest()
    cursor = conn.cursor()
    cursor.execute('SELECT * FROM users WHERE username=%s AND password=%s', (username, hashpassword))
    if cursor.rowcount == 1:
        cursor.close ()
        conn.close ()
        global myDialog
        myDialog.accept()
    else:
        cursor.close ()
        conn.close ()
        msgBox = QtGui.QMessageBox()
        msgBox.setWindowTitle(QtGui.QApplication.translate("loginform", "Σφάλμα", None, QtGui.QApplication.UnicodeUTF8))
        msgBox.setText(QtGui.QApplication.translate("msgBox", "Λάθος Όνομα Χρήστη ή Συνθηματικό", None, QtGui.QApplication.UnicodeUTF8))
        msgBox.setWindowFlags(QtCore.Qt.WindowStaysOnTopHint);
        msgBox.exec_()

class Login (QtGui.QDialog):
    dlg = QtGui.QDialog()
    dlg.setObjectName("loginform")
    dlg.setWindowTitle(QtGui.QApplication.translate("dlg", "Φόρμα σύνδεσης", None, QtGui.QApplication.UnicodeUTF8))
    dlg.setEnabled(True)
    dlg.resize(400, 203)
    dlg.setMaximumSize(QtCore.QSize(400, 203))
    dlg.setSizeGripEnabled(False)
    dlg.setModal(True)
    dlg.setWindowFlags(QtCore.Qt.FramelessWindowHint);
    event = QtGui.QKeyEvent(QtCore.QEvent.KeyPress, QtCore.Qt.Key_Escape, QtCore.Qt.NoModifier)

    font = QtGui.QFont()
    font.setPointSize(10)

    username = QtGui.QLineEdit(dlg)
    username.setEnabled(True)
    username.setGeometry(QtCore.QRect(190, 19, 151, 21))
    username.setFont(font)
    username.setObjectName("username")
    global userText
    userText = username
    userlabel = QtGui.QLabel(dlg)
    userlabel.setEnabled(True)
    userlabel.setGeometry(QtCore.QRect(70, 20, 91, 16))
    userlabel.setFont(font)
    userlabel.setLayoutDirection(QtCore.Qt.LeftToRight)
    userlabel.setObjectName("userlabel")
    userlabel.setText(QtGui.QApplication.translate("dlg", "Όνομα Χρήστη", None, QtGui.QApplication.UnicodeUTF8))
    userlabel.setBuddy(username)

    password = QtGui.QLineEdit(dlg)
    password.setGeometry(QtCore.QRect(190, 59, 151, 21))
    password.setFont(font)
    password.setEchoMode(QtGui.QLineEdit.Password)
    password.setPlaceholderText("")
    password.setObjectName("password")
    global passText
    passText = password
    passlabel = QtGui.QLabel(dlg)
    passlabel.setGeometry(QtCore.QRect(70, 60, 91, 16))
    passlabel.setFont(font)
    passlabel.setMargin(0)
    passlabel.setObjectName("passlabel")
    passlabel.setText(QtGui.QApplication.translate("loginform", "Συνθηματικό", None, QtGui.QApplication.UnicodeUTF8))
    passlabel.setBuddy(password)

    btnLogin = QtGui.QPushButton(dlg)
    btnLogin.setGeometry(QtCore.QRect(160, 170, 75, 23))
    btnLogin.setDefault(False)
    btnLogin.setFlat(False)
    btnLogin.setObjectName("btnLogin")
    btnLogin.setText(QtGui.QApplication.translate("loginform", "Σύνδεση", None, QtGui.QApplication.UnicodeUTF8))
    btnLogin.clicked.connect(authentication)

    global myDialog
    myDialog = dlg

    def keyPressEvent(self, event):
        msg = QtGui.QMessageBox()
        msg.setText("You Pressed: "+ event.text())
        msg.setWindowFlags(QtCore.Qt.WindowStaysOnTopHint);
        msg.exec_()

    dlg.show()


texts: Does installing next QGIS 'build' require uninstalling previous version?

texts:
I am using a Mac OS-10.10.5 ("Yosemite") and have been running QGIS 2.10 "Pisa". Several people have reported issues about processing rasters (http://hub.qgis.org/issues/13162), so I am chomping at the bit to install QGIS 2.12 (supposed to be coming out soon!).

Before I install the new version (QGIS 2.12), do I need to uninstall the old 2.10 version?


texts: how do you erase polygons that don't have a line intersection?

texts:
I am trying to erase any polygons that do not have a line intersection through them and am unsure of what tool to use
using ArcMap 10.1


texts: JTS Vividsolutions API: compute First (direct) geodetic problem

texts:
I'm using the popular JTS Vividsolutions API for Java. I'm struggling to find if I can perform calculation of the First (direct) geodetic problem with the API.
Is there a way to calculate this with this API? It seems that simple one method not exist, but maybe using some combinations of existing methods of that API.

Note - I'm not looking for solutions of implementing the First (direct) geodetic problem algorithm using Java+JTS suite.


texts:
I'm working with a polygon layer in ArcMap (version 10.3) that's currently symbolized as graduated symbols (in this case, teal dots on the map). I'd like to set this layer as slightly transparent to see underlying layers, but when I set a transparent percentage under the Display tab in the layer properties, nothing happens. Anyone know if it's possible to set a transparency for this kind of symbology?

It seems that the graduated symbols symbology has other head-scratching properties, like the fact that it always draws as the top layer, even if it's buried in the TOC. See screenshot below of the layer in the TOC that shows transparency buy the actual layer in the map window is definitely not transparent.


texts: So I am currently creating a land cover map of all of PA. All of my counties will display except for one county. It doesn't even appear in ArcCatalog when the preview of the raster is selected. I reran maximum likelihood with the same signature file as before and got it to appear for a split second and then it went away and won't reappear. When I google this it says to adjust the display quality but that doesn't help. There is also no scale extent either. The only time it will appear is when I'm zoomed in past 1:8000 but again, I don't have a minimum or maximum scale extent set. Was wondering if anyone else has had this problem?


texts: Removing polygon outline that appears when creating multiple maps with ArcPy?

texts:
I am trying to create a series of maps using arcpy (see below my code that creates one map in each iteration). For each map I need to change the definition query and the value field, keeping the same classification break values. The script does that all, but polygon outline appears (see picture on the left).

Does anyone know how to remove the outline or how to perform the same task without having the outline appear (as in the picture on the right)?

I am using ArcGIS for Desktop 10.2 Basic.  My ArcMap Document contains only polygon layers. The polygons are small squares. They are not a raster, though.  I first symbolized the layers manually (Quantities, Graduated Colors, Classification: Natural breaks). The break values from this classification are in list1 in the script below, as I want all maps to have these break values, but in each map different field should be classified. That is why I change the value field in my code, and then use list1 to define break values. However, this causes the polygon outline to appear.

symbologyTemplate = arcpy.GetParameterAsText(1)
print "START"
mxd = arcpy.mapping.MapDocument("path_map.mxd")

prnw_lyr = None
for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.name == "Test_lyr":
        prnw_lyr = lyr
        break

list1 = [10, 30, 50, 70, 100]

for i in range(2, 10, 1):
    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.name == "Test_lyr":
            queryStr = "nrp_" + str(i)
            lyr.definitionQuery = queryStr + "> 0"
            print "new atrribute " + queryStr + " done"
            prnw_lyr.symbology.valueField = "prnw_" + str(i)
            lyr.symbology.classBreakValues =  list1
            lyr.symbology.classBreakLabels = ["Label 1", "Label 2", "Label 3", "Label 4"]

    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToJPEG(mxd, "path_export" + str(i) + ".jpeg")

del mxd
print "Process DONE"


texts: Adding a point at given lat/long Australia - suspected CRS issues

texts:
I am hoping to add a geographical point to a layer in QGIS with a specific latitude and longitude value (from GPS.) I have had no success with a number of methods (trying to edit with Numerical Vertex Edit, and by uploading a csv.) My project is in CRS GDA 94 MGA zone 55, EPSG 28355. "On the fly" is selected. I used the same CRS for my data point when importing the CSV. X field is long, Y field is lat. However, the point is not placed in the right place - it goes off way distant from my other data points, I think in a geographical point that does not even exist on earth.

When I try to use the Numerical Vertex Edit method, I have found no combination of settings that give me the point data X and Y values in Lat Long - they all seem to be in meters.

I am not totally across the difference between CRS systems - but am usually fine by saving all layers and projects as GDA 94 zone 55.


texts: Creating new raster layer based on operations in another layer in ArcMap

texts:
I am using ArcMap 10.7.
I am to do a multi criteria analysis on a tract of land to determine suitable areas for housing development. One of the criteria is "amount of hectare of green land within 1km".
I have the green land data in a vector layer.
I would like to create a raster layer that calculates the amount of hectares of green land within 1km from the cell center. Perhaps (or probably) based off another raster layer with information about the amount of green land.
My skill level is that of a beginner. This is for an assignment where the purpose is to learn how to use ModelBuilder, calculate distance, reclassify and use weighted overlay. It's my first time working with raster analysis.
I already have a simplified solution to my problem (the suggestion is to simply filter out smaller polygons and use euclidean distance from any green land), but out of curiosity/interest are there more ways to solve this issue (just clarifying that this isn't about doing my school work for me)?


texts: Duplicating layer, copies all records, NOT JUST selected records using QGIS

texts:
The requirement is to display only selected records as a New Layer pointing to same source data on disk. If I were selecting using attributes of the current layer it would be simple - I use the filter option. But, I have selected a set of records using a Select by Location. It is these set of records I wish to display as a new layer. Duplicate a layer, copies all records, NOT JUST the selected records into the new layer. I am using QGIS 3.10.6.
OR if there could be a way to get the set of selected records and set it as a filter, this would also do.


texts: Selecting all cells not containing a number in SQL

texts:
I have a field containing some  values. In edit mode, when I click on any one of these  values it gets replaced by 1.#QNAN in the cell.

I would like to query the table for all cells containing 1.#QNAN values. When I try to query for Null values (SELECT FROM .... WHERE .... IS NULL) no rows get selected. So Im assuming SQL doesn't recognize that 1.#QNAN are in fact null values.

How do I get the select by attributes dialog box to select all cells not containing a number?


texts:
I am trying to enable editing of a GeoJSON layer in my Leaflet map using Leaflet Draw. I need to allow only a single edit at a time before saving the results back to the database. I tried to implement the answer here (https://stackoverflow.com/questions/26932125/how-to-allow-only-one-feature-polygon-to-be-edited-at-a-time-with-leaflet). However, each time I activate the draw tool to edit a feature I get the following error in my console:

 Uncaught TypeError: Cannot set property 'previousOptions' of undefined

These are the key bits of code

My drawnItems layer:

drawnItems = new L.FeatureGroup();
map.addLayer(drawnItems)

My getJSON:

$.getJSON(url, function(data){
            geojsonLayer = L.geoJson(data,{
                onEachFeature: function (feature, layer) {
                    drawnItems.addLayer(layer);
                }
            });
            map.addLayer(drawnItems);
        });

My draw control:

var drawControl = new L.Control.Draw({
    draw: false,
    edit: {
        featureGroup: drawnItems,
        remove: false,
        edit: true
   }
});

map.addControl(drawControl);

map.on('draw:edited', function (e) {
        var layers = e.layers;
        layers.eachLayer(function (layer) {
            console.log(layer)
        });
    });

My GeoJSON is definitely valid (checked with geojson.lint) and have tried to just add a single feature rather than a FeatureCollection.

I have tried to follow some other posts (https://github.com/Leaflet/Leaflet.draw/issues/187) but none have solved my problem.


texts: Sequential attribute field calculations with null values using Arcpy

texts:
I have a series of x,y co-ordinates converted from a .ply file. I was able to extract the file into a suitable x,y format for arcmap, the issue I have is that they are co-ordinate points representing 90 polygons. So in arcmap I have the coordinates for polygon 1, then a series of null rows, then co-ordinates for poly 2 etc. I have been manually entering the polygon number using field calculator but with 91 it is taking me a long time, so I thought I could use the sequential calculations script. Each null row is the same number of rows between each poly (7). so what I need is a script that will update a field with polygon 1, skips 7 null field rows, poly 2, etc etc.

The other issue is that the number of points in each polygon is different, so the only way to differentiate each polygon is the null rows in between.

I meant this script in particular

http://help.arcgis.com/EN/ARCGISDESKTOP/10.0/HELP/index.html#//005s0000002m000000

Accumulative and sequential calculations


texts: proj4.js convert projection from 82397 to 4326

texts:
how to add utm84 51n to proj4
So i see something like this

proj4.defs("mySphere", "+proj=longlat +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=12345 +b=12345 +units=m +no_defs");

i want to add

proj4.defs("EPSG:36251", "????");

proj4('EPSG:36251','EPSG:4326',[long,lat]);

i don't understand that +proj.

EDIT
I already found the right projection

   proj4.defs("EPSG:32651", "+proj=utm +zone=51 +datum=WGS84 +units=m +no_defs");


texts: Calculating fields using a list of features

texts:
I've written some basic code which creates a new .dbf and adds two fields to it: 1 for today's date, another which should contain the names of shapefiles located in another folder. What I'd like to do is then write the names of the shapefiles to the newly created .dbf.

env.workspace = r'shapefile folder'

AT = r'location of new .dbf'
appendField = "ProjName"
dateField = "DateProc"
arcpy.AddField_management(AT, appendField, "TEXT")
arcpy.AddField_management(AT, dateField, "DATE")

shpList = arcpy.ListFeatureClasses()
for shp in shpList:
    print shp
    appendExpression = '"'+str(shp).replace('.shp','')+ '"'

    arcpy.CalculateField_management(AT, "ProjName",appendExpression, "PYTHON")
    arcpy.CalculateField_management(AT, dateField, 'time.strftime("%d/%m/%Y)', "PYTHON" )

I believe my issue is a result of having the table and feature classes in two separate workspaces. I am using 10.1. Thanks for any tips which may get this working.


texts:
What is the algorithm that Shapely used to check if two polygons intersect?
from shapely.geometry import Polygon
p1 = Polygon([(0,0), (3,0), (3,1), (1,1), (1,2), (3,2), (3,3), (0,3)])
p2 = Polygon([(4,0), (5,0), (5,1.5), (2,1.5), (2,1.2), (4,1.2)])
print(p1.intersects(p2))

I had read the source code of Shapely but didn't find the implementation.


texts: Recently I have been getting errors in FME Translations from a geodatabase (Reader) to a GT Viewer Format (Writer).  When I update or import feature types I get this error  I have done updates before with same client data from previous months. Not sure what is going on?


texts:
I am trying to process SENTINEL-1 GRD data using python and gdal. I want to get the shapefile from the input data.
I have tried to use gdal.Open() using the following extensions but nothing seems to work:

S1A_IW_GRDH_1SDV_20150109T163157_20150109T163222_004097_004F52_6A00.SAFE/manifest.safe
S1A_IW_GRDH_1SDV_20150109T163157_20150109T163222_004097_004F52_6A00.SAFE
S1A_IW_GRDH_1SDV_20150109T163157_20150109T163222_004097_004F52_6A00.tif
S1A_IW_GRDH_1SDV_20150109T163157_20150109T163222_004097_004F52_6A00.zip

gdal.UseExceptions()
srs = osr.SpatialReference()

src_ds = gdal.Open( "S1A_IW_GRDH_1SDV_20150109T163157_20150109T163222_004097_004F52_6A00.SAFE" )
if src_ds is None:
    print 'Cannot open %s' % src_filename
    sys.exit(1)

try:
    srcband = src_ds.GetRasterBand(1)
except RuntimeError, e:
    print 'Band ( %i ) not found' % band_num
    print e
    sys.exit(1)

dst_layername = "output"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

Does anyone know how to handle these data using gdal?

Best,
Elisavet


texts: Behaviour of Select Feature by Freehand tool in my own plugin

texts: I am calling QGIS "Select Feature by Freehand tool" in my plugin. what it will does is it allows the user to select features in mapcanvas and i made a signal like if selection changed in mapcanvas , will display the number of features selected in mapcanvas in a information box.This is what my plugin does. And what i am facing now is even without using the plugin if i select the features using select freehand tool, it does the same job as what my plugin does.Why does this behaviour in QGIS? What is my solution for this?should i not trigger the Select feature by freehand tool in my own plugin?


texts: Convert Landcover Raster to Vector to Use for Tracking NDVI Over Time

texts:
I want to look at NDVI and other spectral indices over time and compare these time series across different land cover types. I believe I have the code to do this but I'm running into the error that I have too many elements: Error generating chart: Collection query aborted after accumulating over 5000 elements. Would combining all those vectors into a single multi-part vector help? I would really like to do this over a much larger area so I assume I'll need a different approach?

// Function to cloud mask from the pixel_qa band of Landsat 8 SR data.
function maskL8sr(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = 1 << 3;
  var cloudsBitMask = 1 << 5;

  // Get the pixel QA band.
  var qa = image.select('pixel_qa');

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
      .and(qa.bitwiseAnd(cloudsBitMask).eq(0));

  // Return the masked image, scaled to reflectance, without the QA bands.
  return image.updateMask(mask).divide(10000)
      .select("B[0-9]*")
      .copyProperties(image, ["system:time_start"]);
}

// fcn to add NDVI band to each image
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

// Use the NLCD to filter areas by cover type and see how NDVI changed in those areas
// pull in NLCD
var NLCDdataset = ee.ImageCollection('USGS/NLCD');
var landcover = NLCDdataset.select('landcover');
var landcoverVis = {
  min: 11.0,
  max: 95.0,
  palette: [
    '5475A8', 'ffffff', 'E8D1D1', 'E29E8C', 'ff0000', 'B50000', 'D2CDC0',
    '85C77E', '38814E', 'D4E7B0', 'AF963C', 'DCCA8F', 'FDE9AA', 'D1D182',
    'A3CC51', '82BA9E', 'FBF65D', 'CA9146', 'C8E6F8', '64B3D5'
  ],
};
Map.addLayer(landcover, landcoverVis, 'Landcover');

// mask to your land cover class of interest
var landcover2011 = landcover.filter(ee.Filter.eq('system:index', 'NLCD2011')).first()

// so first we need a vector with all the wetland area

// Create a raster that is just the wetland class.
var wetland_raster = landcover2011.eq(95);
// mask the wetland raster to areas that are only wetland
wetland_raster = wetland_raster.updateMask(wetland_raster.neq(0));

// Convert the zones of the wetland areas to vectors.
var wetland_vectors = wetland_raster.addBands(landcover2011).reduceToVectors({
  geometry: roi, // this doesn't need to be just the ROI, could be much larger
  crs: wetland_raster.projection(),
  scale: 30,
  geometryType: 'polygon',
  eightConnected: false,
  labelProperty: 'zone',
  reducer: ee.Reducer.mean()
});

print(wetland_vectors);
//print(ee.Geometry.MultiPolygon(wetland_vectors))

// Make a display image for the vectors, add it to the map.
var display = ee.Image(0).updateMask(0).paint(wetland_vectors, '000000', 3);
Map.addLayer(display, {palette: '000000'}, 'vectors');

// now let's look at NDVI at this section over time

// Create image collection of l8 imagery
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
//filter start and end date
.filterDate('2019-03-01', '2019-04-01')
//filter according to drawn boundary
.filterBounds(roi)
.map(maskL8sr);

// Add NDVI band to image collection
var l8 = l8.map(addNDVI);

// Create a time series chart.
var plotNDVI = ui.Chart.image.seriesByRegion(l8, wetland_vectors, ee.Reducer.mean(),
'NDVI',30,'system:time_start', 'system:index')
              .setChartType('LineChart').setOptions({
                title: 'NDVI time series',
                hAxis: {title: 'Date'},
                vAxis: {title: 'NDVI'}
});

print(plotNDVI);

Here is a link to this code: https://code.earthengine.google.com/64daced9e1c453f1fcd98180196e410d


texts: Join data from ZonalStatistics output to original zone shapefile?

texts:
I am currently using arcpy to compute zonal statistics (mean) for a large number of polygons and a raster file using the following function

in_polys = some_shapefile
in_raster = some_raster

arcpy.ZonalStatistics(in_polys, "some_field", in_raster, statistics_type = "MEAN")

This produces the desired raster with the appropriate values. However, I wish to "join" the average values in the raster to the original zone shapefile to this effect:

ID Year  MEAN
1  2001  3.2
2  2002  1.2
3  2003  7.2

Is there a straightforward way to do this in arcpy?


texts: How customize GeoRSS icon in ArcGis JS API?

texts:
I have this little piece of script to customize the icon of some GeoRSS points, but I have no display and no error in the console.

  var map;
    require(["esri/map", "esri/layers/ArcGISTiledMapServiceLayer",
            "esri/layers/GeoRSSLayer", "esri/InfoTemplate", "dojo/parser",
            "dojo/_base/array", "dojo/dom-style", "dijit/layout/BorderContainer",
            "dijit/layout/ContentPane", "esri/symbols/PictureMarkerSymbol",
            "esri/renderers/SimpleRenderer",
            "dojo/_base/Color","dojo/domReady!"],
    function(Map, Tiled,
            GeoRSSLayer, InfoTemplate,
            parser, arrayUtils, PictureMarkerSymbol,
            SimpleRenderer, Color,  domStyle ) {

        map = new Map("map");
        var carte = new Tiled("http://services.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer");
        map.addLayer(carte);
        var georss = new GeoRSSLayer("http://dl.dropbox.com/u/2654618/data-mixed.xml");
           georss.on("load", function() {
                    domStyle.set("loading", "display", "none");
        // create an info template
                    var template = new InfoTemplate("${name}", "${description}");
  // set the info template for the feature layers that make up the GeoRSS layer
  // the GeoRSS layer contains one feature layer for each geometry type
                    var layers = georss.getFeatureLayers();
                    var picRenderer = new PictureMarkerSymbol("http://static.arcgis.com/images/Symbols/Shapes/YellowCircleLargeB.png", 16, 16);
                    var georssRenderer = new SimpleRenderer(picRenderer);
                    arrayUtils.forEach(layers, function(l) {
                    l.setInfoTemplate(template);
                    l.setRenderer(georssRenderer)
});
});
                    map.addLayer(georss);
                                        //THE END
});


texts: Handle Event When User Selects A Feature

texts:
How do I handle the event when a user selects a feature from a feature layer in ArcObjects?

What I've tried:

    public ESRI.ArcGIS.Carto.IActiveViewEvents_SelectionChangedEventHandler m_ActiveViewEventsSelectionChanged;

    public ArcObjectsController()
    {
        ESRI.ArcGIS.Carto.IActiveViewEvents_Event activeViewEvents = (ESRI.ArcGIS.Carto.IActiveViewEvents_Event)ArcMap.Document.FocusMap;
        m_ActiveViewEventsSelectionChanged = new ESRI.ArcGIS.Carto.IActiveViewEvents_SelectionChangedEventHandler(OnActiveViewEventsSelectionChanged);
        activeViewEvents.SelectionChanged += m_ActiveViewEventsSelectionChanged;
    }

    public void OnActiveViewEventsSelectionChanged()
    {
        MessageBox.Show("Event Fired");
    }


texts: QGIS 3.10: attribute form showing foreign key instead of referencing field

texts:
In a QGIS 3.10 project, I've set up a relationship between two postgis table. In postgis, table A has two unique fields: "id" and "name". "Id" is the primary key but "name" is used as a foreign key for table B as "name A".
In the QGIS project's properties, I've set up the relation as so:
Parent layer: layer A
Referenced field: Name
Child Layer: layer B
Referencing field: Name A
This works as expected, I can then create attribute form for both layers, with the form of layer A having a section for showing layer B's attribute form.
However, I noticed that when I open a the attribute table of layer B, the field which should be showing "name A" is instead showing the "Id" field of layer A.
How can I set up the layer B form so that it shows "name A" as expected? (note I do not have this"error" in Postgis when i check the database in PgAdmin). Maybe I should set up a "reverse" relationship in QGIS?


texts: Converting a Road shp file into a continous Road Raster in QGIS?

texts: I have converted a Road shp file into a Raster in QGIS, however after the conversion the Road is no longer continous (i.e only the cell corners are touching). I have tried increasing the resolution during conversion however that has not solved the problem. The first image is what I need, the 2nd is my current result.


texts:
When geocoding a dataset of about 500 addresses I get a vector point layer. So far everything is fine. The problem is, that the dataset contains multiple entries with the same address. The result is a point stack, in which one point covers one or more points that lie beneath the point on the top. Of course this is a correct behaviour, because: same address = same coordinates. But: When viewing the map, information is lost.

So my question: Is there a possibility to automatically arrange these points not on top of each other but side by side?


texts: An issue with the ee.clip of Google Earth Engine Python API

texts:
I'm new to Google Earth Engine analysis so I can't solve a simple task with the python API.
What I want to do is to clip an image with a specified boundary that corresponds to a reference image then export it to my google drive. When I execute the code, the image is exported, and its appearance is not bad. However, the upper boundary of it is different which I have specified. I've been looking and testing for a solution to this problem for a quite long time, but still, I can't find a good solution for it. The code and the information regarding my problems are attached below. The coordinates of the boundary are acquired by using rasterio.bounds. The target ROI is the coordinates of the reference image.
Please help me!
Product: MODIS MOD14A1
Target ROI: left=97.99955519999997, bottom=50.000064, right=117.99952319999997, top=58.0000512
The boundary of the output file: left=97.9972143445986, bottom=49.991245561251375, right=118.00269572194034, top=58.399476620610095
#The band that I want to use
band=['MaxFRP']

#Set a region of interest
ROI=ee.Geometry.Rectangle(coords=[97.99955519999997,50.000064,117.99952319999997,58.0000512])

#Get a collection of images
output=ee.ImageCollection("MODIS/006/MOD14A1").filterDate('2001-04-01','2001-04-05').select(band)

#Extract the first image of the image collection
output_first=output.toList(output.size().getInfo()).get(0)

#Clip the image by ROI
output_first_image_clipped=ee.Image(output_first).clip(ROI)

#Export the image to google drive
task=ee.batch.Export.image.toDrive(
    description='testImage',
    image=output_first_image_clipped,
    scale=1000,
    crs='EPSG:4326'
)
task.start()

I also attached the reference image over the output image from the GEE. The black area shows the unwanted area due to the wrong top-line coordinate.


texts: Is there an option to remove duplicate polygon labels in QGIS?

texts:
I would like to know if there is an option in QGIS labeling engine which would remove all the duplicate labels.

My polygon layer contains a few features with same name. So when labeled it becomes a cluster which does not look good.


texts: Error generating chart: Remote Request too large in Google Earth Engine

texts:
I was trying to generate to spectral chart for 'VV','VH','angle' using Google Earth Engine. But I am getting the error

  Error generating chart: Remote Request too large

I was using the following code:

    var test1 = data
    .filterDate('2017-06-16', '2017-08-30')
    .select('VV','VH','angle')
    .filterBounds(boundary);

var stack = function(collection)
{
  var bond = boundary
  var first = ee.Image(collection.first()).select([]);
  var appendBands = function(image, previous)
  {

    return ee.Image(previous).addBands(image);
  };
  return ee.Image(collection.iterate(appendBands, first));
};

var stacked = stack(test1)
var mosaic = ee.ImageCollection(stacked).mosaic()
var clip = mosaic.clip(boundary)
print (stacked,'stacked')
print (clip,'clip')
var bandNames = stacked.bandNames();
print (bandNames)

Map.addLayer(clip)
print (clip,'clip')

var chart = ui.Chart.image.series({
  imageCollection: clip,
  region: paddy,
  reducer: ee.Reducer.mean(),
  scale: 10
});

chart.style().set({
  position: 'bottom-right',
  width: '500px',
  height: '300px'
});
Map.add(chart);

What is wrong?

Earth Engine code link https://code.earthengine.google.co.in/82a90af6ef77a68d9e228b8aad76f17a


texts: Apply regression equation to Landsat ETM+ image

texts:
I have sample points in an Excel file with XY coordinates and canopy density. I also have a Landsat ETM+ image. All points are collected from the same area as in the image.

I exported all sample points by first ‘Display XY data’ and then ‘Export data’, thus I got a shapefile.
I used "Sample" extracted all points from ETM+ image, then I got a table and exported it to a DBF file. In the DBF file I have XY coordinates for all sample points and band values.
Then I copied all band values from DBF to original sample points Excel file, then I calculated NDVI and did linear regression (canopy density as dependent variable and NDVI as independent variable) in Excel. Therefore I got a regression equation.
After added the Excel to ArcGIS and joined it to sample points shapefile, I carried out OLS. However, I only got all sample points regressed.

What could be the problem and how should I apply regression equation to the image so I can produce a map?


texts: Modifying the LocateWidget.mxml Locate and Clear functions

texts:
I am using Robert's compiled TOC widget and attempting to modify the LocateWidget.mxml for the flexviewer. The workflow I have in mind is that user performs an address search clicks the locate button and all my defined operational layers turn on. The clear button would then turn them off.

My attempt at the code is posted here:
http://forums.arcgis.com/threads/57851-Modifying-the-LocateWidget.mxml-Locate-and-Clear-functions


texts: How to extract area information from a raster based on a DEM or contour lines

texts:
I have a raster representing the distribution of a specific land cover type - the raster values are 1 for where the land cover type is present and 0 where it isn't.

I also have a DEM with 5m resolution from which I extracted contour lines, and I want to know how much of that land cover type (in terms of area) falls inside each 200m elevation band (0-200, 200-400, etc.).

Can anyone suggest the best way to do this in QGis? I thought about Zonal statistics but I need polygons for that and not lines. Or is it better to use the DEM directly? I'm not sure if that would be a good option because the rasters have different resolutions (2 m vs 5 m) and their grids aren't aligned.


texts: How to populate M values for an existing x,y,m route layer in ArcGIS 10?

texts:
Very surprised to not find this in docs, I swear this was easy in former versions... Whatever, here we go:

Say you have a rivers layer of x,y,m polylines. Each feature is a named river. Each feature was created in the direction of flow, so that the first vertex of each feature is the river's headwater, and the last is its drain (probably where it flows into the next river). I want to assign measures to each vertex of a given feature as the accumulated length of the feature from its drain to that point.

I've only found the 'CreateRoutes' ArcToolbox tool, which not only does the unnecessary work to aggregate features and output a new layer, but it has a wonky idea of 'coordinate priority' to define each feature's measures start point, being that closest to a corner of its bounding box... I want to operate on my layer in-place, and I simply want to start the measures from the start (or end, in my case) of each feature.

One possibility with that tool might be to use its 'TWO_FIELDS' Measure Source, where I use a pair of columns calc'ed with zeros (first column) and feature length (second column)...

I thought I'd check first if anyone else has already been down this path?


texts: How to convert legend to graphics in QGIS?

texts:
In QGIS am trying to convert the legend element in the print composer into graphics, like the option in the ArcGIS "Convert To Graphics", see image below

Is there some way to do it with QGIS?


texts: How to measure the compactness of a set of points

texts:
I'm looking for ideas on how to measure the compactness of a set of points. The idea is similar to identifying clusters, but instead I want to compare two distinct sets of points and return a measure of how compact one set is compared to the other.

So far the only research I've come across is specific to finding clusters within a set, rather than measuring the whole set.

Can anyone point me to useful algorithms for accomplishing this? In particular, it would be great to find an existing tool in QGIS that does this, although that's not necessary.


texts: Making automatic updated date and time field in QGIS

texts:
I wish to create a field which automatically updates with the current timestamp when edits are made to a feature.

I am using QGIS 3.4 and I am testing this using a shapefile.

Following @Joseph's answer for this question I tried updating his script to work in QGIS 3.4 as per his instructions (if I followed them correctly, not guaranteed!) to this:

layer = qgis.utils.iface.activeLayer()

def update():
    field = layer.fields().lookupField('fieldName')
    e = QgsExpression( " $now " )
    e.prepare( layer.pendingFields() )
    for feat in layer.selectedFeatures():
        feat[field] = e.evaluate( feat )
        layer.updateFeature( feat )

layer.attributeValueChanged.connect(update)

However this returns the following error:

2018-11-29T11:13:00     WARNING    Traceback (most recent call last):
              File "", line 4, in update
             AttributeError: 'QgsVectorLayer' object has no attribute 'pendingFields'

What's wrong here?

Also will this work on geometry edits or only with attribute edits?

Will this also work for a SpatiaLite layer?


texts:
I have a PostgreSQL table of students with lat, lon and geom, and a table of schools with the same, and a table of roads that has been converted to nodes.

I'm trying to find the driving distance to the nearest school for each student. But I'm new to PostgreSQL and pgRouting, but I'm not seeing a great way to structure the multiple joins to work with the pgr_drivingDistance function.

Ideally I would end up with a new table that would have the student id, the id of the closest school via driving distance, and the distance.

Any thoughts on how to approach this would be great.

Update addressing Jendrusk's questions:

What do you mean by converting roads to nodes? Maybe you meant network graph?

select pgr_createTopology('roads', 0.0001, 'geom', 'id')

Are you able to find at least one way and the trouble is with many ways?

The challenge is finding the shortest way for multiple students. I could take the lat/lon for one student, find the nearest node and then calculate distances that way, but I need to do it for hundreds.

Nearest school euclidean or by the roads?

Driving distance.

What was the work-flow of this conversion form roads to nodes?

See above.

What are the names of tables and column names in this tables?

Three tables: Kids, roads, schools. Kids and schools have a GID, lat, lon, geom, and a bunch of other non-spatial info. Roads has GID, geom, source, target and length.

How big are your tables with students and schools?

100k students, several dozen schools.

Is it OK to select nearest school euclidean

No, interested in driving distance.

      Also function st_DrivingDistance is not finding route, but it's creating set of edges reachable in specified time (e.g. how far you can go in 15 minutes).

    Also 2 you need distance from student to nearest school but which is nearest if you don't know the distance? Do you have some kind of relation between schools and students? If not you have to count distance from one student to every school to know which is nearest. In most cases nearest euclidean (in straight line) will be also nearest by the roads, but if you have some special cases in your graph such as river with one bridge you can't be sure.<<<

Again, I do not need to find the route. I need to identify the closest school based on driving distance. My impression is that using pgr_drivingDistance it is possible to create a table containing driving distance length from a specified node to all other nodes in a network. You can scope this operation to nodes within a particular distance of the index node, and obviously, in my case, only consider nodes that are the closest node to a school.

I'm basing this info on the following post: http://anitagraser.com/2011/05/13/catchment-areas-with-pgrouting-driving_distance/

I believe pgr_drivingDistance can limit the query based on a number of different costs, including driving time and distance, depending on how the cost is specified both in the query and in the data.


texts: How to make a slider iterate in Google Earth Engine?

texts:
I have a bunch of layers displayed, and I can make a slider so I can click through the layers (see code). So how do I make the slider bar display from 0 to NumLayers and iterate by 1 layer? (Right now I see a slider from 0 to 1 and I have to slide the bar very carefully to end up on a layer of interest).
var slider = ui.Slider();

slider.onSlide(function(value) {

  var int_value = value * (Map.layers().length() - 1) >> 0;

  Map.layers().get(int_value).setOpacity(1);

  for (var i = int_value + 1; i < Map.layers().length(); i++) {

    Map.layers().get(i).setOpacity(0);

}

});

print(slider);


texts: Ways to Speed Up Python Scripts Running As ArcGIS Tools

texts:
This is a pretty general question. I just wondering what tips and tricks GIS programmers have used to speed up arcpy scripts that you import into the toolbox and run.

I work most everyday writing little scripts to help non-GIS users at my office process GIS data. I have found that ArcGIS 10.0 processing in general is slower than 9.3.1 and sometimes it gets even slower when running a python script.

I'm going to list a particular example of a script that takes over 24 hours to run. It's a loop that tabulates the area of a raster in a buffer for each shape in the buffer. The buffer has about 7000 shapes. I don't believe it should run this long. A

while x <= layerRecords:

    arcpy.SetProgressorLabel("Tabulating Row: " + str(x) + " of " + str(ELClayerRecords))
    arcpy.SelectLayerByAttribute_management(Buff,"NEW_SELECTION", "Recno = " + str(x))                                  # Selecting the record
    TabulateArea(Buff, "Recno", MatGRID, "VALUE", ScratchWS + "/tab" + str(z) +".dbf", nMatGRIDc)                          # Tabulate the area of the single row

    arcpy.AddMessage ("          - Row: " + str(x) + " completed")
    x = x + 1
    z = z + 1

Before anyone says it, I have run tabulate area on the entire buffer, but it produces errors if run on more then 1 record. It's a flawed tool, but I have to use it.

Anyways, if anyone has any ideas on how to optimize, or speed up this script, it would be most appreciated. Otherwise, di you have any speed up tricks for python, when used in ArcGIS?


texts:
I have a layer which has one feature. I try to move the feature(multipolygon) using the script below:

layer = iface.activeLayer()

for feat in layer.getFeatures():
    prev_pos = feat.geometry().centroid().asPoint()
    success  = feat.geometry().translate(100, 100) # translation
    next_pos = feat.geometry().centroid().asPoint()

    print(success)
    print(prev_pos)
    print(next_pos)

But I get different results in QGIS 2 and QGIS 3. (The multipolygon has many points, therefore, I printed out only its centroid)

Output in QGIS 2:

0 # means operation succeeded
(552084,433685)
(552184,433785) # changed coordinates

Output in QGIS 3:

0 # means operation succeeded
<QgsPointXY: POINT(552084 433685)>
<QgsPointXY: POINT(552084 433685)> # same coordinates

When I try to move a point feature, I get same results.

Why does not translate method move the feature in QGIS 3? What am I missing?

EDIT: If I create a QgsGeometry  object manually, translate method works correctly in QGIS 3.

points = [QgsPointXY(0,0), QgsPointXY(2,0), QgsPointXY(2,2), QgsPointXY(0,2)]
geom = QgsGeometry.fromMultiPolygonXY( [ [points] ] )
print(geom)

geom.translate(100, 100)
print(geom)

Output in QGIS 3:

<QgsGeometry: MultiPolygon (((0 0, 2 0, 2 2, 0 2, 0 0)))>
<QgsGeometry: MultiPolygon (((100 100, 102 100, 102 102, 100 102, 100 100)))>


texts:
My map is in a HTML page
Lots of markers will be placed on the map
I need to save this map with all the marker in a img (.png or .jpg/.jpeg preferably) format, not the entire HTML page

Is there a plug-in from leaflet which allows us to do such a thing ?


texts:
We are developing a QGIS plugin and would like to override the QGIS Feature Edit Form with our own Dialog.

Have explored using setEditForm and setEditFormInit, but the init process appears to want to use a python file located in the user's project directory.

Ideally we'd use an init method from our Plugin's code and gain access to our Plugin's instance.

Had thought that the new EditWidget V2 might be able to help, but this seems to be focussed on editing individual attributes rather than a feature.

UPDATE

Doh. Managed to resolve it.

In the init method of your plugin:

Firstly add the plugin directory to the system path list (or a subdirectory of it)
Secondly, initialise a static variable with the Plugin's instance (self)

Am yet to add logic to the dialog but looks like we're on the home run.  Thanks Python!


texts: Transforming raster in range of values from 3 to 8 into linear scale from 0 to 1

texts:
I am trying to convert a raster file with values ranging from 3.2 to 8.8 into a linear scale from 0 to 1. The conditions that i want to apply in the transformation are as follows:
3 to 5 gets 0,
5 to 7 Linear variation from o to 1 and gets values based on straight line equation y = (slope)x + c,
And greater than 7 gets 1
How is it possible to create linear transformation as a tool in ArcGIS using python code where i can input those conditions mentioned above as  parameters in the tool?
I would like to distribute this tool to others and they may have different control points instead of 3,5,7, and greater than 7.
import sys
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
env.workspace = r"H:\Test"
arcpy.env.cellSize = 100
Input = arcpy.GetParameterAsText(0)
folder = arcpy.GetParameterAsText(1)
raster = arcpy.GetParameterAsText(2)
arcpy.AddMessage("Folder: {}".format(folder))
arcpy.AddMessage("Raster: {}".format(raster))
save_raster = os.path.join(folder, raster)
arcpy.AddMessage("Save Location: {}".format(save_raster))
output = RescaleByFunction("Elevation", TfLinear(3, 5, "#", "#", "#", "#"), 0, 1)
output.save(save_raster)```


texts: Finding coordinates from nested grid as numpy array?

texts:
I've written a script that allows a user to read in a raster and provide a cell resolution, and outputs the amount of cells necessary to resample the grid (to provide a pre-processing tool to run the WRF weather model). I am still trying to figure out a test to make sure the cell size actually fits the input grid.

I'm having some difficulty figuring out how to nest grids within a 'master' domain.  IDEALLY this would be via user-IO, where multiple grids are selected within themselves. Haven't added this function yet, but the nesting is easy enough with a for loop (though I need to make sure that the nested grids actually fit within the master)

How do I obtain the origin points for each of the nested grids in terms of cells of the master domain (shown in the link above)?

For instance, finding the x,y coordinates of the 0,0 origin of a nested domain WITHIN the larger master-domain.

Here's my code:

#!/bin/python
import numpy as np
from sys import argv

file = argv[1]
grid = np.loadtxt(file)
array_dims = grid.shape
x = array_dims[0]
y = array_dims[1]

cell_size = raw_input("What is the cell size you want? (x,y) \n >")
cell_size = (cell_size.split(","))
cell_size = [int(i) for i in cell_size] #goes through list and turns them to array
print cell_size
cell_x = cell_size[0]
cell_y = cell_size[0]

print "The shape of the grid is %r by %r \n" %(x,y)

e_we = int(round(x/cell_x))
e_sn = int(round(y/cell_y))

print "With that cell size, the WRF-domain is %r by %r cells" %(e_we,e_sn)

center_x = x/2
center_y = y/2

print "The center point is at %r, %r" %(center_x,center_y)
print "This point's value is", grid[center_x,center_y]


texts:
I want to clear the features in a vector layer source , I have tried

features.getFeatures().clear();

but it doesn't work!


texts: Is there any STAC catalog available for Sentinel-2 data on Google Cloud?

texts:
We are looking for a fast and optimum way of discovering and downloading Sentinel-2 data from Google Cloud using Python.
We have tried using index data in BigQuery and it works. However, we have also come across Python packages (such as: sat-search https://github.com/sat-utils/sat-search) for searching and downloading publicly available satellite imagery using STAC compliant APIs that we would prefer to use, but we haven't found any STAC catalogs for Sentinel-2 datasets on Google Cloud or any STAC API endpoints for the same. Is anyone aware of any STAC catalogs for datasets on Google Cloud or if there are plans to have them in the near future?
AWS has both: Sentinel-2-L1C and L2A datasets, and has recently released these in COG (Cloud Optimized Geotiff) format which seems to be better suited to large-scale and distributed computing. Does anyone know if similar datasets will be available on Google Cloud as well?


texts: On the fly/Dynamic measuring not working on Measurement Widget with different Basemap

texts:
I am developing an application with OS Open Carto Map as a Basemap and I have added a Measurement widget in the app.
The measurement widget works fine for when the points are clicked for distance and area, but it doesn't calculate the distance on the fly/dynamically (before you set/click the second point) as in this ESRI Demo.

Here is the code:

<script>
var map;
require([
    "dojo/dom",
    "esri/basemaps",
    "esri/Color",
    "dojo/keys",
    "dojo/parser",
    "esri/config",
    "esri/sniff",
    "esri/map",
    "esri/SnappingManager",
    "esri/dijit/Measurement",
    "esri/layers/FeatureLayer",
    "esri/renderers/SimpleRenderer",
    "esri/tasks/GeometryService",
    "esri/symbols/SimpleLineSymbol",
    "esri/symbols/SimpleFillSymbol",
    "esri/dijit/BasemapLayer", "esri/dijit/Basemap",
    "esri/dijit/Scalebar",
    "dijit/layout/BorderContainer",
    "dijit/layout/ContentPane",
    "dijit/TitlePane",
    "dijit/form/CheckBox",
    "dojo/domReady!"
  ], function(
    dom, esriBasemaps, Color, keys, parser,
    esriConfig, has, Map, SnappingManager, Measurement, FeatureLayer, SimpleRenderer, GeometryService, SimpleLineSymbol, SimpleFillSymbol,
    BasemapLayer, Basemap
  ) {
    parser.parse();

    esriConfig.defaults.io.proxyUrl = "/proxy/";
    esriConfig.defaults.io.alwaysUseProxy = false;

    //This service is for development and testing purposes only. We recommend that you create your own geometry service for use within your applications
    esriConfig.defaults.geometryService = new GeometryService("https://utility.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer");

    var bounds = new esri.geometry.Extent({
        "xmin":1393.0196,
        "ymin":13494.9764,
        "xmax":671196.3657,
        "ymax":1230275.0454,
        "spatialReference":{"wkid":27700}
    });

    var lods = [
        {"level" : 0, "scale" : 5000000, "resolution" : 1322.9193125052918},
        {"level" : 1, "scale" : 2500000, "resolution" : 661.4596562526459},
        {"level" : 2, "scale" : 1000000, "resolution" : 264.5838625010584},
        {"level" : 3, "scale" :  500000, "resolution" : 132.2919312505292},
        {"level" : 4, "scale" :  250000, "resolution" : 66.1459656252646},
        {"level" : 5, "scale" :  100000, "resolution" : 26.458386250105836},
        {"level" : 6, "scale" :   50000, "resolution" : 13.229193125052918},
        {"level" : 7, "scale" :   25000, "resolution" : 6.614596562526459},
        {"level" : 8, "scale" :   10000, "resolution" : 2.6458386250105836},
        {"level" : 9, "scale" :    5000, "resolution" : 1.3229193125052918},
        {"level" :10, "scale" :    2500, "resolution" : 0.6614596562526459},
        {"level" :11, "scale" :    1000, "resolution" : 0.26458386250105836},
        {"level" :12, "scale" :     500, "resolution" : 0.13229193125052918},
        {"level" :13, "scale" :     250, "resolution" : 0.066145965625264591},
        {"level" :14, "scale" :     100, "resolution" : 0.026458386250105836},
    ];

    var myLayer = new BasemapLayer({
       url: "https://tiles.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Carto_2/MapServer"
    });
    var myBasemap = new Basemap({
       layers: [myLayer],
       title: "OS Open Carto",
        //thumbnailUrl: "https://www.arcgis.com/sharing/rest/content/items/86de95d4e0244cba80f0fa2c9403a7b2/info/thumbnail/tempimagery.jpg"
    });

    map = new Map("map", {
        extent: bounds,
        basemap: myBasemap,
        zoom: 1,
        lods: lods,
        showLabels : true
    });

    var div = document.createElement("div");
        div.id ="measurementDiv";
        document.getElementById("measureDiv").appendChild(div);

    var measurement = new Measurement({
      map: map
    }, dom.byId("measurementDiv"));
    measurement.startup();

  });
</script>


texts: NASA SRTM: precipitous contour lines in TileMill

texts:
I was trying to make a nice looking print map of Bali (Indonesia) in TileMill. After creating a terrain layer (using this walkthrough) I noticed that the island's coastline was precipitous.

Is this caused by GDAL or wrong input data from CIGAR? I believe it can't be fixed within TileMill as the hillshading layer seems also affected.

Edit: With a bit CartoCSS the contour lines can easily be hidden. But not the slope scale image...

TileMill 0.10.1 (v0.10.1-134-g31027ed, Windows)
GDAL: OSGeo4W on Win 32bit


texts: Cutting (not clipping) polygons with other polygon layer in QGIS

texts:
I'm working with QGIS 3.16 Hannover under MacOS 10.13 environment. I've got two polygon shapefile layers: one for natural areas and one for regions. I need to divide all natural areas, using the boundaries of the regions, without loosing any part. For example, in the following image, polygon A represents a natural area in shapefile Natural_Areas and polygons B and C represent the boundaries of two regions in shapefile Regions

I want to obtain a layer containing polygons 1, 2 and 3. I've tried with the clipping tool, but it returns a single polygon corresponding to polygons 2 and 3 joined (or, more precisely, not cut).


texts: Creating a memory layer without the CRS dialog in PyQGIS 3

texts:
I am trying to create a memory layer:

tempLayer = QgsVectorLayer('LineString', 'my_temp_layer', 'memory')

But as soon as this line is excecuted (not when I add it to the map!) the CRS dialog pops up. How do I intercept that?

I know I can set the CRS, but that does not help, since the dialog pops up when the first line is executed...

tempLayer = QgsVectorLayer('LineString', 'trassen_tmp_bfp', 'memory')
tempLayer.setCrs(QgsCoordinateReferenceSystem(4326))


texts: Loading Leaflet raster layers asynchronously before adding them to the map?

texts:
I am working on a web application where users can view sequences of 5-10 GeoTiffs overlaid on a map, and played like a slideshow - show the first Geotiff, wait a few seconds, show the second, and so on. Strictly one at the time. The Geotiffs are served using WMS by a Geoserver instance on the same server.

Things can get really slow fetching them (we're talking 4-8 seconds for tiffs in the 2-10Mb range), and users complained, so I thought of loading all the geotiffs in a slideshow before playing it. However, I'm under the impression that Leaflet does not request the raster from Geoserver until its layer is actually added to the map.

What is an efficient technique, if at all possible, to load a group of Geotiffs asynchronously before adding them to the map?


texts: Radiance Calculation of Landsat 7 with Google Earth Engine(search for gain and bias)

texts:
I need to calculate the Radiance of Landsat 7 in GEE
I read in Data Users Handbook under 5.6.5 how to do this. Here it says, these Values are accessible in the Metadata. Now I am searching the Metadata for Gain and Bias in GEE
Link to GEECode.
I made a list of all Property names and can not find it.
Pls help


texts: How do i popup a context menu with custom command items in a Select tool OnMouseDown (Right Click) with VB .Net?

texts:
This should be simple... I am converting VBA customization to VB .NET in VS 2008 Express Edition.
In VBA i have a custom selection tool. Aside from just selecting features, when OnMouseDown receives Button = 2 parameter (right click) I popup a context menu with some macro items that call various other procedures.

I need some help with making this work in .Net...

I used the ESRI .NET sample project to get the feel for this and i am sure i can re-create the select tool but i am having problems with the context menu. I would also like to know how to call procedures when user clicks the context menu item.

Private Sub UIToolControl1_MouseDown(ByVal button As Long, _
    ByVal shift As Long, ByVal x As Long, ByVal y As Long)

Dim pLayer As ILayer
Dim pContextMenu As ICommandBar
Dim pAppPosition As IWindowPosition
Dim pPoint As POINTAPI

If pMxDoc.SelectedLayer Is Nothing Then
MsgBox "Select a layer in the TOC!": GoTo ep
End If

Select Case button
Case 1
Set pLayer = pMxDoc.SelectedLayer
If pLayer Is Nothing Then
MsgBox "You must select 1 feature layer in the TOC!": GoTo ep
End If

If Not TypeOf pLayer Is IGeoFeatureLayer Then
MsgBox "You must select 1 feature layer in the TOC!": GoTo ep
End If

Set m_pPoint = _
    pActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(x, y)
m_bIsMouseDown = True
Case 2
'on right click = popup
Set pContextMenu = CommandBars.Create("Go2Grid", esriCmdBarTypeShortcutMenu)

Set pAppPosition = Application

' Add 3 built in commands to the new context menu using the
' built-in ArcID module.

pContextMenu.CreateMacroItem "New Feature Class", 8, _
    "Normal.Module1.CreateNewFeatureClass"
pContextMenu.CreateMacroItem "Delete Layer", 13, _
    "Normal.Module1.DeleteFeatureClass"
pContextMenu.CreateMacroItem "Copy To...", 20, _
    "Normal.ThisDocument.CopyFeatures"
pContextMenu.CreateMacroItem "Draws Around...", 26, _
    "Normal.ThisDocument.DrawsAroundFeatures"
pContextMenu.CreateMacroItem "Clear Tool", 10, _
    "Normal.ThisDocument.ClearCurrentTool"
pContextMenu.CreateMacroItem "Delete Features", 13, _
    "Normal.Module1.DeleteFeatures"

' Popup the menu.

pPoint.x = x
pPoint.y = y
ClientToScreen pActiveView.ScreenDisplay.hwnd, pPoint

pContextMenu.Popup pPoint.x, pPoint.y
Case Else
GoTo ep
End Select

End Sub


texts: Finding nearest tiepoint in aerial triangulation export from Bentley Context Capture in Cesium

texts:
I am trying to solve the problem of clicking a point on a 3D tileset loaded by Cesium, and finding the nearest tiepoint corresponding to the clicked point

Here are the steps that I have taken:

1) Get the Cartesian coordinates:

cartesian3d = viewer.scene.pickPosition(click.position)

2) Converting the Cartesian3 to WGS84 Lat/Long:

cartographic3d = cartesianToDegrees(viewer.scene.ellipsoid, cartesian3d)

const cartesianToDegrees = (ellipsoid, x) => {
  const carto = ellipsoid.cartesianToCartographic(x);
  const longitude = (carto.longitude * 180) / math.pi;
  const latitude = (carto.latitude * 180) / math.pi;
  return [longitude, latitude, carto.height];
};

I have a list of tiepoints from the aerial triangulation step of Bentley Context Capture, this is the reference system that was used: EPSG:28356

3) Now to convert the tiepoint:

tiepoint = (298886.689734013, 6229317.15626832, 26.6625023260713)

to WGS84 I use pyproj:

in_proj = Proj("EPSG:28356")

out_proj = Proj("+init=EPSG:4326")

wgs_coords = transform(in_proj, out_proj, x=tiepoint[0], y=tiepoint[1], z=tiepoint[2])

Which is:
(150.82094535413248, -34.05758120402254, 26.6625023260713)

I then do this for all the tiepoints, then I am able to query for the nearest tiepoint using a KD tree.

Now I click on my 3D Cesium tileset to get the cartographic3d coordinate, and try to find the nearest matching tiepoint. The problem is that retrieved doesn't match up to the point that I have clicked on the model, it always seems a few meters off.

I'm thinking that I must have done something wrong with the conversion but I can't seem to find out what it is. If I'm not wrong Cesium uses ECEF for Cartesian3d, and ellipsoid.cartesianToCartographic is called I can confirm we are using the WGS84 ellipsoid.

Am I missing anything?


texts: Using output statistics in ArcGIS REST API query?

texts:
I am using an API to pull in some Web Mapping Service Layer data to my script. The API query has a max count limit of 1000 on returned features per request. There are 2 attributes I am trying to use in my query 1) A DFIRM_ID field that has the same value for multiple features (its not unique but there are many DFIRM_ID values) and 2) The OBJECTID -this is a totally unique value. Each feature has its own OBJECTID. My goal is to get the total number of OBJECTID's for a given DFIRM_ID (or the total number of features). So, I want to use Output Statistics in my query. Which also requires 'Group By Fields'. According to the ArcGIS REST API documentation, I believe I should have ``. But the outputStatistics part seems confusing (as to how to put it in my query. The documentation says the following syntax (or something like it) is to be used:

[{
"statisticType": "<count | sum | min | max | avg | stddev | var>",
"onStatisticField": "DFIRM_ID",
"outStatisticFieldName": "DFIRM_COUNT"
}]

I am just a little confused. How would I use that in my params variable below when plugging it into my code?

params = {'f': 'json', 'where': 'DFIRM_ID like \'13245%\'+'and OBJECTID>='+str(minCount)+'and OBJECTID<='+str(maxCount), 'geometryType': 'esriGeometryPolygon', 'spatialRel': 'esriSpatialRelIntersects','outFields': '*', 'returnGeometry': 'false'}
r = requests.get('https://hazards.fema.gov/gis/nfhl/rest/services/CSLF/Prelim_CSLF/MapServer/3/query', params)
cslfJson = r.json()
master.append(cslfJson)


texts: Processing tool to fill all voids within a lines layer

texts:
I'm trying to convert a lines layer from Ordnance Survey which is made up of field boundaries, river banks, woodland boundaries and houses etc into a polygon layer, but only 'polygonise' the fields.

Below is the direction I've taken so far:

Ordnance Survey Data with Estate boundary in red:

I used an expression to select and save only the field boundaries, clipped it against the estate boundary and dissolved the resultant layer:

Highlighted fields have simple geometries and those not highlighted would need manual editing following the processing

Result when I do lines to polygons (same result prior to dissolve):

I tried intersecting the estate boundary and lines but this resulted in no features. I later found the SAGA processing tool 'Polyon-line intersection' which has given me the best result so far but still not great:

How can I  do this?


texts: Why can't I see my globe service in ArcMap?

texts:
I want to explore it in ArcMap or in ArcExplorer but both cannot "see" it. For the other resources like map services everything is working well.

I can add that the service is started and in a preview view everything seems to be ok.


texts: How to get rid of the transformation grid in st_transform?

texts:
I use the last versions of R -4.0.3- and sf -0.9.6- (thus relying on PROJ version 6.3.1).
I want to transform a layer from one system of coordinates (EPSG:27572) to another (EPSG:2154). This can be done with (accurate) or without (less accurate) using a transformation grid file. Once I have copied the appropriate transformation grid file (ntf_r93.gsb for PROJ 6) in the proj directory of sf (located at "C:\Program Files\R\R-4.0.3\library\sf\proj" on my PC), st_transform automatically use it and gives the most accurate results:
st_transform(MyLayer,crs=2154)
My question is:

how to get rid of the transformation grid, coming back to the non-gridded transformation, without removing the file in the proj directory?
Is there any sf command, or any parameter in st_transform to do this simply?
Perhaps using the pipeline parameter?
How?


texts: Not able to progress past loading CSV in New Dataset dialog

texts:
We're getting some strange behavior in the connect dataset dialog, where, once we've selected a CSV file for upload, we get no button or feedback to continue the process.

screenshot

Anyone else seeing this weird behavior, or have a fix? We're using Firefox, several versions, on OSX and Win.


texts: Where can I find a GDAL distribution with DWG support?

texts:
I want to use the GDAL utility ogrinfo with dwg files. Where can I find an GDAL distribution with a support for DWG files (read only)?

DWG support is moved out of the normal GDAL/OGR source distribution. I found this documentation how to include DWG into ogrinfo. But I do not know how to build an executable from a source code. So I am looking for a “ready to use” distribution of gdal (ogrinfo) with DWG support.

Thanks

Jens


texts: Plotting US maps with selected locations (longitude and latitude of locations are known) in R?

texts:
I've been struggling with the following R plotting of the US map and wonder if any of you have suggestions on how the plot should be done.

Here is what I have tried:

US <- map("state",fill=TRUE, plot=FALSE)
US.names <- US$names
US.IDs <- sapply(strsplit(US.names,":"),function(x) x[1])
US_poly_sp <- map2SpatialPolygons(US,IDs=US.IDs,proj4string=CRS("+proj=longlat + datum=wgs84"))
plot(US_poly_sp,col="white",axes=TRUE)

With the above code, I can add points if I need to. For example, suppose I have a location with lat=35 and lon=-110, then I can simply use the next line since the latitude and longitude of the spatial polygon I created from the object US contains information on the location:

points(lon,lat)

How can I get a plot of the map that is curved, as generated from the following line and also add points which are given in terms of their latitude and longitude?

US <- map("state",proj="albers", par=c(30,40), plot=TRUE)

I know it has something to do with the projection of maps, but the moment I projected the map, the corresponding range of the latitude and longitude (US$range) changes.