texts: Connecting to PostGIS database from ArcMap for display and query without ArcSDE?

texts:
How can I connect to a PostGIS database from ArcMap using ArcGIS Desktop 9.3 and later?

I would like to be able to perform spatially enabled queries and receive the results back (e.g. spatial and non-spatial joins, filtering etc.) rather than just dumping the contents of a table.

I don't want to use the ArcSDE spatial extensions, I want to use the PostGIS spatial extensions in ArcGIS Desktop.


texts:
Background: I'm currently trying to dynamically change the ranges for a layer with new style graduated symbology.

Problem: I thought it may be as simple as updating the ranges for the existing symbology using the setLowerValue() and setUpperValue() methods on the QgsRendererRangeV2 class as below:

layer = qgis.util.iface.activeLayer()
rangePoint = 0
rangeIncrement = 100
for range in layer.rendererV2().ranges():
        range.setLowerValue(rangePoint)
        rangePoint += rangeIncrement
        range.setUpperValue(rangePoint)

By my reconning the API documentation suggests that this should be possible.

However, I am getting the error:

AttributeError: 'QgsRendererRangeV2' object has no attribute 'setLowerValue'

I've tried also using:

range.setLowerValue(float(rangePoint))

with no luck.

What am I doing wrong?


texts:
Is there a way to use OSRM with PostGIS (Windows)?

Related: Why mobile routing is so fast?


texts:
I want to change symbology of one layer in ArcGIS 10.3. I want to have two types of symbology, but the field I want to use is string and there are more than 500 combinations, for example:

a1_b1_b2
a2_b2
b1
a1_b2_c3
...
a1_m6
...
a1_x1

I would like to have all the strings containing "x1" with another type of symbolgy, but how can I do it?


texts: How is Band 8 (panchromatic) of Landsat connected to Surface Reflection/Albedo?

texts: I want to work with the surface reflection/Albedo of Landsat 8 and got told to use B8 which is the panchromatic. It catches all the visible light so it makes sense to me that it is connected but I know there is no surflace reflection for the panchromatic band so I was wondering how is it connected and how can I use B8 for my SR work?


texts: QGIS Add PostGIS layer: restricting to single schema?

texts:
My primary PostGIS database has 12 schemas (schemae? TNIRN...), some of which have nasty-sized tables. The DB is on localhost, on my own machine.

For the most part, the actual working tables are in one schema ('aa') - and are at least 1, and often 2 or 3, order[s] of magnitude smaller than the nasty-sized tables.

Every time I want to add a PostGIS layer, the Layer->Add Layer->Add PostGIS layers tool makes me wait while (10 minutes) it adds every.single.schema and every.single.geometry-column to the selection universe, when I really want it to stop when it's finished loading schema 'aa'.

Top marks to QGIS for finding all the geometries in a pretty massive database, but in the ordinary course of events I am only interested in 'aa'.

The name 'aa' is chosen specifically because it does not overlap in any way with names of the 'big-ugly' shemas (schemae? TNIRN). When I use 'Search Options' and tell it 'Search in columns'=Schema for the search term 'aa', it still loads every table.

Is it possible to restrict the domain of 'Add PostGIS Layers' to only look at one schema? I've R'd TFD to no avail.

inb4 "have a different database for your working schema". Not possible, since there is beaucoup interaction between 'aa' and other schemas (schemae? TNIRN) when actual SQL-ish work is to be done.

Versions: PostgreSQL 9.3.5, PostGIS 2.1.3, QGIS 2.6.1-Brighton.


texts: Changing raster resolution using ArcGIS for Desktop?

texts:
My raster is displaying at a very low resolution so it is almost useless.  Last time I used this file it was at a much higher resolution, so I know the information is there.

To change it I first wanted to know what it was so I googled how to display it and read the Esri documentation on 'Displaying the raster's spatial resolution:' right-click for Properties, go to Display tab, then it said to check a box that doesn't exist.

Any ideas?

I'm using ArcGIS 10.2.2 for Desktop.


texts: Software that when input any angle, any CRS, outputs custom CRS that produces maps rotated by the angle

texts:
What open source GIS application, or plugin for same, or separate utility (GUI preferred but command line OK) has as part of its functionality the ability to take any angle (and a center of rotation) and any CRS as an input, and output a custom CRS that creates the same maps only rotated by the angle?  (I have found similar questions here and elsewhere but the answers involved ESRI software and/or were limited to UTM.)

In my specific application I am hoping to create for a QGIS project a custom CRS that produces rotated 67.5 degrees clockwise versions of maps produced with EPSG 102009 (or 102003, both conic).  (If you want to know why see Projections for solar eclipse umbral paths)

(Rotating the canvas would be an alternative approach but there is a bug that throws the labels for the layers way off.)


texts: Point cloud not normalized, how do you normalize with lidR?

texts:
I am using terrestrial LiDAR data, I ran the las_check() function on a file. I get check marks and everything is good until the normalization in the preprocessing. Is there an R code to do the normalization?
 Checking preprocessing already done
  - Checking ground classification... yes
  - Checking normalization... no
  - Checking negative outliers... ✓
  - Checking flightline classification

I am attempting to take terrestrial LiDAR measurements and use package lidR to create a DTM and voxelize my data so I can use package forestr to give me additional metrics I can't get out of lidR package. I have read several papers that seem to create a DTM and voxelize easy. I read the lidR book https://jean-romain.github.io/lidRbook/io.html#asprs-compliance  but I am still new and feel like a dummy when it comes to all this.


texts: Column 'SHAPE' cannot be null in mysql database even if the query string is correct

texts:
I am trying to write an POLYGON into an column named SHAPE which is a geometry type. And created query which looks

$sql = "INSERT INTO rokka(zone,district,location,area,SHAPE)
VALUES ('$zone','$district','$location','$area',PolygonFromText('$polywrite'))";

The variable $polywrite comes as

var polygon = 'Polygon((' + northeastlngLat[0] + ' ' + northeastlngLat[1] + ',' + northwestlngLat[0] + ' ' + northwestlngLat[1] + ',' + southwestlngLat[0] + ' ' + southwestlngLat[1] + ',' + southeastlngLat[0] + ' ' + southeastlngLat[1] + ',' + northeastlngLat[0] + ' ' + northeastlngLat[1] + '))';
$polywrite = "<script>document.write(polygon)</script>";

The arrays contain the coordinates.
And when executing

if ($connect->query($sql) === TRUE) {
    session_start();
    header('location: admin.html');
    session_destroy();
} else {
        echo "Error: " . $sql . "<br>" . $connect->error;
    }

The error it shows is

Column 'SHAPE' cannot be nullError: INSERT INTO rokka(zone, district,location,area,SHAPE) VALUES ('Rapti', 'Rolpa','Liwang','603750', PolygonFromText('Polygon((85.06290725100276 28.00935569163021,85.06392404943365 28.0093478254886,85.06393735099435 28.01070145242192,85.06292053985693 28.010709319009717,85.06290725100276 28.00935569163021))'))
Column 'SHAPE' cannot be null

But whenever I run the query

INSERT INTO rokka(zone, district,location,area,SHAPE)
VALUES ('Rapti', 'Rolpa','Liwang','603750', PolygonFromText('Polygon((85.06290725100276 28.00935569163021,85.06392404943365 28.0093478254886,85.06393735099435 28.01070145242192,85.06292053985693 28.010709319009717,85.06290725100276 28.00935569163021))'))

in my database directly It works exactly as it should. All the values are written into the respective columns. Please help me figure out what i am doing wrong.


texts: QGIS 3.4 feature filter MAX(), GROUP BY does not work

texts:
I thought my problem was solved with the question I asked bevor (QGIS 3.4 feature filter maximum() does not work). But there is still a problem.
I still want to filter a layer with the maximum value of field_a, grouped by field_b. I used this filter: "field_a" IN SELECT MAX("field_a") FROM 'layer_test' GROUP BY "field_b") The result is a reduced number of features, but there are still several features from one value in field_b with different values in field_a. Here a picture from the wrong result:

But I want only the maximum value of field_a, grouped by field_b. For example I want for value '12' in field_b only the feature with value '6' in field_a. If I use the expression without GROUP BY, then it works like I expected, I get only the total maximum value. How can I get the right GROUP BY result? Unfortunately I don't know SQL...


texts: Should I be using the OpenGeo Suite SDK or GeoEXT?

texts:
Im new to the OpenGeo Suite, and was wondering if I should create my apps using the SDK Link  Here Or I should use openlayers and GeoEXT Link here

I believe the SDK uses openlayers / geoext but makes the creation of map applications user friendly? I just need to know what to use each for in the terms of a new user

EDITS:

I want to create an app that displays layers, I have done that with the SDK pretty easy with GXP. However I also want to add graphs/charts to the maps, this is where I am unsure which route to go down. The user would click on a point and would be able to view a chart of that points change over time.


texts: Cannot Find Vertices Counter in Plugin QGIS 3

texts: I have recently installed the new version of QGIS and can't seem to find the vertices counter in the plugin menu to install. Does anyone know where I can find it?


texts: QGIS is crashing when I click on the SVG fill function

texts:
This problem is new from the last two or three days. When I click on the SVG fill in the Layer Styling tool bar

I get this message in upper left corner

and after few seconds I get this other message

and then QGIS crushes. I tried to re-install, but it did not change this problem. The QGIS software used to crush from time to time, but never this often and never often. I am working with the 2.18.9 (64 bit) version that I download from the QGIS web site.

Can anyone help me with this topic?


texts: WFS Request against a layer with a foreign table join

texts:
I have a Geoserver layer that is pulling from a Postgresql view that contains a left join to a foreign table (data from a SQL Server database).  If I query against the view in PGAdmin (ie. SELECT * FROM myView WHERE myID = 5) the correct results are returned.  But if I make a WFS request against the layer (ie. ...geoserver/READ/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=myViewLayer&CQL_FILTER=myID=5) the results include ALL records from the joined foreign table, as if the view contained a CROSS JOIN.

The foreign table in question is pulling from a SQL Server instance set up with a tds_fdw foreign data wrapper.  If I actually insert the records from the foreign table into a postgresql table and use this new table in the join instead then my WFS request returns the correct results.

Here are my software versions:
GeoServer:  2.15.2;
Postgresql:  9.6.18;
PostGIS:  2.54;
tds_fdw (foreign data wrapper): 2.0.1

Any thoughts?


texts:
Which .NET mapping components are currently available? Focus is on desktop use and fairly basic mapping of SHP files (thematics etc.), not the more advanced GIS functionality (reprojections, clipping, point-in-polygon etc). Royalty free distribution is needed.

I am aware of TatukGIS, AspMap, Map Suite and MapXtreme. I believe ESRI has nothing to offer in this area any longer, now that mapobjects is gone.

Other options?
Any stable open source offerings available?


texts:
I have a problem with creating - inserting GPS coordinates into QGIS. I have my coordinates in this format: N 45.66209; E 013.99382 and I want to insert them. I am using the function Numerical Digitize.

The problem is that when I put them into the function, the location that I get is wrong.
I think the problem is in the format of the coordinates, because when I used the function Location (X,Y) I obtained coordinates in this format: 421949.975,58031.353.

How can I transform into this format my coordinates? How is this format called?


texts: How to let the right toolbar of CartoDB show out? My page window is different from that in the guiding video.


texts: Is it advisable to limit privileges in `public` schema of PostGIS database?

texts:
When a PostgreSQL database is spatially enabled with PostGIS, PostGIS creates metadata tables (spatial_ref_sys), views (geometry_columns), and many, many functions in the public schema. As is usual practice, I keep actual data tables outside of the public schema.

I use PostGIS to teach a course on spatial databases. The students necessarily have privileges to much of the database, although they can only alter data or create objects in certain schemas. By default, the public role has full access to the public schema. It occurs to me that this should be curtailed. Not only do I want students to not make inappropriate changes or deletions to data and objects in this schema, I would also like their "experiments" to be contained within particular schemas (named for a course code, or to have one scratch schema for all courses), so that I can easily drop an entire schema at the end of the semester.

However, I have not seen discussion before in PostGIS or PostgreSQL mailing lists or forums about limiting privileges in the public schema, and I am concerned about the impact. My questions are:

Is it advisable to limit privileges in public schema of PostGIS database?
If I do so, what would be the correct set of privileges and/or default privileges that would keep all (existing and future) tables and functions in the public schema available to all users, but unable to be altered in their data or their definition?

On various servers I am running PostGIS 2.1 or 2.2 on various flavors of Linux, but the question applies to all versions and platforms.


texts: What's the specific meaning of ArcMap "FDO error: -2147154947"

texts:
ArcMap 10.3.1 is giving us an error of "Unable to save edits. FDO error: -2147154947".  This is not contained in the ESRI list of FDO errors: http://resources.esri.com/help/9.3/ArcGISEngine/ArcObjects/esriGeoDatabase/fdoError.htm

Here is some more information about our environment and under what circumstances the error occurs:

Once a user receives this error, they will continue to receive this error on every subsequent attempt to Save Edits.  The only way we know to fix the problem is to restart ArcMap and lose all edits since the last save.  The database we are using is an ArcGIS 10.3.1 enterprise multi-user geodatabase hosted on MS SQL Server 2012.  The client software we are using to access the database is ArcMap 10.3.1.  There are multiple versions in the database, and each version is being edited by a single user.  We can have as many 16 users editing different versions simultaneously.  Users will typically save edits every 15-20 minutes, and then reconcile and post their edits a couple times a day.  Archiving is enabled on the GDB feature classes, as well as editor tracking.  We have determined that this error number is probably an HRESULT, but we are unable to find any documentation from Esri or Microsoft for this particular error.

What is the specific meaning of this error code and message. In other words, what specifically is preventing ArcMap from saving the edits?

Edit:
The error happened again and this time we were able to retrieve the SDE log file from the client machine. I have substituted "client_machine_name", "database_name" and "feature_class_name" for the actual values for clarity.

  [Fri Feb 24 13:57:07 2017] [619] [client_machine_name] DB_registered_table_trim_state (database_name.DBO.feature_class_name,38063) Execute Immediate error (-1021).
  [Fri Feb 24 13:57:07 2017] [619] [client_machine_name] Error when updating kept rows in geometry table.
  [Fri Feb 24 13:57:07 2017] [619] [client_machine_name] Warning: Unable to trim state 38056->38063 for table database_name.DBO.feature_class_name, error = -1021


texts: How to remove subgroup labels from legend?

texts:
In print composer legend is it possible to turn off the display of the group name in the legend?

So it shouldn't show the DSC_Roads and other items marked with X


texts: Spatial overlays: Union between multiline and multipolygon with geopandas

texts:
I try to overlays (apply union) my multipolygon on the multiline with GeoPandas, but it seems like the GeoPandas overlay function work only with (Multi)polygon. Is there any way to make that with multiline and multipolygon with GeoPandas. Here is my code:

my_shp = gpd.read_file('my_shp.shp') # it's MultiLines shapefile
shp_buf = my_shp.buffer(5)
# as shp_buf it geoSeries i need to convert it to GeoDataFrame
df_shp_buf = gpd.GeoDataFrame({'geometry' : shp_buf})
union_1 = overlay(my_shp, df_shp_buf, how = 'union')

And I got this error:

TypeError: overlay only takes GeoDataFrames with (multi)polygon geometries


texts: Any web mapping solution for mapping points and manage 1:n relationship?

texts:
I need to show on a web map some points (> 12.000) ... I know that there are many solution for that, also for free (CartoDB, MapBox, ArcGIS On Line, ecc ... ), starting from some shapefiles. That is not my problem

My problem is the follow: I've a table with this structure

TABLE_1: key, x-coord, y-coord, attr-1, attr-2 .....attr-n

and another table with this structure

TABLE_2: key, description

Selecting a value for KEY in TABLE_2 there are several rows in TABLE_1 for that KEY value (1:n relationship).

I'd like to select a value in TABLE_2 and show (or better, point out ...) the points of related rows in TABLE_1 on my map.

I'd like to find a solution without writing "ad hoc" code because my customer is really entry level on GIS or web mapping and I've to keep as more simple as possible my solution.

Any suggestion or example about a situation like this one?


texts:
I don't know what's going wrong when I try to join a shapefile with an .xls or CSV. I do have the same identifier in both files. The joined file has all the expected columns but they are all empty.

What am I doing wrong?


texts: Why does the area of overlap with st_intersection not match my expectation?

texts:
I'm using the protected area shape files from Protected Planet for Swaziland/ Eswatini.

I calculate a home range (95% kernel density) and want to measure the overlap in area between this and the protected areas. When I do so the area of overlap is greater than the area of the home range which is not possible. I don't see where I'm going wrong.

Here's my code:

# load packages
library(tidyverse)
library(amt)
library(sf)
library(spex)

#' load in the protected areas
merged_Africa = read_sf("shape//eswatini//WDPA_Apr2020_SWZ-shapefile-polygons.shp")
st_crs(merged_Africa) <- 4326
st_crs(merged_Africa)
merged_Africa_tranform <- st_transform(merged_Africa, "ESRI:102022")

#' location data
x_ <- c(707692, 707589, 707998, 708407, 708916, 709415)
y_ <- c(-3030991,-3031423,-3031640,-3031750,-3032508,-3037158)
mydata <- data.frame(x_, y_)

# transform to trk
trk <-
  mk_track(mydata,
           .x = x_,
           .y = y_,
           crs = CRS("+proj=aea +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))

#' calculate home range area
kde_shape_1 <- trk %>%
  hr_kde(., levels = c(0.95))
hr_area(kde_shape_1) / 1e6 # 52.00846

#' convert raster to polygon
psf <- qm_rasterToPolygons(kde_shape_1$ud, na.rm = TRUE)

#' measure the intersection
intersection <-
  st_intersection(psf$geometry, merged_Africa_tranform$geometry)
sum(st_area(intersection)) / 1e6 # 61.21044


texts: How to fix Cannot create PoolableConnectionFactory (FATAL: password authentication failed for user )

texts:
I am working on project using GeoServer, PostgreSQL and docker. I am newbie in all these technologies.
I publish new layer in GeoServer and when I try to preview, my browser tries to download instead of showing preview as it does before.
I had the idea to delete container and recreate with pulling the last image again.
Once done, nothing not working anymore. And I lost some data. I check I get missing styles, and I recreate some of them. The issue now is I do not be able to see my preview. I can see only the last layer I added, but when I check in GeoServer workspace data folder, I see my layers folders but it seems that something is missing. I re-check and I get this error.
java.lang.RuntimeException: Error getting FeatureType, this should never happen!
Error getting FeatureType, this should never happen!Unable to obtain connection:
Cannot create PoolableConnectionFactory (FATAL: password authentication failed for user "root")
Cannot create PoolableConnectionFactory
(FATAL: password authentication failed for user "root")
FATAL: password authentication failed for user "root"

So I am in trouble, I do not understand why I am getting this. Is it because I remove the container? I was thinking that I must get my data because there is a volume option in docker-compose file for my system?
So can I fix everything if I copy data from another container?


texts: Executing Python script for GRASS from PHP

texts:
I wrote the following Python script:
#!/usr/bin/env python
import os
import sys

gisbase = os.environ['GISBASE'] = "/usr/lib/grass64"
gisdbase=os.path.join("PATH")
location=LOCATION NAME
mapset=MAPSETNAME
sys.path.append(os.path.join(os.environ['GISBASE'], "etc", "python"))

import grass.script as grass
import grass.script.setup as gsetup
import grass.script as grass

gsetup.init(gisbase, gisdbase, location, mapset)
print "one"

grass.run_command('v.in.ascii','-o', input='/opt/lampp/htdocs/project/basic_statistics/sites_grass.csv', output='sites',fs=",",skip='1', x='9', y='10')
print "two"

grass.run_command('v.db.addcol', map='sites', columns="height INTEGER")
grass.run_command('g.region', rast="dem21")
grass.run_command('v.what.rast',vector='sites', raster='dem21', column='height')
print "three"

grass.run_command('v.out.ascii','-o', input="sites",output="sites_geography_stats.csv", columns='height')
print "four"

It runs without any problems when I simply execute it from UNIX terminal (/.FILENAME.py). However, my goal is to execute it from a PHP code:
<?php
...
echo system("./FILENAME.py");
...
?>

When I run the PHP code the following output is printed to the website:
one two three four
Sadly, the sites_geography_stats.csv file is not created. It looks that the python script is executed from php but its GRASS-related part is ignored.
Does anybody have any ideas how to fix it?


texts: Copying coded value for field and not raw value?

texts: I am working with ArcGIS 10 and I have added a field in a feature class and am attempting to copy the text from another field's coded values, but the field calculator is simply copying the "code" text from the source field and not the "value" text from the domain rule.


texts: QGIS 2.12 Edit rendering styles for outputs not listing all outputs

texts:
I have built a model in QGIS 2.12 and have been using the Edit rendering styles for outputs to set the styles for the models output data. But I have come across an issue.

When I am configuring a particular model one of the outputs are not listed. As shown in the image, there is a line for the particular output but nothing else. Has anyone come across this issue before? If so what was your fix? if not does anyone have any pointers to get to the bottom of this?


texts: Hardwire definition query on layers in mxd VB.net

texts:
I am creating a tool that applies definition query to certain layers of the mxd. The user clicks the button and it should hardwire the definition query into the layer. This is useful if there are many layers you want to apply the same def query onto. When I click the button I get the following error message

Obviously my code has flaws in it that would make such drastic error message.

Public Sub DefQUpdate_HardWired()
    Try
        Dim pMxd As IMxDocument
        pMxd = My.ArcMap.Application.Document
        Dim pMap As IMap
        pMap = pMxd.FocusMap
        Dim pEnumLayer As IEnumLayer = pMap.Layers
        pEnumLayer.Reset()
        Dim pLayer As ILayer = pEnumLayer.Next

        Dim pFlyr As IFeatureLayer, pFlyrDef As IFeatureLayerDefinition
        pFlyr = Nothing
        pFlyr.FeatureClass = pLayer

        Do While Not pLayer Is Nothing
            If TypeOf pLayer Is IFeatureLayer Then
                pFlyr.Name = pFlyr.FeatureClass.AliasName
                pFlyrDef = pFlyr
                If pFlyr.Name = "UPL" Then
                    pFlyrDef.DefinitionExpression = "LASTEDITED < 'TO_DATE('2013-12-03 09:09:29','YYYY-MM-DD HH24:MI:SS')'"
                End If
            End If
            pLayer = pEnumLayer.Next
        Loop
        pMxd.ActiveView.Refresh()
    Catch ex As Exception
        MsgBox(ex.Message)
    End Try
End Sub


texts: Cloudless days count with MODIS in Google Earth Engine

texts:
I'm trying to estimate the total of clear days in a year with the help of state_1km band in the MOD09GA product. But I'm not sure if it is right, also the possible counts are higher than expected when compared with cloud cover observation at weather stations.

  // Load morning (Terra) MODIS data.
        var morning = ee.ImageCollection('MODIS/006/MOD09GA');

        // A function to extract cloudiness data from a MODIS image.
        function clear(image) {
          // The cloud state is in the low two bits of the 'state_1km' band.
      var cloudState = image.select('state_1km').bitwiseAnd(0x03);
      // 0:clear, 1:cloudy, 2 means mixed, and 3:missing data.
      return cloudState.float()
         .remap([0,1],[1,0])
         .updateMask(cloudState.neq(3))
        .updateMask(cloudState.neq(2));
     }

    var select_morning = morning.filter(ee.Filter.calendarRange(2015,2015,'year'));
    var clear_days = select_morning.map(clear);

    var count = clear_days.reduce(ee.Reducer.sum());


texts:
I am using VueLeaflet to create a map like in this example: https://vue2-leaflet.netlify.com/examples/geo-json.html

I like to add some colors to the polygons like in this example:
https://leafletjs.com/examples/choropleth/

This is my try here: https://codesandbox.io/s/old-framework-mlc7b

 getcolorFunction(d) {
      return d > 1000
        ? "#800026"
        : d > 500
        ? "#BD0026"
        : d > 200
        ? "#E31A1C"
        : d > 100
        ? "#FC4E2A"
        : d > 50
        ? "#FD8D3C"
        : d > 20
        ? "#FEB24C"
        : d > 10
        ? "#FED976"
        : "#FFEDA0";
    },

    styleFunction() {
      const fillColor = this.fillColor; // important! need touch fillColor in computed for re-calculate when change fillColor
      return () => {
        return {
          weight: 2,
          color: "#ECEFF1",
          opacity: 1,
          fillColor: getColor(feature.properties.sales),
          fillOpacity: 1
        };
      };
    },

    onEachFeatureFunction() {
      if (!this.enableTooltip) {
        return () => {};
      }
      return (feature, layer) => {
        layer.bindTooltip(
          "
Name Province:" +
            feature.properties.varname_1 +
            "

Sales: " +
            feature.properties.sales +
            "
",
          { permanent: false, sticky: true }
        );
      };
    }
  },

How could I Get work that function getColor in Vue.js ?


texts: Free download of Soviet topographic maps (also known as Genshtab maps)

texts:
In the past, http://poehali.org/maps/ offered a free-of-charge download of Soviet military topographic maps of many areas in the world. Those maps are in the public domain. Now, this site redirects to http://mapstor.com , a website offering a low-cost archived download of those maps, but not free-of-charge. http://www.topomapper.com/ allows for interactive looking at Soviet military maps, but only for parts of the world and only up to a scale of 1:100.000, even for areas mapped at 1:50.000.

Where can I download free-of-charge public domain Soviet military maps?


texts:
I'm quite okay at R programming, but fairly new in GIS (especially in R).

I have a set of coordinates and I'd like to create a custom shape around them.
What I want to do step-by-step:

Create a circle around the points and slice it to couple of sectors. The radius and the number of wedges are given manually at this point.
Group the points by the wedges
Calculate the radius of the given wedge based on the points characteristic of the points in the sector (I have the model for this)
Set the radius sector by sector.
Export the wedges as different regions of a polygon shapefile

Unfortunately, I'm stuck at the very beginning. I created a standard pie chart by using add.pie {mapplots}, but I wasn't able to convert its sectors to spatial objects.


texts: How to add more vertexes to a LineString

texts:
I created a map of LineStrings using QGis and then it was exported to a Postgis DB, then I used the query below in order to create the vertexes in the map so it could be used with the functions of PgRouting.

select pgr_createTopology('ways', 0.00000001, 'geom', 'gid');

The following map is the result of above procedure.

How can I add more vertexes in the LineStrings?


texts: What signal can be used to get the selection change on a feature?

texts:
Right now I am using:

self.iface.mapCanvas().selectionChanged.connect(FunctionThatReactsToTheSignal)

But this one sends two signals, one for the deselection of the previous feature and one for te selection of the new one, so the FunctionThatReactsToTheSignal is executed twice, is there a way to make it send just one signal for the new selection or another signal function to avoid repetitions?


texts:
How do I get all values of a specific field of a feature class in ArcObjects (c#.net)?
Is there any way to select all features using IQueryFilter interface?


texts: QGIS android issues with editing points on lines and polygons

texts:
I am using QGIS on an android tablet and cannot delete mistaken points in a line or polygon. Does anyone know how to do this?

I also have an issue with the first point of a line/polygons starting not where I put the first tap when I start a polygon, but where the last point was on the previous polygon - resulting in the first section being incorrect. How to get around this issue?


texts: How to Create Custom ZoomIn Button that is Windows Tools Button with Custom Cursor?

texts:
I want to create large zoomin Button without toolbarbutton. I created a 'Zoomin Script' on axMapContorol. It worked is fine but I can't set cursor.
My code is as follows:

IMapControl3 mapControl = (IMapControl3)FMap.Instance.axMapControl1.Object;
MemoryStream ms = new MemoryStream(LargeZoominCursor);
Icon icon = new Icon(ms);
mapControl.MouseIcon = (stdole.IPictureDisp)ESRI.ArcGIS.ADF.COMSupport.OLE.GetIPictureDispFromIcon(icon);
mapControl.MousePointer = ESRI.ArcGIS.Controls.esriControlsMousePointer.esriPointerCustom;

This code doesn't call errors during the first time only. After checking SelectFeatureTool on ToolbarControl and pushing my custom zoomin botton, the above script doesn't work anymore. The cursor does not change.

The 'SelectFeatureTool' and 'Zoomin' still work with 'SelectFeatureTool' cursor.

How to create change cursor?


texts: Automatically discard nodata area of raster

texts:
Given a raster with an overly large nodata surround, how can I automatically discard the surround? Said another way: how does one shrink the raster extent to only occupy the area containing values?

I know one can interactively draw a box and clip the raster, or zoom and export only visible extent, or set an analysis mask. I'm after a geoprocessing tool or simple python script. Tools at hand are ArcGIS with Spatial Analyst, QGIS, and GDAl/OGR, but really anything is fair game.

The example raster is 16bit signed integer, with -32768 defined as NoData, but the null value or bit level really shouldn't matter.


texts: Adding tree height values to raster / vector data in ArcGIS 3D Analyst?

texts:
I am new to working with raster data and I am trying to teach myself a few things. I have what I think is a simple problem but I can't seem to figure this one out after a few hours of searching Arc Help etc.

I have raster with ground elevations to which I want to add tree heights. I have a Forest Cover polygon vector layer which has tree heights for my study area. I would like to add the heights of trees to the elevation data.

My first guess was to convert the Forest Cover polygon to a raster based on the tree height field and then Merge/Combine the two rasters?

I have esri's 3d Analyst so I think I would have all the available tools?
I am using ArcMap 10.0


texts: Why the line in the shapefile does not align with the base map reference Runway?

texts:
I created a shapefile and set the coordinate system to WGS 1984 UTM zone 45N, and I added the base map Imagery with Labels.

I created a line feature with Coordinates of two ends of the runway

(27°41'1.694"N  85°21'12.202"E)
and (27°42’30.058’’ N  085°21’52.283’’ E)

The line correctly aligns the runway center line, but the direction of the line does not match with the actual direction of the runway .

Example: When I draw a line with the runway length 1339.5m and heading 22°  from  point(27°41'1.694"N  85°21'12.202"E), the line does not aligns with the runway centerline.

the direction of the runway is 22° from the (27°41'1.694"N  85°21'12.202"E)


texts: Writing Python code to create intersection loop with geodatabase?

texts:
I am trying create a loop that will intersect each individual township with a statewide section layer. I am having trouble with the syntax of the intersect tool I think, but maybe it is something else.

# Import system modules
import sys, string, os, arcgisscripting
import arcpy, os, sys, string
from arcpy import env
from arcpy.sa import *
gp = arcgisscripting.create()

arcpy.OverWriteOutput = 1
arcpy.env.overwriteOutput = True

input_folder = "C:\\state\\towns.gdb"
arcpy.env.workspace = input_folder
out_wksp = "C:\\state\\township_sections.gdb\\"

env.workspace = input_folder

fcl = arcpy.ListFeatureClasses()
print fcl

infeatures = ["C:\\state\\towns.gdb\\sections", fcl]
print infeatures

#ExtractByMask_sa (in_raster, in_mask_data, out_raster)
for fc in input_folder:
    print fc
    arcpy.Intersect_analysis(infeatures, out_wksp + fcl, "ALL", "", "INPUT")
    #arcpy.Intersect_analysis(["C:\\state\\towns.gdb\\sections" fcl], out_wksp + fcl, "ALL", "", "INPUT")
#    fc = fclist.Next()
    print 'clipped', fcl

print 'done'


texts: How to convert all shapefiles in folder into KML using ArcPy?

texts:
I was trying to convert all the shp in a folder into kml.

featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
 # Set Local Variables

 composite = 'COMPOSITE'

 pixels = 1024

 dpi = 96

 clamped = 'CLAMPED_TO_GROUND'

 scale = 1

 outKML = fc[:-4] + ".kmz"

 arcpy.LayerToKML_conversion(fc,outKML, scale, composite,'', pixels, dpi, clamped)

It always says
Failed to execute. Parameters are not valid.
ERROR 000732: Layer: Dataset ZZZ.shp does not exist or is not supported
Failed to execute (LayerToKML).
But I can manually do it within ArcMap 10.1 Desktop...


texts: Create event handler for "x" button on type: 'floating' widget?

texts:
How do you create event handler for "x" button on type: 'floating' widget?

"OnClose" fires every time: this.parentWidget.hide(); is called.

I want to intercept the event when the "x" button in the title bar is clicked.


texts:
I am new to QGIS and most of the tutorials I have found are for 1.x and the overhaul has me confused.

I would like to create R scripts and run them in QGIS. I have read here that it's possible using the SEXTANTE plugin.

It seems that SEXTANTE is part of the core install now and that it's been renamed to processing. I've found out how to make python scripts but not R.

Is it  possible to create and run custom R scripts in QGIS 2.0? If so, where do I go in the new processing function?


texts: Bypassing schema lock for a join on a feature class published to service?

texts:
I would like to programmatically join a table to a feature class that is published as a service on ArcGIS Server 10.4.1.  The problem is that the schema is locked.  How can I bypass this error?  When schema locking is disabled in the server manager, the error still appears.

It is not an option to programatically turn the service off while geoprocessing tasks are executed on the service.

The code below outlines where I would like to execute my join; the schema lock error is thrown at the join line.

def calculate_overdue_tickets(self):

    self.logger.info("calculating overdue days for feature classes '{0}' and '{1}'".format(self.line_ticket_feature_class, self.poly_ticket_feature_class))

    expression = "getDaysOverdue(!{0}!)".format(self.begin_date_field)
    code_block = """def getDaysOverdue(begindate):
import datetime
if begindate is not None:
    begin = datetime.datetime.strptime(begindate, '%m/%d/%Y').date()
    return (datetime.date.today() - begin).days
else:
    return None"""
    arcpy.management.CalculateField(in_table=self.line_ticket_feature_class, field=self.days_overdue_field, expression=expression, expression_type='PYTHON_9.3', code_block=code_block)
    arcpy.management.CalculateField(in_table=self.poly_ticket_feature_class, field=self.days_overdue_field, expression=expression, expression_type='PYTHON_9.3', code_block=code_block)

    self.logger.info("calculating overdue status for feature classes '{0}' and '{1}'".format(self.line_ticket_feature_class, self.poly_ticket_feature_class))

    expression = "getOverdueStatus(!{0}!, !{1}!, !{2}!, !{3}!)".format(self.work_status_field, self.days_overdue_field, self.origpriority_field, self.nextjobid_field)
    code_block = """def getOverdueStatus(workstatus, days, origpriority, nextjobid):
if days is None:
    return '{5}'
elif days < 0:
    return '{0}'
elif days == 0:
    return '{1}'
else:
    if origpriority == '{2}':
        return '{3}'
    elif workstatus in ('N','P','S','A','H') and (nextjobid is not None):
        return '{3}'
    elif workstatus in ('N') and (nextjobid is None):
        return '{4}'
    else:
        return '{5}'
""".format(self.status_upcoming, self.status_today, self.origpriority_value, self.status_inprogress, self.status_overdue, self.status_none)
    arcpy.management.CalculateField(in_table=self.line_ticket_feature_class, field=self.overdue_status_field, expression=expression, expression_type='PYTHON_9.3', code_block=code_block)
    arcpy.management.CalculateField(in_table=self.poly_ticket_feature_class, field=self.overdue_status_field, expression=expression, expression_type='PYTHON_9.3', code_block=code_block)
    arcpy.JoinField_management(in_data=self.line_ticket_feature_class, in_field=self.poly_line_mobile_id, join_table= self.mobile_info_table, join_field=self.mobile_table_id)

Error:

EDIT: After a bit of testing, the join was executed correctly with the schema locking disabled in Server Manager.


texts: Using zonal statistics code from Python GDAL Cookbook?

texts:
I am new when it comes to GIS and relatively new when it comes to programming.

I was trying the “zonal statistics” recipe from the Python GDAL/OGR Cookbook, but I cannot get a decent result. For the code: click here. I know I can also use rasterstats, but I would like to build on this code for some other calculations.

There are 2 problems I am encountering with this code:

It is not taking into account all the features.
It is only calculating the average (and doesn't get any values for mean, median, etc.)

For my test, I am using a raster image of a DEM and a shapefile containing 3 polygons.

Here is my output:

{0: (25.960679824827782, masked, masked, masked, masked), 1: (25.960679824827782, masked, masked, masked, masked), 2: (25.960679824827782, masked, masked, masked, masked)}

I tried printing the top corners and the offsets of every feature it should take into account during the loop, but the problem is that it doesn't change feature (which it should with the GetFeature() method).

I found a similar question  as Zonal computations: computing raster pixel counts within polygons using GDAL versus rasterstats, but there was no real answer there.

I am not using arcMap or arcPy

Can anyone help me out?


texts: How to build Virtual Raster using grouped DEM layers as input with QGIS 3?

texts: I have compiled about 700 DEM tiles (each in its own layer) for some islands. I would like to create a virtual raster for each island. I grouped the DEM layers by the island name in the layer panel. But when I try to choose input layers, the groups are not available to choose from and I am forced to choose from the original 700 layers. Is there a workaround for this or a clever way to export the layers within the group so I can load them as files?


texts: Is it possible to modify the default style for imported polygons (etc) on Google Earth Pro?

texts: Google Earth has a default style setting for imported polygons that drives me crazy. They are set to 100% opacity, "white", with the random selector enabled, meaning they show up by default after being imported as a rainbow of randomly colored, opaque polygons. Is there any way to change this default style setting?


texts:
ArcGIS 10 SP3

I have read the other posts in this forum about updating metadata programmatically, and have examined Esri's help topics regarding xslt stylesheets to update metadata, but am hoping that someone has produced a more efficient solution in the last ten months:

I would like to perform something like the following examples:

for all layer files which contain the word "road" in their name, add the following list of words to their tags: road, street
for all SDE feature classes which do not currently have an Item Description, add the feature class name and tags into the Item Description
etc.

My question is: can I do this with python or VBA, or must I use xml/xslt stylesheets to do this?


texts: Can I do a downgrade just by over-installing , i.e.  deploying the war file of the older version, without losing settings?


texts:
I have several images from different times and for every time I have three types of image (temperatureA, temperatureB, pressure). Should I store them as

every time as band in a single raster with N bands => I will have (3) rows (one for every data type)
every time as a single raster with 1 band => I will have (times count * 3) rows (one for every data type and for every time)
all 3 types in bands of a single raster => I will have (times count) rows (one row for every time)

What is tbe best option?

I will mostly read from the database and update all data hourly. Users will be reading from PHP website values at [x, y] from a single time, but all three variables.


texts: insert update delete wfs layer to gesoerver on openlayers3

texts:
i am using geoserver 2.5.2 and openlayers3
i can show geojson wfs data from geoserver with geojson and edit with interactions on map

(function () {

    showData= function () {
        var format = new ol.format.GeoJSON();
        vectorLayer = new ol.layer.Vector({
            source: new ol.source.StaticVector({
                format: format,
                projection: 'EPSG:3857'
            }),
            //style: mapping.customStyleFunction
        });
        var select = new ol.interaction.Select();
        var modify = new ol.interaction.Modify({
            features: select.getFeatures()
        });
        // Create the map
        map = new ol.Map({
            interactions: ol.interaction.defaults().extend([select, modify]),
            target: 'map',  // The DOM element that will contains the map
            renderer: 'canvas', // Force the renderer to be used
            layers: [
                // Add a new Tile layer getting tiles from OpenStreetMap source
                new ol.layer.Tile({ source: new ol.source.MapQuest({ layer: 'osm' }) }),
                vectorLayer
            ],
            view: new ol.View({
                center: ol.proj.transform([28.86820725, 41.14670434], 'EPSG:4326', 'EPSG:3857'),
                zoom: 10,
            })
        });
        //map.on("")
    },
    setData= function (geojsondata) {
        var source = vectorLayer.getSource();
        var text = geojsondata;

        if (text === '') {
            source.clear();
            return;
        }
        var json = JSON.parse(text);
        var features = source.readFeatures(json);
        source.addFeatures(features);

    }
})

but I do not know how to save the updated data on geoserver using openlayers3
in other words how to post updated data  to geoserver?


texts: No change in graduated symbol size in map, the change occur only in layer panel

texts:
I want to make a map of India with the population in some cities, with symbol size ranging from smallest to largest according to the population value. I used graduated symbol. I have 13 categories of population and I customised the symbolize from 0.8mm to 15 mm. I want to make the size of large population symbol a bit larger than this. But when I change the largest symbol size from 15 to 25 or 35 no change occur in the map symbols. Only the size of symbols in layer panel is changing.
Is there any way to enlarge the symbol?
I am using QGIS version 3.18.1.
I am adding the screenshot of how I changed the size of symbol.
the layer styling panel looks like this.


texts: ESRI Clarity vs Google Maps satellite images, distortion due to elevation?

texts:
I have a couple of satellite images, each about 1 km square, of the same rock climbing area. One is from ESRI Clarity, the other from Google Maps. The Clarity image is projected to UTM and has a UTM grid superimposed, which is extremely convenient for my purposes, but it was taken in winter when a lot of the rock was in shadow or covered with snow. The Google Maps one is much better for showing the features I want to see and locate.

The Google image didn't register correctly when framed on top of the Clarity one, so I thought, OK, this is just because the projection Google is using isn't UTM. I decided to try to handle this by locating control points that I could identify on the two images, and distorting the Google image to make these points match. Surprisingly (to me, at least), it became clear that this could not be done accurately using an affine transformation, i.e., one that sends lines to lines. Below is my best fit with 6 control points and a polynomial fit. On the left is Clarity, on the right is Clarity with Google superimposed as an inset:

It seems to me that the registration of the two images at the edges is rather good. However, some features still seem extremely distorted in one image compared to the other. The easiest to see is the long dark line marked by the arrows, which is the shadow of a natural wall running down the north face of the rock. In the left-hand image, it's the dark line that you see if you trace out the red arrow and then continue on in the same direction. In the right-hand image, it's a little harder to see at this scale, but it's the same idea.

Now the two lines are not parallel at all. They differ by about 10 degrees. This is not just because I distorted the image with some goofy polynomial fit that might oscillate in an unrealistic way. In fact, the rotation of this feature is even more severe in the original Google maps image -- more like 20 degrees.

How should I interpret this? My guess is that ESRI made an effort to correct for perspective and the motion of the satellite as it scanned over the landscape, while Google did not, or didn't do it thoroughly enough. This might not be a problem in most cases, but this rock sticks up about 200 meters above the surrounding landscape and has steep walls. This would explain why I'm able to make things mesh fairly well on the surrounding terrain with a relatively gentle curvilinear transformation, but I still get crazy disagreement in the central area of steep terrain.


texts: What is the relationship between the SOC and the DBMS?

texts:
I have read lots of documentation about ArcSDE 9.3 recently but there are still are few things I don´t quite understand because I have not installed the infrastructure myself. If the SOM manages the connections to services and the SOCs host the services, what is the relationship between the SOC and the DBMS?

Is this Workflow correct?

Client application requests features from AGS
AGS communicates to the SOM that a service has been requested
SOM finds the SOC which is runs the service
SOC pulls data out of the DBMS and relays the data to SOM and AGS
AGS provides visualization and analysis of data


texts: Checking current user privileges for ArcSDE dataset with ArcPy

texts:
I am managing data distribution via an ArcSDE Oracle gdb to a number of users (schemas?), each of which can only have read only access to datasets they have requested access to, and have then been licensed to see - due to the confidential nature of the data.

Currently I keep these privileges in a spreadsheet as 0, 1, or 2 for each user, and with a Python program, I then update each dataset in the SDE to match the privileges in the spreadsheet using the ArcGIS "Change Privileges" tool.

However, this means I am currently running the change privileges on all SDE datasets, when in reality only a handful may actually need be changed, which is then taking a lot longer than it really needs to. So I am looking for a way to find out what privileges each dataset already has, to know if they do need to be updated or can be skipped, probably with a custom function (nothing Oracle please), as in a simplified example of my existing code below:

import arcpy

#  New function here
def CheckPrivileges(sde_fc, user):
    etc..
    return privileges_value  #  e.g. 'Select'

#  Input SDE and spreadsheet paths
sde = r'SDEConnectionPath'
in_xl = r'SpreadsheetSheetPath'

#  Dictionary to convert 0, 1, 2 values from spreadsheet to privilege lists
new_privileges = {"0": ["REVOKE",  "REVOKE"], "1": ["GRANT",  "AS_IS"], "2": ["GRANT",  "GRANT"]}

#  A second dictionary to check against the returned value from the CheckPrivileges function,
#  but using the same values from the spreadsheet e.g.
privileges_check = {"0": "", "1": "Select", "2": "All"}  #  (This is almost certainly incorrect)

#  In the spreadsheet, each user has a column with the 0, 1, 2 privilege values e.g.
user = 'USER1'

#  Search cursor with each dataset and the users privileges.
rows = arcpy.da.SearchCursor(in_xl, ['SDE_DATASET', user])
for row in rows:

    #  Path to SDE dataset having privileges altered.
    sde_fc = sde + '\\' + row[0]

    #  Pulls out the correct privilege values from the dictionaries.
    user_privs = new_privileges[str(row[1])]

    check_privs = privileges_check[str(row[1])]

    #  This is where I would call the new function I need e.g.
    current_privs = CheckPrivileges(sde_fc, user)

    #  And then check against what it should be according to the dictionary.
    if current_privs != check_privs:

        #  If not equal, then amend privileges, otherwise continue to the next row
        arcpy.ChangePrivileges_management(sde_fc, user, user_privs[0], user_privs[1])

How do I construct the missing CheckPrivileges function?


texts: How can ee.Image.stratifiedSample() be used in Google Earth Engine to select training samples for a binary classifier?

texts:
This is a follow-on question to How can one create a binary image from a FeatureCollection in Google Earth Engine?

The corresponding minimal example can be found here:
https://code.earthengine.google.com/3cbd9c722f8ceca84281327cecd8cdb8

I don't quite get how a raster image containing class labels for every pixel can be used to select training samples using GEE's in-built function ee.Image.stratifiedSample(). In the minimal example, the following error occurs:
"FeatureCollection (Error) Remote request too large (134217728 > 83886080) for output: [< Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >, < Object >]."

Although my question and the minimal example refer to a simply binary classifer (which could theoretically also be modelled as a one-class problem, as I am only interested in the foreground), the problem should be directly transferable to classical multi-class classification problems.

Maybe I am also misunderstanding the concept of ee.Image.stratifiedSample()? This possibly also relates to my other question What's the difference between sample, sampleRegions, and stratifiedSample in Google Earth Engine?


texts: Sum the values of a single raster and output the value as a number?

texts:
I have a floating point raster and have applied a conditional to it in order to only extract the positive values, so there are 'no data' "holes" in the raster. I simply need to sum the values from this raster and output it as a number.

I came across this thread, which I'm sure some of you will recognize. Has anyone developed a working solution? I'm working in Arc10 and making use of ArcObjects if necessary.

EDIT: I've included a VB.NET version of @Kirk's C# code. It's provided as is, but it should be a good starting point if it doesn't already completely satisfy its purpose.

    Public Shared Sub TestSumRaster()
        Dim path As String = "D:\topo_modeling\topo_modeling\gis\usgs\ned_10m.img"
        Dim rds = OpenRaster(path)
        Dim sum As Double = SumRaster(rds)
        Debug.Print("sum : {0}", sum)
    End Sub

    Public Shared Function SumRaster(rasDs As IRasterDataset) As Double
        ' sum up each band in a raster dataset
        Dim rasBandCol As IRasterBandCollection = DirectCast(rasDs, IRasterBandCollection)
        Dim sum As Double = 0.0

        ' Iterate through each band of the dataset.
        For i As Integer = 0 To rasBandCol.Count - 1
            sum += SumBand(rasBandCol.Item(i))
        Next
        Return sum
    End Function

    Private Shared Function SumBand(rBand As IRasterBand) As Double

' sum up each block in a band
        Dim sum As Double = 0.0
        ' QI to IRawBlocks from IRasterBandCollection
        Dim rawBlocks = DirectCast(rBand, IRawBlocks)

        ' Determine the tiling scheme for the raster dataset.
        Dim tile = New Tile(rawBlocks.RasterInfo)

        ' NoData is array, but usually just has one value in it.
        Debug.Print(rawBlocks.RasterInfo.NoData.[GetType]().ToString())

        ' Create the pixel block.
        Dim pb = rawBlocks.CreatePixelBlock()

        ' Iterate through the pixel blocks.
        For pbYcursor As Integer = tile.startY To tile.endY - 1
            For pbXcursor As Integer = tile.startX To tile.endX - 1
                ' Get the pixel block.
                rawBlocks.ReadBlock(pbXcursor, pbYcursor, 0, pb)
                sum += SumBlock(pb, pbXcursor, pbYcursor, rawBlocks.RasterInfo.NoData)
            Next
        Next
        Return sum
    End Function

    Private Shared Function SumBlock(pb As IPixelBlock, pbXcursor As Integer, pbYcursor As Integer, noDataValues As Object) As Double
        ' sum up each pixel in the pixelblock
        Dim sum As Double = 0.0

        ' Put the pixel block into a SafeArray for manipulation.
        Dim safeArray = TryCast(DirectCast(pb.get_SafeArray(0), System.Array), System.Array)

        Dim noDataList = New List(Of [Single])(DirectCast(noDataValues, [Single]()))

        ' Iterate through the pixels in the pixel block.
        For row As Integer = 0 To pb.Height - 1
            For col As Integer = 0 To pb.Width - 1
                Dim val As Object = safeArray.GetValue(row, col)
                Dim singVal = DirectCast(val, [Single])

                'Debug.Print(val.GetType().Name);
                If Not noDataList.Contains(singVal) Then
                    sum += singVal
                End If
            Next
        Next
        Return sum
    End Function


texts: Hang up when trying to create clipped .bmp from merged rasters .tiff

texts:
I'm a software guy (but don't know much about GIS).  I'm also a glider pilot in upstate NY.  A lot of the glider folks use a simulator called Condor and QGIS is used to help create photo realistic ground scenery....

They seem to have a hard time at it (just blindly following instructions)  I just started using it and was trying to create a step-by-step tutorial for them.  I'm using QGIS 2.18.28 on Windows 10/64 bit, 16 Gig of Ram.  The tutorial is at http://www.murtari.org/downloads/FLSC_SceneryMultiR1.pdf

It has more details, but here are the basic steps up to the hang. Starting with a brand new project:

Layer -> Add Layer -> Add raster layer. I select 6 tiff files of satellite imagery (each about 1.5 gig in size)
Layer -> Add Layer -> Add delimited text.  I have a .csv file with 4 calibration points, one in each corner.  The selected area is a rectangle 14 miles wide and 52 miles long.
I select the proper coord system: WGS84 EPSG:4326 - OK
Raster -> Misc -> Build virtual raster.  I choose the same tiff files.  OK (no load into canvas)
Raster -> Projections -> Warp(Reproject).  Chose the .vrt created above, Source SRS EPSG:4326, Target SRS our UTM zone 18N, resampling cubic spline (no load into canvas) - OK
Raster -> Extraction -> Clipper.  Chose the file from above. Clipping mode 'extend'. Dragged an area inside the calibration dots (load into canvas).  OK .  I see a red select area inside my dots.
Right click in Layers Panel over name of image above.
Select Save as...   Rendered Image, Create VRT, format GTiff, CRS to our zone 18N, uncheck saved file to map. No change to extent.
NOTE - Wanted hi res using 8192.
Resolution (current layer): Columns 8192, Rows 24576
Tiles:  8192 and 8192
Clicked OK - immediately see a small dialog box pop up, the title says "Saving raster...  ?" -- nothing else is visible.  QGIS hangs, windows says not responding, about 1.5 gig mem usage and stays there....

If the previous stepped work, we would get a tiff of the selected area which with other tools is converted to a .bmp and then merged into a landscape over textured terrain...


texts: Uncertain about the "n" in Getis Ord Gi* statistics

texts:
The formula of Getis Ord Gi* statistics is here. n is supposed to be the total number of features.

In the numerator of the Gi* statistics, it looks as though it is comparing (by subtracting) the weighted sum of the feature values around i and its neighbors against the overall weighted sum of all features from a map.

I don't understand why the notation in the first summation where j goes from 1 to n, shouldn't j goes from 1 to however many neighbors that are surrounding i instead of going through all features?


texts: Access to custom GeoJSON members in OpenLayers 2

texts:
I need to read a layer-wide attribute, from a GeoJSON file, set by a custom member to the object:

{"type": "FeatureCollection", "features": [...], "search_polygon": "POLYGON (...)"}

Since it seems neither the Vector Layer nor the Layer classes have any attribute or method to expose this data, I tried to use the callback from HTTP Protocol with no success. (I put an alert to the callback function, but no response.)

My question is: Is there any way to access that search_polygon custom member using OpenLayers 2?


texts: Multiplying multiple specific areas in NetCDF file with constants

texts:
I have wind speed NetCDF files for the UK (land and sea areas) I wanted to multiply a constant (say 1.5) by the wind variable for all land areas and a constant (say 2.4) for all sea areas. How can I do that, preferably in Python? If not in Python, what else could I use, as I have a lot of these .nc files and I will loop over them to apply the constants? An example of a .nc file is in the link.
https://www.dropbox.com/l/scl/AABKha2vlU90gbR25TgzVE9Wkq1kDvn1ntM


texts: Programatically export map clusters in QGIS

texts:
I'm trying to export each microcluster of a large map into a composed map that shows the ff on an A4 sheet:

1) Cluster/grid box of the map

2) Larger index map with the pertinent region highlighted

Workflow for each cluster:

Shade the intended region in canvas
Reflect the shading in the index map (map composer)
Adjust the view of the cluster map (map composer, same page)
(layers are locked for both map items i.e. index map and cluster map)
Save the entire file as PDF

Repeat 50 times. :(

Now, there's a LOT of clusters. I've spent a whole day on a region alone, and there are 5 other regions.

Is there any programatic way to do this?


texts: installing GEOS, PROJ, GDAL/OGR into a python virtualenv on Mac OS X

texts:
I am trying to setup an isolated python virtualenv to work on GIS projects on my Mac OS X. Sounds like I will need to install the GEOS, PROJ, GDAL/OGR from the kyngchaos site here http://www.kyngchaos.com/software/frameworks and not easy_install it into my virtualenv. Is that accurate?

The downloads available on kyngchaos are pkg installers and not just an egg or bdist of python modules. So, I am not very sure as to what else the installer is doing to the environment other than just copying files to /Library/Framework/geos.framework folder.

How would I go about installing a specific version of GEOS, PROJ, GDAL/OGR into my new python virtualenv? I would appreciate your help/suggestions on this issue.


texts: Rasterize points by 'most frequent' in R

texts:
I have a series of spatial data points that I want to convert to a raster using the R function rasterize. The data I want represented is a factor value, so I want to represent the most frequently occuring value for each raster cell.

I have tried to write a funtion to get the mode, but it fails.

getmode <- function(v) {
uniqv <- unique(v)
uniqv[which.max(tabulate(match(v, uniqv)))]
}

run rasterize

test <- rasterize(best, rast, best$val, fun = getmode)

results in this error:

Error in FUN(X[[i]], ...) : unused argument (na.rm = na.rm)


texts: Changing 1000 values in attribute table using ArcPy?

texts:
I work on a point layer which contains 2000 features. I'm trying to change all the values "industry" in field "name" only (it has 1000 features of industry values) into "landUse". The FID of the those features start from FID 500 til FID 1500 in running order.

In this code, continually to Change specific rows in attribute table using Arcpy in ArcGIS 10.3

import arcpy

lyrName = r"G:\PROJECTS\gis\layers\New File Geodatabase.gdb\digum2"
with arcpy.da.UpdateCursor(lyrName, ["name1", "OID@"]) as cursor:
    for row in cursor:
        if row[1] in (70, 71, 72):
           row[0] = 'no land use'
           cursor.updateRow(row)
           print """Updated row {}""".format(row[1])

I can change several rows, but I can't write 1000 number in the brackets.

How can I change those 1000 features values only using Arcpy?


texts: How to generate geoserver ImageMosaic properties not by geoserver?

texts: How can I generate properties for ImageMosaic not by geoserver? I can do index *.shp with gdal as here, but if geoserver not see .properties file it remove this index file and generate it again.


texts: Viewing components of complex ENVI .hdr Labelled Raster

texts: An output of my InSAR analysis is an ENVI .hdr complex raster with data type 6. The data is stored as a real-imaginary pair of single-precision floats. QGIS can load the raster, but displays nan for every pixel. Is there a way to view complex float rasters in QGIS?


texts: Problem with sampling points vectors on r.sim.water Grass GIS

texts:
I'm trying to run r.sim.water on grass GIS. it works but i don’t get any result on the text file i saved on the variable "logfile". I have tried different vector point for sampling locations. Does the sampling location vector point’s map need any specific characteristics?

What kind of info is supposed to be stored on this text file?


texts: Using FeatureClassToFeatureClass_conversion to convert tab to shapefiles?

texts:
I'm trying to use FeatureClassToFeatureClass_conversion() to convert tab files to shapefiles. When I run the code below (no errors reported) it prints the list of all tab files in the workspace like this.

C:/temp1\circle.tab
C:/temp1\square.tab
C:/temp1\triangle.tab

import arcpy
from arcpy import env
import os

path = "C:/temp1"
  for root, dirs, files in os.walk(path,topdown=False):
    for file in files:
      if file.endswith(".TAB"):
         tablist = os.path.join(root, file)
         print(tablist)
         **#Unsure what to put here next to add the tab files above into a feature class/layer**

Since FeatureClassToFeatureClass_conversion() requires an input feature class or layer, how can I add the list of tab files above into a feature class/layer so FeatureClassToFeatureClass_conversion() can run?


texts: How to sync/update an out-of-sync replica?

texts:
We have ArcSDE replicas (both one-way and two-way) of our databases between a number of different servers/locations, and have found that some of the feature classes are out-of-sync even though the replicas have completed successfully.  As a result there are more features in one FC than in its replica, and potentially also attribute values that are also different.

You can see in this screenshot there are 39 more records in one than in the other:

I have tried to manually run the sync (no-change) and also attempted to copy the missing features manually, which resulted in the new copies being replicated (and duplicated) so I had to remove them again.

How can I get my two feature classes back in sync without messing about with features as described above, and without having to recreate the child replica?  Is this even possible?

Windows Server 2012 R2
SQL Server 2012
ArcGIS for Server (SDE) 10.3
ArcGIS for Desktop 10.3.1


texts: Count points by non-continuous polygon in QGIS

texts:
I'm very new to QGIS and I'm learning by practice and reading both QGIS documentation and StackExchange answers. However, I think one of the main tasks has been naming the concepts correctly and look for the correct question, so I want to know which concepts are involved in this problem and how to solve it:

I want to count points by polygon, it may sound easy, the thing is that I don't have well-defined polygons, I have a categorical variable, this is a sneek peak:

I have created a different layer for each categorical variable (that's why each of them have different colors).

So my question is: How do I combine all the polygons of one layer in one non-continuous polygon? The final goal is counting all the points of another layer.

For example if I want to count points in a USA map, Alaska, Hawaii and another non-contiguous territories would count as well.

I would really appreciate if you have another idea to reaching what I want. However, this is how I think I could do it for the moment.


texts: Mars MOLA .csv layer in QGIS misplaced and very small

texts:
I'm trying to add Mars MOLA .csv data on top of CTX and HiRise DEMs but when I load them all up in QGIS the MOLA data is at a completely different location and very small in comparison. I downloaded the MOLA data from Mola PEDR Query and then edited it according to the CTX to MOLA section here
This is how my .csv table looks:


texts: I have a JGW file to add some cropmarks.  I have tried adding it as a vector  or raster file but it comes up with unsupported data each time.  Does anyone recognise this file type and how to load it?


texts: Locking geometry of ArcSDE feature class from ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2.

Does anyone know how I can allow editors to select features and update attributes of an ArcSDE feature class but prevent accidental edits (shifts/deletions) to the geometry?

I have tried making the SHAPE field read-only (via field properties in attribute table) but I was still able to move/delete cells.

I asked here about the "Editable" property for fields, which I had seen in ArcDiagrammer and Enterprise Architect, but you can see Vince's response below.


texts: laspy ValueError: buffer size must be a multiple of element size

texts:
I am trying to learn LASPY and am running into road blocks off the bat. I first ran into and Error 8 (Not enough storage is available to process this command) using my full sized LAS file (1.7GB). I was able to get around it by using FUSION's polyclipdata to create multiple smaller sized LAS files; just frustrating. The second road block I have run into, and I can't find out much about it, is:

ValueError: buffer size must be a multiple of element size

when I try:

import laspy
import numpy as np
from laspy import file
infile = file.File("F:/P834/temp/pointfile3_2.las", mode="r")

Error occurs on line 4.

I am trying to work through the tutorial http://laspy.readthedocs.io/en/latest/tut_part_1.html

What am I missing in my early start?

Full error printout:

Traceback (most recent call last):
  File "F:/P834/laspy/laspy_tests.py", line 2, in <module>
    infile = laspy.file.File('F:/P834/temp/pointfile3_2.las', mode='r')
  File "C:\Python27\ArcGIS10.3\lib\site-packages\laspy\file.py", line 64, in __init__
    self.open()
  File "C:\Python27\ArcGIS10.3\lib\site-packages\laspy\file.py", line 75, in open
    self._reader = base.Reader(self.filename, mode=self._mode)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\laspy\base.py", line 273, in __init__
    self.setup_read_write(vlrs,evlrs, read_only=True)
  File "C:\Python27\ArcGIS10.3\lib\site-packages\laspy\base.py", line 321, in setup_read_write
    self.data_provider.point_map()
  File "C:\Python27\ArcGIS10.3\lib\site-packages\laspy\base.py", line 163, in point_map
    offset=self.manager.header.data_offset)
ValueError: buffer size must be a multiple of element size


texts: QGIS 1.8 tries to use UTM in spite of geo projection file

texts:
I'm wrestling with a problem that I am sure I've created for myself, but cannot resolve.
I received a set of shapefiles that were expressed with UTM, and I want them to use lat/lon. QGIS will do the translation, of course, but there are hundreds of files and I didn't want to do this manually.
So I obtained the code for converting the UTM (N, Zone 14) to lat/lon and re-wrote the .shp files, taking care to change the bounding box values in the main header, the bounding box values in the record headers, and all the point values (these are polygons).
I replaced the .prj file with this, which reflects the constants I used in the conversion:

GEOGCS["GCS_North_American_1983",   DATUM["D_North_American_1983",    SPHEROID["GRS_1980",6378137,298.257222101]],    PRIMEM["Greenwich",0],    UNIT["Degree",0.0174532925199433]]

Yet QGIS treats the file as if it's UTM -- it displays nothing on the screen, and the cursor position is expressed as easting/northing.

If I look at the layer properties, in the Metadata tab it shows this:

Extents:
In layer spatial reference system units : xMin,yMin 677384.09,4068257.57 : xMax,yMax   732884.09,4351757.57
Layer Spatial Reference System:
+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs

Am I correct that there appears to be dissonance between the min/max units and the projection?

I cannot find anything that is telling QGIS to display this using UTM, but obviously,
I'm overlooking something.
GlobalMapper displays the file without complaint.

Any thoughts on what I'm overlooking?

Thanks.


texts: Application of GIS in passive and agent defence

texts: I am working on passive defense modeling, using GIS. Modeling passive defense can be better implemented using GIS and I've come up with my own methods. It would be far better for me if I know there are any previous works on this area. I want to be informed about any application of GIS in passive and agent defense that are publicly available. For example, where to locate the military bases, so that we have the minimum casualties in the case of external attack.


texts: How can I create a choropleth from the 2010 census ZCTA shapefile using R and ggplot2?

texts:
I have created state and county choropleths in R using ggplot2.  The method I have used for this came from Winston Chang's "R Graphics Cookbook".  For those choropleths the map was already represented as a data.frame in R (as opposed to a shapefile).

I now need to create choropleths of the continental US at the ZIP code level.  The boundaries that I would like to use are the US Census ZCTA boundaries which you can download as a shapefile here (select "ZIP Code Tabulation Areas" from the dropdown).  The shapefile is large, at over 500MB zipped.

I do not have a background in GIS, and attempted to create a choropleth from the shapefile in R using ggplot2 by following the advice on the ggplot2 wiki: Plotting Polygon Shapefiles.

However, when I attempt to follow the advice listed there the program just hangs:

require("rgdal") # requires sp, will use proj.4 if installed
require("maptools")
require("ggplot2")
require("plyr")

setwd("~/Desktop/tl_2013_us_zcta510")
zcta = readOGR(dsn=".", layer="tl_2013_us_zcta510")
zcta@data$id = rownames(zcta@data)

zcta.points = fortify(zcta, region="id")

The call to fortify hangs - even running it overnight it did not complete.  I have dug into the code for fortify a bit and can see the line that it causing it to hang (coords <- ldply(model@polygons,fortify)), but am not sure how to fix it.  I have also tried to follow the advice I have seen on numerous other websites about reading in shapefiles into R and ggplot2, but they don't seem to work for me either.  I can't tell if it is because of the size of this shapefile or for some other reason which I do not understand because of my lack of knowledge about GIS.

It would be great if someone could help me with this.


texts: Moving points 5 feet "downstream" along polyline using ArcPy

texts:
I've got a bunch of points that currently sit at the origins of a bunch of different lines. I need to move every single point 5 feet "downstream". I'm using ArcMap without an advanced license so I can't use Generate Near Table or anything in the Proximity toolbox, which seems like it would be helpful.
I'm new to Linear Referencing but it seems like something using Transform Route Events might be what I'm after?
I'd also really like to write a python script to accomplish this for a large number of points, but I'm just stuck on where to start.


texts:
I hear people complain:

  my biggest gripe with google spatial apis is that they hide the fact that latitude/longitude is not a universal coord system

My question is Is latitude / longitude a universal coordinate system? What are the implications if it is not?


texts: Calculate values for all lat/long in another table based on where calculation

texts:
I'm working with data in Postgresql+PostGIS and I've stumped myself (probably from overthinking it). For each lat/long in a site table, I want to calculate an average of a field for all points within a certain distance of that lat/long. For example:

SELECT
       site.id,
       ring1.avg as r1_avg,
       ring5.avg as r5_avg
FROM
       (SELECT id, lon, lat FROM sites) as site,
       (SELECT avg(sales) FROM revenue WHERE ST_DWithin(geom::geography, ST_SetSRID(ST_MakePoint(site.lon, site.lat)::geography, 4326), (1.0*1609.34))) as ring1,
       (SELECT avg(sales) FROM revenue WHERE ST_DWithin(geom::geography, ST_SetSRID(ST_MakePoint(site.lon, site.lat)::geography, 4326), (5.0*1609.34))) as ring5;

Yep, the implicit/explicit syntax is bad form and wrong. That's what I'm struggling with: how to get the subqueries to pickup the lat/long from the site table and use it for the WHERE calculation?


texts: How to join two tables using date, time, lat and lon as common field?


texts: Comparing 2 GPS location data files to find if any location from one file is within a specified elliptical region of a location from another

texts:
I have 2 CSV files with GPS coordinates. I want to check if there are any locations in the first file which are inside a specified elliptical region surrounded by each location in the 2nd file. The major axis , minor axis and major axis bearing(angle) is given in the 2nd file for each location coordinates. For start I have plotted the points from both files in an ArcMap document where I can see red dots as first file locations and green dots as 2nd file locations.

For example if the coordinates of one location in 2nd file is 39.6145,-89.7823 (latitude(h), longitude(k)) and the major axis(a) 1200 and minor axis(b) 400 and bearing(A) 114 degree, then I want the list of the instances in first file which falls inside that region

N.B its is basically solving (x,y) in following equation,

The image is basically the ellipse equation if you can't view it.
〖((x-h)  cosA+(y-k)  sin〖A)〗〗^2/a^2 +〖((x-h)  sinA-(y-k)  cosA)〗^2/b^2 ≤ 1

I am totally new in ArcGIS. I am also trying to do it in Python but want to know how to implement it in ArcGIS too. And both the files are pretty large, with millions of entries in both. So I would also like to know if it is too much time consuming in ArcGIS.


texts:
I'm new to GIS. I'm working on a project right now that is about assisted living facilities in my state. I have a map of my state that I've divided into six regions by creating Thiessen polygons. The Thiessen polygons are based off of state office building locations. I also have a table of assisted living facilities that I have geocoded into a point layer. Each of these points has a population attribute that represents the population of that facility. What I would like to do is divide each Thiessen polygon into two or three regions, such that the sum of the population for the points in a region is as equal to the other regions as possible. I'm using ArcGIS Desktop 10.6.

Is there a tool that I can use for this?


texts: Are ECEF coordinates determined by PPP post-processing valid between different receivers and antennae?

texts:
I was reading the tutorial for PPP using gpsd, and it emphasized how the position and altitude produced by different receivers can produce wildly different results depending on the datum and ellipsoid used.

Am I correct in interpreting this to mean that the raw ECEF coordinate returned by a PPP service will be applicable between receivers and antennae? For example, if I used a ZED-F9P receiver with a ANN-MB antenna to obtain an ECEF position accurate to within a few centimeters, could I then swap that out with a NEO-M8T and an ANN-MS antenna centered on the exact same position, using the previously obtained position as the Fixed Time Mode True Position (ECEF) value on the M8T?

To put it more simply, are the post-processed PPP ECEF coordinates the “true” location of a point?


texts: ArcGIS Split Line at Point Tool splitting at intersections - faulty output?

texts:
The screenshot below shows a 'brown_FC' of brown lines generated by this script which i want to be split at the intersection points with the 'green_FC'.
This works just fine if i take the 'brown_FC' as the script outputs it. However I need to erase a polygon FC off the 'brown_FC' before and if done so the Split Line at Point Tool fed with 'brown_FC_Erase' outputs this:

The most upper split is correct - same as with 'brown_FC' without the Erase before, the others split within the 'brown_FC_Erase' apparently. Why and how to avoid?

What I've tried so far in different combinations:

Using the Integrate Tool as suggested in this post.
Converting 'brown_FC_Erase' from multipart to singlepart.
Tinkering with the search radius which only seems to worsen accuracy of the output as oppesd to when left blank.

The faultiness of the output suggests that this problem isn't just of geometrical but also of systemic nature: Instead with the intersection-points FC  (as defined in the toolbox) the lines get split where they intersect with themselves somewhere near the intersection-points FC (no matter if search radius is set to zero or left blank).

This problem leads me to the conclusion that either the tool might be somewhat corrupted or it's me overseeing some detail here...

Hopefully the latter is the case. I'm on 10.2.2 btw.


texts: Export whole Projekt to multiple geotiff

texts:
I want to export my whole Projekt (with the current activated Layers) to multiple geoTIFF's. Each TIFF images should be have a dimension of 1km x 1km. My project consists polygons out of a shapefile. So in principle like many screenshots with 1km x 1km saved in geoTIFF's.

How does it work?


texts:
I'm using ArcMap 10.5.

I have two individual days of precip data rasters from NOAA/NWS.

I want to combine the two days with the cumulative total storm precip displayed in a new raster.

Many of the cells contain no data, and display as a value of -10000.

I've tried the raster calculator but the missing data cells cause a strange display.


texts: How to fix org.mapfish.print.JsonMissingException: attribute [spec.pages[0].mapTitle] missing?

texts:
I've installed GeoServer's map printing module. I've tested that it works by following the instructions here. Now, I'm using GeoExt.PrintMapPanel, and my code is as shown below:

var mapPanel = new GeoExt.MapPanel({
        renderTo: "mappanel",
        map: map,
        layers: [wms, gLayer, vector],
        center: new OpenLayers.LonLat(21.0936111111, 7.18805555556),
        zoom: 2.5,
        tbar: toolbarItems,
        bbar: [{
            icon: "img/printer.png",
            tooltip: "Print map",
            handler: function() {
                var printDialog = new Ext.Window({
                    autoHeight: true,
                    width: 350,
                    items: [new GeoExt.PrintMapPanel({
                        sourceMap: mapPanel,
                        printProvider: {
                            capabilities: printCapabilities
                        }
                    })],
                    bbar: [{
                        text: "Create PDF",
                        handler: function() {
                            printDialog.items.get(0).print();
                        }
                    }]
                });
                printDialog.show();
            }
        }]
    });

where printCapabilities is http://localhost:8081/geoserver/pdf/info.json?var=printCapabilities, my local installation. I've written this by following the example here.

However, when I try printing, the following error is shown in GeoServer log:
org.mapfish.print.JsonMissingException: attribute [spec.pages[0].mapTitle] missing. Any idea on how I can solve this will be appreciated.


texts: Changing mask value for image in Google Earth Engine

texts:
Is it possible to change the mask of each pixel by calculating mean value across the time?

var IC = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
            .filterBounds(AOI)
            .filterDate('2017-01-01', '2017-12-31')
            .select('B5');

var ChangeMaskAcrossTime = function(image) {
???
};

var output = IC.map(ChangeMaskAcrossTime);


texts: Centroids: Based on travel distance/time and not straight line distances

texts:
I am working on a Sales Person Deployment project and quite often come across requests for calculating centroid of a territory where work is distributed across 1000s of ZIP codes. Increasingly we get requests to calculate centroids based on travel distance and not straight line (or crow fly) distance.

What is the most optimal way to calculating centroid based on travel time (i guess it is similar if not same to point of minimum distance)?

I have attempted to calculate distance (using Google API) between all the ZIP codes with workload and then take the approach of travelling Sales Man problem to come up with a starting point. But there are couple of challenges here:

A sales person may cover multiple targets if s/he goes to a cluster and
Using this method is time and resource intensive

Excuse any digressions from Forum rules as I am a novice here and catching up with the protocols. Looking fwd to be part of this interesting community.


texts: Viewing ST_MULTIPOLYGON data from ArcSDE DB2

texts:
We have inherited some tables that were created with ArcGIS/ArcSDE, and the data in these is still being used by our main web application to perform some spatial lookups.  The queries are pretty simple, we just pass in some x/y coords and the query returns a code uniquely identifying the shape that the coords exist in.

The majority of the mapping solution has now been decommissioned as we are no longer providing live maps, and thus we no longer have the Arc toolset.  The basic lookups remain however, and the customer wants to provide us with a new shapefile for import.  We used to use the tools to check that the new shapefile fit in with the existing ones visually, and to assign the unique code to the shape.  Without the tools I am unsure how to proceed.

There are three tables, each of which is associated with a different lookup, and each table has a column that looks like:

Column     Type      Type
name       schema    name               Length   Scale Nulls
---------- --------- ------------------ -------- ----- ------
SHAPE      DB2GSE    ST_MULTIPOLYGON           0     0 Yes

along with some other columns.  Can anyone recommend any free tools we can use to view this data visually?  Is it even possible?

Sorry if this comes off as very basic or n00bish, I'm a Java & Web developer and have no real experience with GIS solutions.  If you need any other info please let me know.


texts: How to do a Closest Facility Analysis using OSM?

texts:
I have  locations of about 1100 Branches across the country, and I want to find the closest branch for any user-provided point.

I'm looking at OSM data, and supporting tools to do this, since I don't have access to ArcGIS Network Analysis. I was planning to set up an instance of a tool like GraphHopper or something similar.

I looked at the tools mentioned here, and none of them mention Closest Facility solving.

Is there any tool which does closest Facility on OSM data, or am I going to have to calculate the distances from Multiple facilities to my input find, to find the closest tool?


texts: ArcMap: How to change mouse pointer to an arrow in Mapview - VB.NET

texts:
I'm writing a add-in for ArcMap 10.5 in VB.NET. I've chosen the item "Button". So when I click on the button I need the mouse pointer, when in MapView, to turn into an arrow (like the Select Feature arrow). Right now it only turns into the Pan (hand).

This is how my current code looks like:

Public Class btnMatareIdentify
    Inherits ESRI.ArcGIS.Desktop.AddIns.Button

    Public Sub OnClick(ByVal cursorID As Object)

        Dim app As IApplication
        Dim appCursor As IMouseCursor = New MouseCursorClass

        Dim pMxDoc As IMxDocument
        Dim pActiveView As IActiveView

        pMxDoc = My.Document
        pActiveView = pMxDoc.ActiveView

        If TypeOf pActiveView IsNot IPageLayout Then
            appCursor.SetCursor(0)
            myMethod()
        Else
            MsgBox("You need to be in MapMode")
        End If
    End Sub
End Class


texts: How to write features in json file after drawing point on map in openlayers3

texts:
Can i save the features on selecting points into same json file?

`function addInteraction() {
        draw = new ol.interaction.Draw({
          source: floodAlertLayer.getSource(),
          type: typeSelect.value
        });
        map.addInteraction(draw);

    draw.on('drawend',
      function(evt) {
        console.log(evt.feature);
        var parser = new ol.format.GeoJSON();
        var features = floodAlertLayer.getSource().getFeatures();
        var featuresGeoJSON = parser.writeFeatures(features);
        var text = JSON.stringify(featuresGeoJSON);
     //   document.getElementById('info').innerHTML = text;
       // alert(features.get("COUNTY"));
        $.ajax({
            url: 'FloodAlert.json',
            type: 'POST',

            data: text

        }).then(function(response) {
            console.log(response);
        });
      },
    this);
  }`


texts: Explaining different Standard Deviation results from same data in ArcGIS Desktop and MS Excel?

texts:
I have done some interpolation in Geostatistical Analyst in ArcGIS and I got the Standard Deviation (SD) a bit different from the SD that I calculated in Excel Spreadsheet.

What algorithm does each package use? Why are the results different?


texts:
I am new to GIS data.

Our organization uses ArcGIS to map parcels of land used for agriculture. The lowest addressable level is currently a single polygon, representing a parcel. Throughout the year, these parcels may change in size or name.

We are importing shapefile data to extract parcel data, and all information associated with each parcel. We are storing this information in SQL Server, merging it with work data, and tracking changes using a data warehouse framework and slowly changing dimensions, but we have run into a problem I can't seem to solve.

Consider this scenario.

Yesterday we had a five acre parcel named A, and work was performed in said parcel. Today, management changed its name to B, and more work is performed. We track those changes, and show parcel A as being in existence until today, when parcel B sprang to life. When we aggregate acres for parcels in which work was performed, we have a total of 10 acres, when in reality, it's the same five acre parcel.

Currently, the parcel name is the key, and needs to be visible for the users to make sense of their data. They know it was named A, but is now B, and want to see the cost per acre of both A and B. The challenging bit is they want the current acres to be accurately represented as 5; not 10.

From a data perspective, I have no way to know this is the same land.

Is there a way to make the physical land be the key instead of an arbitrary name?

I'm envisioning a data store filled with lat/long coordinates for each square foot (or whatever is the smallest addressable unit), and extracting each addressable unit from the parcel polygon. That way, the parcel name is simply an attribute of each addressable unit, and can change as business needs evolve. With this approach, I can see tracking changes over time being easier, as we would track the land rather than the arbitrary labels.


texts: Ogr2ogr - shape to postgis - is there a way to use and preserve an existing primary key?

texts:
I've been struggling with this for a while, tried a few different ideas inspired by these pages:

http://www.gdal.org/ogr2ogr.html
http://www.gdal.org/ogr/drv_pg.html
http://www.gdal.org/ogr/drv_pg_advanced.html
http://trac.osgeo.org/gdal/wiki/ConfigOptions

But I always end up with ogc_fid as the primary key. Anyone know of a way to get ogr2ogr to use an existing primary key for the fid?

Alternatively, does anyone know how I can get this existing key to be an autoincrement sequence after the ogr2ogr load?

I've tried:

CREATE SEQUENCE my_table_my_pk_seq
    START WITH 1
    INCREMENT BY 1
    NO MINVALUE
    NO MAXVALUE
    CACHE 1;

ALTER TABLE public.my_table_my_pk_seq OWNER TO postgres;
ALTER SEQUENCE my_table_my_pk_seq OWNED BY my_table.my_pk;
ALTER TABLE ONLY my_table ALTER COLUMN my_pk SET DEFAULT nextval('my_table_my_pk_seq'::regclass);
ALTER TABLE my_table ALTER COLUMN my_pk SET NOT NULL;

SELECT pg_catalog.setval('my_table_ogc_fid_seq', NUMBER_OF_FEATURES_IMPORTED, true);`

But I get database errors when trying to insert features.

So, to get it to work with wfs-t, I deleted ogc_fid and made my_pk the primary key


texts: Script for removing unwanted images from imagecollection based on masked values

texts:
I have trouble handling a dataset (ArcticDEM) on google earth engine. The images in the ArcticDEM dataset come already masked out, therefore an imagecollection sometimes contains images which have very limited pixels, where many pixels have the value of "masked". I want to create a script, which would automatically remove an image from an imagecollection, if it contains less than a certain amount of unmasked pixels within the region of interest.
I also tried instead of removing images from imagecollection, to append the "good" images to a new imagecollection, however that didn't work for me either.
Does anyone know how to I should proceed?
The image below is an example of an image containing too many masked pixels.
Script: https://code.earthengine.google.com/b629febdaa1448b2a6847233526348b7


texts: Is it possible to auto rotate images in the photo widget?

texts:
If I use the photo widget in the default QGIS form images are not auto rotated according to their exif rotation.
Is there a way to show rotated images?

The correct way to display the image according the the exif rotation would look like this:


texts: Why is osm road data showing up with lines radiating from a point?

texts:
My openstreetmap data has artifacts that I can't seem to get rid of. It is showing the roads underneath the artifacts so I know that the query is working at least partially.

The data is from geofabrik.de for the D.R./Haiti region, downloaded via the .bz2 file into a PostgreSQL database as per the instructions here: http://wiki.openstreetmap.org/wiki/Osm2pgsql.

I'm displaying it in TileMill via the query: ( SELECT * FROM planet_osm_line
        WHERE "highway" IN ('motorway','trunk')
        ORDER BY z_order asc ) AS route


texts: How to make Geological Cross Sections in Adobe Illustrator, and any suggestions for GIS software?

texts:
I have been tasked with producing a geological cross section in Adobe Illustrator. It is a lengthy process, as I am new to the software. I never learned it in my GIS studies. I know that there is plenty of GIS software out there to produce x sections automatically (we have Petrel for example), but they want it done in Adobe Illustrator as it's not a data driven x section.

Does anyone have any familiarity with doing such a thing? I am having to basically use the cutting tool to divide up a box into the sedimentary layers. Then I fill the cut shapes with the correct colour. Complications arise, however:

I have to show wavy lines between some sedimentary layers to indicate a gap, or layers which aren't shown. I am using the distortion tool to automatically produce the wavy line, which means I need to:

a. copy/paste the original area
b. change the appearance so there is no fill, delete all strokes except the one which will turn into the wavy line, apply distortion
c. copy/paste the original area
d. change the appearance so there is no fill, delete the stroke that will turn into the wavy line, change stroke to black line

This means I have 3 elements whenever there is a part of the cross section showing a wavy line: the area which has a fill but no stroke, the wavy line (derived from the area by deleting the fill and strokes which aren't wavy), the rest of the stroke around the area.

Does anyone know of a faster way to do such a thing?

I wish I could just draw a shape using some sort of line tool and AI would detect when there is a closed off area so I could then select a colour fill.

Another problem is the cutting tool does not make it easy to produce matching curves:

In this image you can see that there is a wrinkle across many sedimentary layers, but my cutting skills in Illustrator mean that the lines do not follow a nice uniform shape. Of course this is the case in the real world as well, although at the moment the image has a bit too much of a 'hand drawn' look to it. Part of the problem is that the cutting tool, usually  helpfully, smooths out the cut line as one goes (digitising a smooth curve with a mouse is impossible!). Sometimes this smoothing has undesired effects. It would be good to know if I can apply an offset to another line when cutting.

Is it right that I should be using the cutting tool as my primary way of setting up a x section, anyway?

I'm thinking a better way to do all of this would be to somehow have a tool that enables me to select vertices across multiple lines in sedimentary layers, then drag them up or down. This would then apply the same amount of curve across the layers.


texts: How to delete a non-empty store in geoserver

texts: I have a workspace in with many stores, while i'm trying to delete the store an error shows up saying that it is unable to delete non-empty store. I tried to delete it after adding shape file the store, but only the shapefile is deleted.


texts: Merge polygon without attribute (in field) with nearest polygon with attribute

texts:
In the image you can see districts in BLUE with attributes (name of district) and small pieces in RED. The red polygons belong to the blue and I want to merge them. The RED are only towards the outside and never surrounded by BLUE.

How can I add/merge/union/eliminate the RED with the nearest/closest/longest-sharing-border of a BLUE.

I am using ArcMAP 10.2 on Win7 and I have tried the Eliminate Tool but that does not give the desired result, because the area of the RED is sometimes larger and would absorb a nearby BLUE = attributes are lost. Also small BLUE districts surrounded by other BLUE are getting merged = attributes and shape are lost.


texts: TypeError: L.Proj is undefined in leaflet

texts:
I am using leaflet 0.7.3. I am trying to load a wms tile:

var crs = new L.Proj.CRS('EPSG:2236', '+proj=tmerc +lat_0=24.33333333333333 +lon_0=-81 +k=0.999941177 +x_0=200000.0001016002 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs',
{
    resolutions: [8192, 4096, 2048],
});

// pass crs into map constructor
var map = new L.map('map',
{
    crs: crs
});

...

But am getting an error:

  TypeError: L.Proj is undefined

Has anybody received this error, and if so, how did you resolve this issue?


texts: Issue merging multiple polygons using ST_UNION

texts:
I've got regional data loaded into a PostGIS table and want to merge regions that share the same parent name (and ultimately output the data as KML.)

This is my query (I'm using php):

SELECT ST_AsTEXT(geom) FROM ST_UNION(ARRAY(SELECT geom FROM RegionsRaw WHERE (ParentName = '" . $parentName . "'))) AS geom;

When I run this script, I'm getting strange results:

It's pretty obvious that certain areas are a little messed up and I have a hunch that it's because some of the polygons I'm merging are next to each other (where I'm seeing success) while others are islands and separated from the rest (getting weird results).

Looking for some guidance from the community on this, thanks!


texts:
I have a SRGB image from which I need to get the NDVI value.

I know how to calculate the NDVI, but it does not work with this SRGB image.


texts:
I am relatively new to GIS and am facing some difficulties trying to summarize/pool data together.

I've attached the file to make it easier. FID_pc11_1 (Refers to postcode identifier) and count refers to the number of each postcode. Basically, there are 2107 objectid_1 (which is the different buffer zones around schools), code represents the different lower super output areas in UK.

I would like to summarize my data, such that it would give me data in the format of; in zone 1 (eg.objectid_1= 1),what is the summated count (of postcodes) for each of the code/lower super output layer (eg:E101002711).