texts: Can I get the bounding box if I have just the shp file?

texts:
Is .shp file, without the dbf, .shx, .prj or others, enough to get the bounding box?

I am trying to get user to upload their shapefiles in a web page. I want to obtain longitude and latitude from the point shapefile, and bounding box (long,lat long,lat) from the polygon shapefile (both the shapefile will contain single shape/feature). I would like to think .shp file is enough to obtain the shape, but when I use fiona, it gives me an error. Is there any package that might help?

The code looks like this:

import fiona
from shapely.geometry import shape

shp_file = 'a_polygon_shapefile.shp'
c = fiona.open(shp_file)

# first record
first_shape = c.next()

box_topY = shape(first_shape['geometry']).bounds[3]
box_bottomY = shape(first_shape['geometry']).bounds[1]
box_rightX = shape(first_shape['geometry']).bounds[2]
box_leftX = shape(first_shape['geometry']).bounds[0]

The error I get is:

Unable to open a_polygon_shapefile.shx or a_polygon_shapefile.SHX.Try --config SHAPE_RESTORE_SHX true to restore or create it


texts:
I created isochrones using Iso_area as polygons algorithm of QGIS. I created the green polygon feeding the algorithm with the speed field of the road network. Then I created the blue polygon using the same start point and the same road network, but I only added a new attribute to the road network regarding the speed and fed the algorithm with this new field. Briefly, two polygons were created using the same start point, the same road network, and the same parameters, and the only difference is that the blue one was created by applying a higher speed limit.
I supposed that their shape should be the same but larger in the polygon that was created considering higher speed limits. But, the results did not prove this. These are the created polygons:

The blue and green polygons were created as I explained above, but why, for example, in the places shown in the orange squares, the shapes are not the same. As it shows in this place, the polygon created with lower speed limits (green) can cover an area that the blue one (with the higher speed) cannot cover!


texts: Can I install two versions of QGIS on the same computer?

texts: I thought a few weeks ago someone asked a question about running multiple versions of qgis on the same computer.  But I can't find it back.


texts:
I'm newbie in the GIS field.
My project is a web mapping application that depends on OpenLayers, Geoserver, and PostGIS.
I have a spatial tables called "features-table", that has a one to many relationship to a non-spatial table called "tags-table", so that every feature has many tags.
I want to update a feature with all its related tags; that means the two tables need to be updated. I'm not able to do that.

I've tried the following:

1- Created a view with the two tables, and used a WFS-T reuqest to update the attributes of that view, but found that geoserver treats views as read only.

2- Created a stored-procedure in PostGIS that updates the two tables, but when calling it from geoserver's "Configure new SQL view" it makes no effect.

Can you help me with how to do that?


texts:
In ArcGIS 10.1, I have a shapefile/layer with labels displayed on a map. Can I modify any of the formatting of the labels using Arcpy? Either directly in code or by importing style from another layer or the style manager would be fine. Specifically I want to add a halo around the text.

I have tried using arcpy.mapping.UpdateLayer with a source layer that has the label style that I want, but the only options are to update everything (including changing the source data, which I don't want) and to change symbology only (which doesn't cover the label style). Any way around this? Or another solution?


texts: How to identify multiple points with multiple buffers

texts:
I have 40,000 house points that is mapped in GIS and another layer with all the amenities around the area. I am trying to find the number points of interests (POI) within a buffer of 500m from each house point.

How do I do that analysis it in ArcGIS so that it is exported into a table format - showing each house ID and the number of POI from a 500m buffer.


texts: Is there a browser/management view in QGIS?

texts:
Is there a browser/management view in QGIS 1.8 so that I can manipulate shapefiles and others (such as geo tiffs, etc.) like in ArcCatalog? Maybe a plugin? Right now I have to use my system's file manager to manage them, but is error prone because I have to select five files for a shapefile, instead of just on in ArcCatalog.

BTW, I found the builtin browser dock, but can't get it to do any manipulation other than drag it into the TOC.

Thanks!


texts: Centering hexagon grid layer on point in MMQGIS?

texts:
I would like to use GPS coordinate of one point and create a hexagonal grid layer around it where it is at the exact middle of a hexagon.

When using MMQGIS it produces a hexagonal grid on the current window and produces something like the attached.

The green point is a layer with one GPS coordinate.

I can get it roughly to the centre by moving the window around but I was wondering if there is a more precise way of doing this?

From the other options in the grid layer window (see attached) custom area seems the most promising.

However I am not sure how to use this to get the green point in the middle?

And it obviously depends on how large a grid I what.

UPDATE:

I have solved this when using a square grid but the Hexagon grid is still giving me trouble - mainly because I do not understand what the X spacing term relates to?

First I will go step by step through a rectangle grid:

In the grid window I change the units to degrees
For this example I set X spacing (Longitude) at 0.0001 and Y spacing (Latitude) at 0.0001 - this sets up the size of each rectangle or in this case square.
Select Custom Area from Extent drop down
My point coordinates that I want to base the grid around are:
S -6.58361, E 38.84678
Say we want 10 grid squares width in every direction from the point. We know  my grid squares are 0.0001 degrees (WxL) thus ten grid squares will be 0.001. Thus the custom area can be set as follows:

Top Y = -6.58361 + 0.001
Bottom Y = -6.58361 - 0.001
Left Y = 38.84678 - 0.001
Right Y = 38.84678 + 0.001

We are not quite there as the middle of a grid is the intersection of four grid squares so we need to offset everything by half a square (see attached for result):

Top Y = -6.58361 + 0.001 + 0.00005
Bottom Y = -6.58361 - 0.001 + 0.00005
Left Y = 38.84678 - 0.001 + 0.00005
Right Y = 38.84678 + 0.001 + 0.00005

For the Hexagon I thought I could do a similar thing however there are a couple of problems. Because to of the shape Y and X spacing are related and if you change one the other changes. If I set Y to 0.0001 the grids come out as the the same hight as the previous example (this makes sense) however at this setting X spacing is 0.0000866 - I have no idea what this relates to? I have tried to measure various things on the Hexagon but nothing matches. The width of the Hexagon is larger the the hight not smaller as indicated by the settings.


texts: How to Change the buffer distance from degrees to meters in QGIS 2.14 Essen?

texts:
I am new to QGIS and the whole field of GIS in general.

My project and all its layers are set to the CRS of EPSG:26913, NAD83/UTM zone 13N.

I am trying to use the geoprocessing tool Buffer(s) on a layer to dissolve any overlapping polygons by a certain buffer distance. I want to do it in meters, but it seems the "buffer distance" is not in meters (from what I've seen/read I think its in degrees), as all my polygons turn into one giant polygon across the whole map.

I understand that the buffer tool always uses the layer units, but I am not very familiar with the different projections/datums and I had thought that NAD83/UTM zone 13N worked with meters. I know the project is set to use meters as its distance of measurement, and everything else I've done so far has worked fine using meters.

I have tried to make sure my layer is set to the correct CRS by using "save as" and selecting the NAD83/UTM zone 13 CRS, saving it with a new name, and loading the new SHP.file to the map, but it doesn't seem to have made a difference.


texts: LonLat values returning from Openlayers 2

texts:
I have developed an app using Openlayers 2 to retrieve the latitude and longitude values based on the clicked position.

The map display the area of Great Britain and Northern Ireland.

Below the code I wrote for that.

var lonlat = self._map.getLonLatFromViewPortPx(evt.xy);

However Lonlat value it return is
      Latitude : 515881.5 and Longitude : 172199

I'm not sure whether these values are correct or not. Because Google map gives a different LonLat values for same area. [link]

The values from google maps are Latitude : 51.438677  Longitude : -0.335115

I want to know how to retrieve the LatLon from openlayers 2 if this is wrong.

Edit : This is the object returning from Openlayers end

initialize
lat
:
256465
lon
:
385906.5
__proto__
:
Object
CLASS_NAME
:
"OpenLayers.LonLat"
add
:
(lon, lat)
arguments
:
null
caller
:
null
length
:
2
name
:
"add"
prototype
:
Object
constructor
:
(lon, lat)
arguments
:
null
caller
:
null
length
:
2
name
:
"add"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1867
[[Scopes]]
:
Scopes[1]
0
:
Global
__proto__
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1867
[[Scopes]]
:
Scopes[1]
clone
:
()
arguments
:
null
caller
:
null
length
:
0
name
:
"clone"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1852
[[Scopes]]
:
Scopes[1]
equals
:
(ll)
arguments
:
null
caller
:
null
length
:
1
name
:
"equals"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1887
[[Scopes]]
:
Scopes[1]
initialize
:
(lon, lat)
arguments
:
null
caller
:
null
fromArray
:
(arr)
fromString
:
(str)
length
:
2
name
:
"initialize"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1813
[[Scopes]]
:
Scopes[1]
lat
:
0
lon
:
0
toShortString
:
()
arguments
:
null
caller
:
null
length
:
0
name
:
"toShortString"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1841
[[Scopes]]
:
Scopes[1]
toString
:
()
arguments
:
null
caller
:
null
length
:
0
name
:
"toString"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1830
[[Scopes]]
:
Scopes[1]
transform
:
(source, dest)
arguments
:
null
caller
:
null
length
:
2
name
:
"transform"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1908
[[Scopes]]
:
Scopes[1]
0
:
Global
wrapDateLine
:
(maxExtent)
arguments
:
null
caller
:
null
length
:
1
name
:
"wrapDateLine"
prototype
:
Object
__proto__
:
()
[[FunctionLocation]]
:
OpenLayers.debug.js:1927
[[Scopes]]
:
Scopes[1]
__proto__
:
Object

Edit 2: These are my grid references
"gridReference": "386170,256495",

Edit 3

projection: "EPSG:27700",
        units: "m",
        maxExtent: {
            left: -521611.7877215,
            bottom: -99060.2018315805,
            right: 1107515.6274664556,
            top: 1234593.14792
        },

        scales: [1.025694, 0.533333, 0.238194, 0.139583, 0.090278, 1.35, 1.17, 1.08, 1.04, 1.02],
        resolutions: [1,0.5,0.25,0.1,0.05,0.025,0.0125,0.00625,0.003125,0.0015625,0.00078125],

        controls: ["navigation", "panZoomBar", "scaleLine", "mousePosition", "attribution", "loadingPanel", "navigationHistory"],
        controlPanel: {
            controls: ["zoomBox", "measure", "getFeature"],
            defaultControl: "panPanel"
        }


texts: Saving extent of satellite imagery in ArcGIS Desktop?

texts: I'm using 50cm imagery for the whole country of South Africa, which make zooming and panning painfully slow. Can I save imagery for just the area I am working with (ie. data frame, city, province)? I'm using ArcMap 10 Editors license.


texts: QGIS2 Plugin Replacement for Join Attributes by Location

texts:
I am using QGIS (ver2.18.22) to join 2 layers (polygons) at intersections to derive some summary statistics into the join layer ("Grid" layer). The "Grid" layer is assembled by using Create Grid, and the "Merged" layer is assembled by Merge vector layers tool using 8 different layers that all have the same polygons and attributes, just different values for the attributes.

In PyQGIS I do the joining part by the following:

      processing.runalg("qgis:joinattributesbylocation", "Grid", "Merged",
        u'intersects', 0, 1, "mean,min,max,median,q1,q3", 1, OUTPUT_PATH)

(I have also done this using the QGIS GUI) I am getting incorrect values for the Q1 and Q3 (Quartiles) statistics. For example, Q1 = Q3, furthermore, Q1 > median. It seems QGIS 3 does this correctly though. My question is 2 parts:

Is there a plugin for QGIS 2.18 that corrects this behavior in the core Join Attributes By Location tool for getting quartiles?
What would be the equivalent PyQGIS code for QGIS 3? The api for 2 and 3 are totally different. For example, processing.runalg, processing.alglist are not found in QGIS3.


texts: Extracting information from KML files in R, beyond coordinates and date/time

texts:
I have GPS-telemetry data from animal movement monitoring that records, besides Longitude-Latitude coordinates and date/time, other information (GPS HDOP, elevation and other important stuff).

I can read them in R (using rgdal::readOGR or sf::st_read functions, for example), but they do not recognize this additional information as separate columns. This additional data is read into a column called description that contains all additional information together, in a format that seems like XML. In each row, I have something like (in an example including only HDOP and elevation):

-21.54328, -49.85085, hdop:  1.0, elev: 898.78

The question is: is there an implemented way of reading kml/kmz files with this additional information already organized as separate columns?
And, if not, what would be the easiest solution?


texts: I would like to extract 'landuse' information from OSM data using postgis

texts:
I am trying to extract OSM land use information for a given bounding box in PostGIS. Using this example I wrote my code:

SELECT *
FROM planet_osm_polygon
WHERE
  planet_osm_polygon.landuse &&
  ST_Transform(ST_MakeEnvelope(1.49464, 49.1963130, 1.628030, 49.287031, 4326), 32647);

Unfortunately it returns me an error:

  ERROR:  transform: couldn't project point (1.49464 49.1963 0):
  latitude or longitude exceeded limits (-14)


texts: Unable to serve .terrain files in cesium sandcastle

texts: I am learning 3d terrain generation using CesiumJS. I generated .terrain files usind Cesium terrain builder and kept them in 'cesium/apps' directory for testing purposes and to avoid CORS issues. Whenever, I try to generate terrain I get error that tile at X:0 Y:0 level 0 as well as tile at X:1 Y:0 level 0 were not found even though I added empty files at the specified locations.


texts: Find maximum theoretical join count (spatial statistics)

texts:
In addition to calculating the joint count statistic, I would like to compare the actual count of joins to the theoretical maximum count of joins. I know how to compute the statistic and find the actual count of joins for a given group (under conditions of two group), but how can I find the theoretical maximum given the known count of members in each group? I would preferably solve this using Python or PostGIS. I am currently using the Pysal module to calculate the joint count statistics. The input data are census tract polygon shapefiles with arbitrary boundaries, making it impossible to simply agglomerate from east to west or north to south, etc.

EDIT: I am asking about the maximum possible count of within-group joins.


texts:
I am working with series of raster data.

I am checking temporal outlier of each pixel by using image before and after, then replace current pixel value (t) by average of image before (t-1) and after (t+1).
I am able to run with test of 2 rows x 2 columns, but with the real data (NA value in background of image). It can be run but very slow. It cost more than one week to complete one image (979x526 pixels)

I am a beginner R user.

Following is the script

# rm(list=ls())
# setwd ("D:\\5Validation\\92NDVIValid\\GapTestNDVI")
# I've commented these out because these could catch people out

Files <- list.files(pattern="*.tif")

library(raster)
library(stats)

#### Need to difine variables in loop
Ndvit       <- as.matrix(raster(Files[2]))                # The second image

Ndvit_1     <- as.matrix(raster(Files[1]))                 # The first image

Ndvit_3     <- as.matrix(raster(Files[3]))                 # The third image

Con1        <- (Ndvit-Ndvit_1) <(-0.01*Ndvit_1)           # Condition 1
Con2        <- (Ndvit-Ndvit_3) <(-0.01*Ndvit_3)           # Condition 2

Con         <-  Con1 & Con2                               # Conbine both Condition 1&2 (T==1, F==0)
summary(Con1)

## Dimension of image in series
NRow       <- nrow(raster(Files[1]))

NCol       <- ncol(raster(Files[1]))

#Create blank Table to merge all output
NdviReli  <- matrix(nrow = 979, ncol = 526)

Out       <- matrix(nrow = 979, ncol = 526)

## Need image before and after to fill,
## So it will process from second image to lenghth of image - 1)

for (t in 2:(length(Files)-1)) {
  for (i in 1:NRow) {
    for (j in 1:NCol) {
        x <- t-1         ## Set variabe t-1
        y <- t+1         ## ## Set variabe t+1

        ## Assigned input image as matrix
        Ndvit       <- as.matrix(raster(Files[t]))
        Ndvit_x     <- as.matrix(raster(Files[x]))
        Ndvit_y     <- as.matrix(raster(Files[y]))

        ## Condition
        Con1     <- (Ndvit-Ndvit_x) <(-0.01*Ndvit_x)
        Con2     <- (Ndvit-Ndvit_y) <(-0.01*Ndvit_y)
        Con      <- Con1 & Con2
        if(is.na(Con[i,j])) {

           NdviReli[i,j] <- NA

          } else if  (Con[i,j] == 1) {
          ## If Condition satify
          NdviReli[i,j] <- (Ndvit_x[i,j] + Ndvit_y[i,j])/2

          } else  {
          ## If Condition satify
          NdviReli[i,j] <- Ndvit[i,j]
          }

          Reli <- print(paste(NdviReli[i,j])) ## Using this if want to check the value one by one
          Out <- raster(NdviReli)                                       ## Assign matrix as raster
          Ex       <- extent(c(487124.6, 618624.6, 1902496, 2147246))   ## Assign coordination, projection base on image
          extent(Out) <- Ex
          projection(Out) <- CRS("+proj=utm +zone=47 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
          writeRaster(Out,filename= paste("Reli",Files[t]),format="GTiff", overwrite=T)
    }
  }
}


texts: Finding distance between positions of same point on two different layers using QGIS?

texts:
I am trying to calculate the difference between the positions of the same point on two different layers. One layer is a shapefile, other layer is a text file, and each point has a unique identifier, Stationval.
Found something similar, using a virtual layer, however my python skills are not good enough to convert this.
Was also looking for the Points2One plugin, with no avail.


texts: QGIS: Creating average of a precipitation multiple raster set?

texts: I have a set of 12 precipitation raster layers downloaded from WorldClim, each one representing monthly precipitation data on a global scale. I've clipped all 12 layers with the area of interest. How can I obtain on QGIS a unique raster representing the average of annual precipitation?


texts:
I would like to expand my knowledge coordinate systems, projections and reprojecting data, with a view of writing a reprojection engine.

Could you provide me with some resources (books, websites, scrolls) on the maths behind coordinate systems and reprojecting data?


texts: Updating point and line data attribute from polygons using PostGIS

texts:
I have a number of point and line tables in postgres in which I want to update an attribute according to which area of a 5km grid they fall into. I use the following postgis query to do this:

UPDATE point_table AS B
SET tile = A.tile_name
FROM 5km_grid_polygons AS A
WHERE ST_Within (B.geom, A.geom);

This worked find for most of my data. However, when a point or line falls exactly on a grid boundary (as in screenshot below), then no attribute is copied over.

This has happened hundreds of times, so manually fixing is not an option. I am wondering if anyone can think of a way to create a rule that assigns one grid square to the point or line that is 'in limbo'?

The rule/query would have to be consistent, for instance, always assigning the attribute from the lower or left grid square to the point or line. It cannot pick one of the adjacent grid squares at random.


texts: Is there any modify feature functions like rotate,resize and drag for vector layer in OpenLayers 3 or 4?

texts:
I am looking for a sample code for rotate and resize vectors in Openlayers 4, but I can not find anything about it. I found this sample code http://dev.openlayers.org/examples/modify-feature.html but its version Openlayers 2.

I've done this so far

   var modify = new ol.interaction.Modify({source: source});
  map.addInteraction(modify);

  var catch;
  var typeSelect = document.getElementById('type');

  function addInteractions() {
    catch = new ol.interaction.Snap({source: source});
    map.addInteraction(catch);
  }

I can't find a function like this

OpenLayers.Control.ModifyFeature.ROTATE;


texts:
I'm trying to run a maximum likelihood classification on a Landsat 8 mosaic to assess land cover, but I keep getting the same error during signature generation:

I never got this error before, and looking into the python scripts listed in the error message (i.e. scpdock.py and utils.py), it seems like it is a matter of ROI size:

However, no matter what I do to the ROIs (enlarge, shrink, change the place I sample from in the image, draw different shapes), I always get the same error. So, I come asking for help in deciphering what's wrong. I don't know if this helps, but I'm using QGIS 3.4.4 with GRASS 7.6.0 on Windows 10 Home (v1607) 64-bit with python 3.6.


texts: Counting number of intersections for N linestrings

texts:
I have a database with D entries of polygons representing buildings. Each entry/row has the following information: id, min_height, max_height, geom (min_height is for buildings that have some kind of overpass). Details aside, I am trying to solve the following problem:
(1) Given N LineStrings (similar to lidar beams) with a shared origin and spanning 360 degrees (for simplicity you can solve the N=2 problem, i.e., two LineStrings in pointing opposite directions, or N=4 for a cross) how can I count the number of LineStrings that intersected a building geometry, hence n/N where n <= N.


texts: How to extract values from .tiff file using R?

texts:
Maybe a rather simple question, but I could not find an answer. I am trying to extract long, lat, value from a .tiff file of a map of Africa with long/lat and soil acidity as the values.

Here is some information:

class       : RasterLayer
dimensions  : 9578, 8990, 86106220  (nrow, ncol, ncell)
resolution  : 1000, 1000  (x, y)
extent      : -4108621, 4881379, -5076086, 4501914  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=5 +lon_0=20 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : C:\Users\M17GPPALT.tif
names       : M17GPPALT
values      : 0, 65535  (min, max)

I want to extract simple long, lat, ph value for each point as an array.


texts: Issues synchronizing QField into PostGIS

texts:
I'm having issues synchronizing my QField project into a PostGIS database.
Nothing is synchronized, and I have the following error in the PostgreSQL logs:

2019-08-29 13:06:34 HST INSTRUCTION :  UPDATE "secosud"."elec_moa_supports" SET "_cree"='2001-10-06',"n_rec_trvx"='',"dphoto_1"='DCIM/20190829180916580.jpg',"dphoto_2"='DCIM/20190829181715420.jpg',"dphoto_3"='DCIM/20190829181757880.jpg',"dphoto_4"='DCIM/20190829181816755.jpg',"dtype_pot"='COMPOSITE',"dinc_pot"='OUI',"ddate"='',"dheure"='',"dutilisateur"='',"trvx_cnfrm"='OUI',"prof_pose"='CORRECT',"haut_pot"='9M',"com_haut"='CONFORME',"resis_pot"='150daN',"com_resis"='CONFORME',"depo_exis"='OUI',"com_deppot"='CONFORME',"mass_renf"='NON',"com_mass"='CONFORME',"dist_pot"='NON CONFORME',"com_dist"='PORTEE VERS TRNSFO 42M',"armem"='CONFORME',"com_armem"='',"nat_cable"='T70',"com_cable"='CONFORME',"depo_cable"='OUI',"com_depcab"='CONFORME',"haut_surpl"='CONFORME',"dphoto_5"='DCIM/20190829181833737.jpg',"com_typpot"='CONFORME',"com_inclin"='LEGEREMENT INCLINE VERS PLUS LONGUE PORTEE',"geom"=st_multi(st_geomfromwkb($1::bytea,3297)) WHERE "gid"=5197

 2019-08-29 13:30:16 HST ERREUR:  syntaxe en entrée invalide pour le type date : «  » au caractère 288

Note: this error disappears if I do step 3.

I have done the following, without success:

closed project by opening new project in qfield restarting tablet
manually added "ddate" and "dheure" fields (date, time
respectively; updated using trigger upon db update)
updating "ddate" and "dheure" manually on one point resulted in all added
points to synchronize (new points not synchronized before this
step), existing points with modified attributes not synchronized. I think the
problem lies here, but why, and how to fix? I have been able to synchronize on
occasionally even with these existing fields.
set fid column to hidden

The date fields I have are all nullable. When I query the number of nulls via QGIS, I get the same number as via pgadmin. However, it seems like the NULL values are being interpreted as an empty string only when synchronizing. Is this an issue with the sqlite packaging? How to preserve NULL values instead of ''?

I believe my issue is related this thread, but I cannot figure out at which level the issue is coming from: QField 1.1.0, QField sync 3.2.1, QGIS 3.8, Postgresql 9.6, PostGIS?


texts: Scheduled/emailed reports using ArcGIS Notebook Server?

texts:
Is it possible to send scheduled/emailed reports via ArcGIS Notebook Server scheduled tasks?

An example of a scheduled report would be:

The report/task runs weekly

Emailed to a list of recipients (via an STMP server)

The report would be emailed out if the resultset/count meets certain conditions (i.e. send the report if the resultset has >0 records).

I only have ArcGIS Enterprise 10.7 right now, so it's not possible for me to test this (notebook server scheduled tasks are only available in 10.8 +).
A related question here: Scheduled reports in ArcGIS Enterprise (emailed; condition-based)


texts: Creating a drop-down menu selection with preset options in Python Script Tool - ArcGIS 10.3

texts:
I'm trying to add a simple drop down menu in my python script tool that always has the same two options to choose from. From the list of parameter types, I haven't found a type that allows this kind of selection.

I've seen many questions answered about using the arcpy.SearchCursor function to search through a field and then use those attributes as the drop down menu options, but I just want to simply hard code two options.

I imagine I'd still use arcpy.GetParameterAsText for setting the menu selection to a variable in my script. But in terms of either hard coding the options into my script, or adding them right in the tool shell(?) I'm lost and haven't been able to find any info on this particular task.


texts: How to extract pixel value counts from a raster in QGIS?

texts:
How does one go about obtaining a summary of the counts of cell values within a raster in QGIS? I know QGIS is doing this, in order to provide the summary statistics and histogram in the properties of the layer. I do not need the summary statistics nor a histogram image, but rather counts of individual values.

It looks like there was at one point an attempt to add a special feature for this (Feature #4321, http://hub.qgis.org/issues/4321).

Am I missing something really obvious?


texts: Counting number of vertices in GeoPandas


texts: Manually editing a LiDAR cloud to remove non ground points

texts:
I've recently been tasked with cleaning up a LiDAR file of a river where I need to remove trees, shrubs and power lines that run right across the river.

I was looking into LAStools but it has a limit of 1 million points before you require a license (the file has ~4.7mil points distance ranging from 0.001 to 0.020). Its a relatively small area, but as you can imagine pretty detailed.

lasview was awesome for this job, getting the views just right and entering edit mode to remove points.  Is there an alternative to lasview for this type of editing? What other kind of software are available?

I've already contacted Martin Isenburg about the possibility of purchasing a lasview license but I understand that he is a pretty busy guy and getting a reply by mail might take some time that I don't necessarily have.


texts: Solving routes to offset points in ArcGIS Network Analyst?

texts:
The scenario is that there are points on a grid with intervals of 1000m covering an area of hundreds of square kilometres. I want to find routes on a road network to any one of these points. A road network built with Network Analyst works fine and the points are also built into the network. However the points do not intersect with the roads and sometimes offset 10s to 1000s of metres so no solution/route can be found to these points.

How can routes be solved to the closest part of a road nearest to these points if they are loaded as stops?


texts: Split No Longer Working When Creating Enclosed Polygons

texts:
I have been using ArcGIS Pro to create a large polygon feature class.  Some of the polygons have enclosed areas I need to remove.  I have been using the Split tool under Modify Features and it has been working for weeks. I updated Pro yesterday to 1.2 and didn't have any problems.  Today when I try to use Split to run the same process I have been for weeks I get a message saying

  Split failed. New geometry must intersect polylines and intersect
  polygons in two or more places.

The Split area I draw is completely enclosed so I am not sure why it is not working.


texts:
I am doing some raster subtraction to get a canopy height model and I am struggling to write the output as a raster. The original datasets were .tif files in ESPG 27700 created from LiDAR. I used masked arrays as I know there are some no data values (-3.048...) in the original tif files. Ideally I would like to return the result of the raster calculation as a geotif with the same crs. The code I have used is as follows...
grnd_DTM = rasterio.open(ground_raster_full_path) #path variable to original bare earth tif
ground_array = grnd_DTM.read(1, masked = True) # reading bare earth raster as a masked numpy array
sfc_DTM = rasterio.open(surface_model_full_path) #path variable to surface model tif
surface_array= sfc_DTM.read(1, masked = True) #reading surface model raster as masked numpy array
canopy_ht = surface_array-ground_array #

profile = ground_raster_full_path.profile
band_number = 1
# to update the dtype
profile.update(dtype=canopy_ht.dtype)

with rasterio.open(raster_folder+'canopy_ht.tif', 'w', **profile) as dst:
    dst.write(target, band_number)

Unfortunately this is giving me the following error
AttributeError                            Traceback (most recent call last)
<ipython-input-71-4914f600fe68> in <module>
      5 canopy_ht = surface_array-ground_array
      6
----> 7 profile = ground_raster_full_path.profile
      8 band_number = 1
      9 # to update the dtype

AttributeError: 'str' object has no attribute 'profile'

I'm very new to all of this.


texts:
I've been struggling for a while now with pseudo-nodes.

Is there a tool that removes them automatically?

I have 1400 pseudo-nodes to correct and I'm wondering if there is a way to automate the process. I searched and found something about sextante plugin but that one is deprecated.
I tried to dissolve everything and then snap them at intersections with v.clean but topology check still gives me 1400 errors of pseudo nodes after dissolving.


texts: I'm working on a project which uses GPS coordinates. I have a UAV that flies, for example, over a meridian or a parallel. I suppose that the coordinates I'm reading from UAV are always increasing or always decreasing, but it isn't like that. In details, I can have for example 3 decreasing coordinates and 1 increasing coordinate, then (always for example) 10 decreasing coordinate and 1 or 2 increasing coordinates. Is that normal? It could be the precision of the sensor?


texts:
I have created a simple script tool which has a single parameter of type "String".  I have written the tool validation to disable the input field (simply trying to learn tool validation).  The validator code is as follows:

class ToolValidator:

  def __init__(self):
    import arcpy
    self.params = arcpy.GetParameterInfo()

  def initializeParameters(self):
    self.params[0].Enabled = 0
    return

  def updateParameters(self):
    return

  def updateMessages(self):
    return

Unfortunately, when I open the tool, the field is not disabled (as I expect it should be). Any ideas what I'm missing?

SOLUTION: In case it's not clear in the threads below.  The solution to the issue was to change Enabled to enabled, which is a case change between ArcGIS 9.3 and 10.0.  This capitalization works correctly both in the initializeParameters method and the updateMessages method.


texts: OpenLayers TextStyle background fill only uses one color at a time when rendering

texts:
For my vector tile layer, I have several Point geometry features with text. I'm setting the text style's background fill property for these features to different colors based on the text:

if(label === "blue")
{
    textStyle.setBackgroundFill(new Fill({ color: "blue" }));
}
else
{
    textStyle.setBackgroundFill(new Fill({ color: "red" }));
}

This works when only a single type of label is rendered and, therefore, only one background fill color is used:

However, whenever I have different types of labels (meaning multiple colors for the background fills) it only applies one of the colours to the all the features:

Notice that the labels that were once red (and should still be red) are now blue. The weird thing is, this will change if I pan the map a little:

The styles for the features are set when the feature is created in the tileLoadFunction for the custom source I'm using:

tileLoadFunction = (tile, url) =>
{
    //retrieve the labels using XHR request
    ...

    var labelPoint = new Point(Projection.fromLonLat([labels[i].x, labels[i].y]));
    var label = new Feature();

    //set the text style for this feature
    label.setStyle(_this.getLabelVectorFeatureStyle(tile.getTileCoord()[0], labels[i].type, labels[i].label, labels[i].angle));

    label.setGeometry(labelPoint);

    features.push(label);
    ...
}


texts:
Does anyone know where I can find a sample shapefile that has null entries in it (preferably a polygon shapefile)?

I'm not sure how to create one with null entries.


texts: Simplifying adjacent polygons on subset of vertices only

texts:
We have many adjacent polygons which are very complex and we need to simplify only some portions of them.

There are no gaps nor overlaps between any two adjacent polygons.

What we'd like to do is to select a subset of the vertices delimiting two adjacent polygons together, and then simplify only this subset, while maintaining no gaps or overlaps between two adjacent resulting polygons:

Is it possible to achieve that in QGIS for instance (or any other tool)?

We've looked into some options to simplify adjacent polygons, but can't find a way to do that on only a subset of the polygons.


texts:
I would like to re-calculate large geotiff files (10,000 x 10,000px) with the help of gdal2tiles into little (256 x 256px) tiles.  The process works fine. However, the result is not perfect.  There is a slight shift in the tiles (see the picture below).  I've already tried different resampling methods and get the same result every time.

Here is my statement:

  Gdal2tiles.py -s EPSG: 25832 -z 11-20 /home/input/file.tif output_path /

I think there might be a problem with transforming EPSG:25832 to EPSG:3857, but if I try to transform with gdalwarp I get the same problem...

Does anyone have any ideas?


texts: Choosing Field in Python Script Tool using ArcGIS Pro

texts:
I am writing a script to quickly split large FCs and make them shareable. So I need the tool needs to know what the FC name is and What the Field Name is to split the FC into manageable segments.

To make the Field Name a dropdown would I need to have the user choose a FC and run the tool then have a second tool pop up with the Field Names available?


texts: Any reason to use TileMill if Mapbox Studio is meant to replace it?

texts:
I'm new to using TileMill and it still seems like a popular cross-platform package, though I keep running into common bugs and errors with it (on Mac OSX). On Mapbox's page, it says that Mapbox Studio is meant to "fully replace and improve upon TileMill in functionality".

When learning Mapbox, the service still refers to using TileMill in its tutorials. So am I right in assuming that Mapbox Studio is the way going forward, it's just that they haven't updated their documentation yet?


texts: How to migrate data directory of opengeo geoserver?

texts: How to migrate data directory of geoserver in  opengeosuite 3.1 to stand alone geoserver version 2.5.1. Is it to copy data directory and restart the server only or any thing else..


texts:
I want to do vegetation, water, and soil analysis to an area using satellite imagery. My area of interest is East Java, Indonesia. Here's the map of East Java, Indonesia.

I got 6 scenes of LANDSAT 8 OLI/TIRS that cover the AOI. I want to combine  all the scenes of all bands (I mosaic Band 1 from 6 scenes, and so on) will shorter my time for analysis because I can analysis all the scenes at once. I build a ModelBuilder for combining the six scenes, band by band. Here's my ModelBuilder.

The tools I use in the model is Copy Raster to remove the 'black area' (if you know what is it called, please tell) in the georeferenced scene so it doesn't overlap when I combine all the scenes. I set the pixel type 16-bit unsigned, Ignore Value Data 0, and No Data Value 0.

Then I use Mosaic To New Raster tool to combine all the scenes. I use itu for combining scene in band 1, band 2, band 3, and so on. The model works for band 1 until band 11. But when I use the Model Builder for quality assessment band (BQA), it didn't work.

This is mosaic of six scenes in band number 3 (also my expectation for BQA too), using the ModelBuilder.

This is the mosaic of the scenes in BQA. The black area isn't removed even when I did the Copy Raster without the ModelBuilder.

Is it possible to use Copy Raster tool to BQA to remove the 'black area' and then combining all the scenes?

I only find in the internet that Mosaic To New Raster Tool is used for combining composite bands, I don't know if it's okay to combine single band from all scenes, combine to make composite bands, then use it for analysis. I thought is it okay until I find this problem.


texts: Removing no data values from multiple raster files in PyQGIS?

texts:
I have 180 raster files, which have values in the range of 100 - 900. However, there are no data values in each file designated as 32766. Sometimes, it varies around 32700. I want to write a code that says remove the values that are greater than 30000 and save the file. With the help of GIS StackExchange, I was able to figure out how to save files. But, I want to know how to remove values for all raster files at once.
myDir = 'E:/'
layers = [layer for layer in QgsProject.instance().mapLayers().values()]
pipe = QgsRasterPipe()
for layer in layers:
   extent = layer.extent()
   width, height = layer.width(), layer.height()
   renderer = layer.renderer()
   provider=layer.dataProvider()
   crs = layer.crs().toWkt()
   pipe.set(provider.clone())
   pipe.set(renderer.clone())
   file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
   file_writer.writeRaster(pipe,
                           width,
                           height,
                           extent,
                           layer.crs())```


texts: ArcMap crashes when trying to merge training samples during supervised classification?

texts:
I have a recurrent problem. When I select training samples from a subset of a landsat mosaic, arcmap crashes ALWAYS! Specifically it crashes when i try to merge training samples, sometimes when i try to visualize the scatter plot. My training samples dont have so many pixels.

Some background:

I did a landsat multi-band mosaic of the whole coast of Peru (14 landsat scenes).
I then kept only the coast by clipping with a mask that covered the first 5 km of the coast.
I divided the coast strip in three parts (North, center, south) and finally in each part I identified and clipped all the wetlands (they are not so many). This last multiband tiff containing patches of wetlands is what i have been using to perform the supervised classification.

I tried ArcGIS Desktop 10.3 in two computers and the same happened, I reinstalled the Arcmap, i optimized my computer and I also re-did the mosaic. I am not doing very complicated procedures and arcmap always crashes.


texts:
I think the issue outlined in this question PostgreSQL default values in QGIS 2.8.1 and ticket #12391 is not fully solved yet. Negative numeric default values, and dates are not used by QGIS 2.8.2 correctly. When I work on them and save the layer, the default values were not computed, and the NOT NULL constraints block me. 2.6 works.

Example:

CREATE TABLE bla
(
  vres double precision NOT NULL DEFAULT (-32768.0),
  startvalid date NOT NULL DEFAULT '2001-01-01'::date
)


texts:
I'm trying to load a relatively big Geopackage into my database with ogr2ogr and Windows Subsystem for Linux:

ogr2ogr -f PostgreSQL "PG:user=postgres password=#### dbname=db host=localhost port=5432" de_grid_etrs89_utrm32_100m.gpkg -lco LAUNDER=NO

This syntax works for other gpkg's but not for this one. Error Message:

  ERROR 1: CREATE TABLE "public"."de_grid_utm32_100m" ( "rowid" SERIAL,
  PRIMARY KEY ("rowid"), "featuretype_name" VARCHAR(200), "dataset_name"
  VARCHAR(200), "OBJECTID" INTEGER, "id" VARCHAR(254), "x_sw" FLOAT8,
  "y_sw" FLOAT8, "x_mp" FLOAT8, "y_mp" FLOAT8, "f_staat" FLOAT8,
  "f_land" FLOAT8, "f_wasser" FLOAT8, "p_staat" FLOAT8, "p_land" FLOAT8,
  "p_wasser" FLOAT8, "Shape_Length" FLOAT8, "Shape_Area" FLOAT8, "ags"
  VARCHAR(50), "geom" geometry(GEOMETRYZM,1) )
  FEHLER:  aktuelle Transaktion wurde abgebrochen, Befehle werden bis
  zum Ende der Transaktion ignoriert
  ERROR 1: Unable to write feature 1 from layer de_grid_utm32_100m.
  ERROR 1: Terminating translation prematurely after failed
  translation of layer de_grid_utm32_100m (use -skipfailures to skip
  errors)

My guess is that it's some encoding issue with the german language. Does anyone have a clue on what I should try next?

Gdal version: GDAL 2.4.2
PostgreSQL version: 12

Edit:
Ogrinfo screenshot


texts: Rasterize (vector to raster) errors in QGIS

texts:
I have a point shapefile that I want to convert to a raster in QGIS, where I want to have the values from the points into pixels with a 30m x 30m resolution. But when I run the command an error happens and the following message shows up:
2020-11-28T12:14:10     INFO    gdal_rasterize -l None -burn 0.0 -ts 1.0 1.0 -a_nodata 0.0 -te -60.16838598138009 -4.000042531481377 -57.62922248328133 -2.0386070269701575 -ot Float32 -of GTiff file:///Volumes/New%20Volume/Thesis/Era5/Clipped_area.csv?type=csv&delimiter=,%5Ct&detectTypes=yes&xField=longitude&yField=latitude&zField=Temp_C&crs=EPSG:4326&spatialIndex=no&subsetIndex=no&watchFile=no "/Volumes/New Volume/Thesis/Era5/testing_rasterize.tif"
2020-11-28T12:14:10     INFO    GDAL execution console output
             /bin/sh: /Volumes/New Volume/Thesis/Era5/testing_rasterize.tif: No such file or directory

 FAILURE: No output file specified.

 Usage: gdal_rasterize [-b band]* [-i] [-at]

  {[-burn value]* | [-a attribute_name] | [-3d]} [-add]

  [-l layername]* [-where expression] [-sql select_statement]

  [-dialect dialect] [-of format] [-a_srs srs_def] [-to "NAME=VALUE"]*

  [-co "NAME=VALUE"]* [-a_nodata value] [-init value]*

  [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]

  [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

  CInt16/CInt32/CFloat32/CFloat64}] [-optim {[AUTO]/VECTOR/RASTER}] [-q]

I'm not sure what they mean by NO OUTPUT FILE SPECIFIED.
Then I try changing the resolution of pixels, and making the Layer extent larger, but the following shows up:
2020-11-28T12:24:21     INFO    gdal_rasterize -l None -burn 0.0 -ts 0.1 0.1 -a_nodata 0.0 -te -6835376.729709029 -482088.2328243251 -6308769.931712861 -143164.8265015317 -ot Float32 -of GTiff file:///Volumes/New%20Volume/Thesis/Era5/Clipped_area.csv?type=csv&delimiter=,%5Ct&detectTypes=yes&xField=longitude&yField=latitude&zField=Temp_C&crs=EPSG:4326&spatialIndex=no&subsetIndex=no&watchFile=no "/Volumes/New Volume/Thesis/Era5/testing_rasterize.tif"
2020-11-28T12:24:21     INFO    GDAL execution console output
             /bin/sh: /Volumes/New Volume/Thesis/Era5/testing_rasterize.tif: No such file or directory

 ERROR 1: Wrong value for -outsize parameter.

 Usage: gdal_rasterize [-b band]* [-i] [-at]

  {[-burn value]* | [-a attribute_name] | [-3d]} [-add]

  [-l layername]* [-where expression] [-sql select_statement]

  [-dialect dialect] [-of format] [-a_srs srs_def] [-to "NAME=VALUE"]*

  [-co "NAME=VALUE"]* [-a_nodata value] [-init value]*

  [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]

  [-ot {Byte/Int16/UInt16/UInt32/Int32/Float32/Float64/

  CInt16/CInt32/CFloat32/CFloat64}] [-optim {[AUTO]/VECTOR/RASTER}] [-q]


texts:
I am working in an ArcMap session, I've checked the option to get an ID for each of the layers in my document, so anyone has it.

Now, what I need is a way to get these values within a python script (to copy that automatically in a table): there's no such information in the layer's properties reached by a describe, nor from the .mapping module.


texts:
How can I get the average width of a region or inside a curve?
or
How can I calculate minimum,maximum, and average width of an irregular polygon?
I have attached a figure, I want to know
how can I calculate width of yellow areas which are inside the polygons that exist in the main shape file.


texts: What datum (reference ellipsoid) does Google Earth use?

texts:
What datum (reference ellipsoid) does Google Earth[1] use?

I can't seem to find where Google Earth tells explicitly what datum they use. I need some sort of valid link.

[1] The default free version that you download from Google.


texts: Proper configuration of proxy for deegree

texts:
I am successfully running deegree 3.4.3 in
Apache Tomcat/9.0.14 on windows 10.

However, the imports are not working. It is unclear what I should be doing to setup the proxy info. I run a localhost proxy which sets up all the information for the corporate host. Applications that look for http_proxy environment variable can use: http://localhost:3128  (There is a local service at port 3128 that adds all the credentials and forwards http/s requests to the company proxy).

Here is what the default is, but using that will allow me to connect to repo.deegree.org

  <ProxyHost>http://proxy.deegree.org</ProxyHost>
  <ProxyPort>8000</ProxyPort>
  <ProxyUser>deegree</ProxyUser>
  <ProxyPassword>pwd</ProxyPassword>
  <NonProxyHosts>127.0.0.1|localhost</NonProxyHosts>

I have tried techniques that usually work for other programs here:

  <ProxyHost>http://localhost</ProxyHost>
  <ProxyPort>3128</ProxyPort>
  <ProxyUser></ProxyUser>
  <ProxyPassword></ProxyPassword>
  <NonProxyHosts>127.0.0.1|localhost</NonProxyHosts>


texts:
I'm pretty new to Quantum GIS, and I was looking to add zeros to a field I was having some trouble with.  I learned through browsing a few posts here that existing fields can be modified by defining actions.  My question is, how does this process work?  How do I go about defining actions and what is the language (for example here) used to define these actions?

Example: substr( ( '000' || tostring ("myfield")) , ( length ( ( '000' || tostring("myfield")) ) -2) , 3 )


texts: How to remove points from Google Streets in QGIS?

texts:
I'm working on a map that has already some points directly from the map.

For example: 'Birdland Golf and Country Club'

Is there any solution to remove it from my map? I only need roads, towns and their names.


texts: geomwithin is inconsistently returning features from another layer

texts:
I have two layers, one with more specificity applied than the other (territorial subdivision within subdivision). There are unique keys fields that are equal between both layers. For any feature, it returns null or the correct value consistently for the same features. I rechecked that the key fields are exactly the same strings.
case
  when geomwithin('dom_admbnda_adm2_2020-iso', 'adm2_pcode') = adm2_pcode
    then geomwithin('dom_admbnda_adm2_2020-iso', 'ISO3306-2')
  else -1
end

Maybe I'm misunderstanding something, but isn't geomwithin supposed to behave like a sql join?


texts: I need to display a line-type layer differently according to scale. Do I necessarily have to use a rule-based representation and make several different rules (one per scale interval) or is it possible to achieve that within a single symbol representation? I guess the answer is the first one, but just in case it could be done easier...


texts: Problem with ST_SummaryStats, ST_Buffer, and casting between geometry/geography types

texts:
I have a test raster that covers the world ((-180 90), (180 -90)). The SRID is WGS 84 (EPSG:4326). Values for cells in the north-western quadrant ((-180 90), (0 0)) are random integers in the range 1 to 10 (inclusive). Values for other cells are -99; this is the explicit NODATA value.

The test raster is stored in a PostGIS database (v.2.1). I'd like to compute summary statistics for a buffered point where the buffer is specified in either decimal degrees or in meters.

The Problem: My query using decimal degrees and the geometry type works as expected. However, the same query, using meters and the geography type, which I then cast back to the geometry type, does not. Can anyone suggest why this is the case?

Here's my query, where the buffer is specified in decimal degrees:

SELECT
    filename,
    (stats).*
FROM (
    SELECT
        filename,
        ST_SummaryStats(
            ST_Clip(
                rast,
                1,
                ST_Buffer(
                    ST_SetSRID(
                        ST_Point(-90, 45),
                        4326
                    ),
                    1
                ),
                TRUE
            ),
            1,
            TRUE
        ) AS stats
    FROM
        rasters
    WHERE
        filename IN ('raster_grid_1_1_nodata.tif')
) AS stats_tbl;

As you can see, there are filename and rast columns that contain the name of the input file and the input raster, respectively. The result is as expected.

Here's my query, where the buffer is specified in meters:

SELECT
    (summarystats).*
FROM (
    SELECT
        ST_SummaryStats(
            ST_Clip(
                rast,
                1,
                CAST(
                    ST_Buffer(
                        CAST(
                            ST_SetSRID(
                                ST_Point(-90, 45),
                                4326
                            ) AS geography
                        ),
                        1000
                    ) AS geometry
                ),
                TRUE
            ),
            1,
            TRUE
        ) AS summarystats
    FROM
        rasters
    WHERE
        filename IN ('raster_grid_1_1_nodata.tif')
) AS summarystats_tbl;

As you can see, I pass a geography to ST_Buffer so that the radius is specified in meters. The result, however, is not as expected: PostGIS reports that all pixels have the NODATA value and, consequently, does not compute summary statistics. Can anyone suggest why this is the case?


texts: Using a Raster image in QGIS report linked to a field automatically switches back to SVG and doesn't display the image

texts:
I have a layer setup in QGIS which stores all the information from field collection including images.  So I have image 1 to image 15.  In the QGIS Report composer you select a layer then it will print a report for all the items in that layer.  I want this report to display images that the file paths are stored in the fields.

Soon as you click that it defaults to a SVG image not a raster and give an error saying it can't display the image and you cannot select raster.  I'm using QGIS 3.16.4 and the images are .jpg.
Any suggestions on how to make this work?


texts: How can I add markers and attributes contained in JSON data to an OpenLayers 3 vector layer?

texts:
My question is similar to a previous question I asked (Adding attributes contained in JSON data to OpenLayers popup?).  However, since I originally asked this question, I have upgraded to OpenLayers 3 (OL3) and am having some difficulties.

I am currently displaying an OL3 map using the following snippet code:

    var sourceBingMaps = new ol.source.BingMaps({
    key: '<My API Key>',
    imagerySet: 'Road',
    culture: 'en-us'
   });
    // Microsoft Bing tile source for roads
   var bingMapsRoad = new ol.layer.Tile({
    preload: Infinity,
    source: sourceBingMaps
   });
    var bingMapsAerial = new ol.layer.Tile({
    preload: Infinity,
    source: new ol.source.BingMaps({
    key: '<My API Key>',
    imagerySet: 'AerialWithLabels',
    })
   });

   var map = new ol.Map({
    renderer: 'canvas',
    layers: [bingMapsRoad,bingMapsAerial],
    controls: [],
    target: 'map',
    view: new ol.View({
    center: ol.proj.transform([0,0], 'EPSG:4326', 'EPSG:3857'),
    zoom: 3,
    minZoom: 3,
    maxZoom: 19
    })
   });

The code displays the BING layer without issue.  I have a JSON file, with a sample of the data below:

"LOCATION_A": {
      "latitude": 10.2070,
      "longitude": 25.3215,
      "altitude": 1.2,
      "status": "CLOSED",
      "serverAddress": "192.168.0.1"
 }

I am using the code below with jQuery to cycle through the JSON file and place the markers on the BING map layer.

$.getJSON('json/XMIT.json', function(data) {
      $.each(data.markers, function(key, val) {
        var pointFeatures = [];
        var markerName = this.name;
        var markerLocation = this.location;
        var markerLatitude = this.latitude;
        var markerLongitude = this.longitude;
        var markerAltitude = this.altitude;
        var markerStatus = this.commsState;
        var pointGeometry = new ol.geom.Point(ol.proj.transform([markerLongitude, markerLatitude], 'EPSG:4326', 'EPSG:3857'));
        var pointAttributes = {'location': markerLocation, 'latitude': markerLatitude,
                               'longitude': markerLongitude, 'altitude': markerAltitude,
                              };

        var pointFeature = new ol.Feature(pointGeometry, pointAttributes, {
            title: markerLocation,
            pointRadius: 16,
            });

        var vectorSource = new.ol.source.Vector({
            projection: 'EPSG:4326'
        });
        vectorSource.addFeatures([pointFeature]);

        var vectorLayer = new ol.layer.Vector({
            source: vectorSource
        });
       });
      });

When this code is executed, the browser console provides the following error:

expected expression, got '.'

In regards to this variable:

var vectorSource = new.ol.source.Vector({
        projection: 'EPSG:4326'
        });


texts: QGIS and PYQT5: Getting QgsAttributeForm.modeChanged signal

texts:
I wonder how I can catch the signal, that attributes in a form are editable?

E.G. I open a form to see the values of the attributes of my geodata. Then I click on the "pencil" in QGIS (to make data editable). And right now a message should occure.

I tried:

def test(): qgis.utils.iface.messageBar().pushMessage("Status: ", u"Status: ")

    myDialog.modeChanged.connect(test)

...but that`s obviously not the correct way. Could anybody help?


texts:
In my python plugin ,i select features from a vector layer dynamically using a loop and executing :
    VectorLayer.select(feature.id())

when i start the loop, the selection become very slow because the rendering .
can i suspend and reactivate  rendering in qgis programatically in python.


texts: Deleted layers still appear while loading project in QGIS?

texts:
I'm working on a project where I digitized a few maps (georeferenced) and I also added a couple of csv files. I deleted the raster and cvs layers after I included them into the project as shapefiles but when I re-load the project to continue working on it these phantom layers still appear in the layer panel until the project is actually loaded and ready to work.
It doesn't make any problem other than these layers being listed as options to choose, for instance doing interpolation or other spatial options. The problem occurs when I delete the original csv file from the original location, even though it is not visible in the layer panel. If I do that, QGIS closes with error immediately after opening. I cannot move the csv or raster files from the folder they were originally saved because I cannot load the project.

Same thing was discussed in qgis.org but the discussion was closed apparently because it was solved, which is not my case. Any clues?

Windows Vista - QGIS Lisboa 1.8.0


texts: Bit of a theoretical one this, but I was wondering if anyone had tried to use CartoDB and Torque with any big(ish) data. What I am talking about is having the location of 75,000 people at 1 minute resolutions for a 24 hour period.  Would torque 'work' or would it just crash? I've done it successfully with 200 people, but am now getting to the stage where I'll want to test this with much larger datasets.


texts: Using Reclass in Arcpy, Error 010005: Unable to Allocate Memory

texts:
I've derived a raster of slopes from a 10 m DEM resampled to 15 m for an area of approximately 20,000 mi^2. I want to use Python to reclassify those slopes into categories according to range, but whenever I run the code (and I have tried on 3 different machines) I get the following error:

  "ExecuteError: ERROR 010005: Unable to allocate memory. Failed to
  execute (Reclassify)."

My understanding based on this description (http://pro.arcgis.com/en/pro-app/tool-reference/tool-errors-and-warnings/010001-020000/tool-errors-and-warnings-10001-10025-010005.htm) is that the machine I'm working on does not have enough system memory. However, considering that the area is not very large and I have tried it on multiple machines, this doesn't make sense to me. Working on a 16 GB machine with 15.9 GB available. Am I missing something/does anyone have suggestions for why my code is not working? The code is listed below

import os
import sys
import arcpy

arcpy.env.workspace = "D:\\SanDiskSecureAccess\\ProgrammingProject\\mudslideanalysisdata.gdb"

slope = "Slope_DEM_ma1"

arcpy.env.mask = "DEM_mask"
arcpy.env.cellSize = 15

slopeRemapRange= RemapRange([[0, 4.5, 1], [4.5, 9, 2], [9, 15, 3], [15, 25, 4], [25, 35, 5], [35, 80.816727, 6]])
outReclassSlope = Reclassify(slope, "VALUE", slopeRemapRange)
outReclassSlope.save("SlopeReclass")


texts: What is the equal of ArcGIS Raster calculator "Mod" function (Finds the remainder of the first raster when divided by the second raster on a cell-by-cell basis.) in r. mapcalc in GRASS GIS?


texts: Extracting (x,y) of cells overlapped with polyline in ArcGIS for Desktop

texts:
I have a grid which is made up of 100*100 cells. I want each cell to be identified by its x (1 to 100) and y (1 to 100) value . For example the cell on the top left corner must be Identified as (1,1) and the cell on the bottom right corner be identified as  cell (100,100). There is a polyline which overlap some of these grid cells. I want to extract the (x,y) of each overlapped cell.

How can I do it in ArcMap?(Without python)


texts:
I have a system that uses shortest_path function pgrouting 1.x. At the moment, I want to migrate to pgrouting 2. And I assume the equivalent function to shortest_path is pgr_dijkstra. (Correct?)

When running this, (Hibernate Query Language)

SELECT SUM(t2.length) as totalLength FROM (
    SELECT * FROM pgr_dijkstra(
        'SELECT gid AS id, source::int4, target::int4,
                cost::double precision AS cost,
                reverse_cost::double precision as reverse_cost
            FROM openstreetmap',5374,5364, false, true
    )) AS t1, openstreetmap AS t2 WHERE t1.edge_id = t2.gid

I get error
ERROR: column t1.edge_id does not exist

What is the edge_id in the output of pgr_dijkstra ?


texts: Using the field calculator in rule based symbology

texts:
I am trying to categorise my site points by species group with rule based symbology;

e.g. bat uses the formula

"Group"='Bat'

But I would also like to have a different symbol for where there are no bats recorded at that site (i.e. based on the abundance in a separate column which is taken from a join with another table). I have tried this below but it hasn't worked, does anyone know how to fix this?

"Group"='Bat' AND "records_Abundance" = '<= 0' (for sites where species not present)

"Group"='Bat' AND "records_Abundance" = '>= 1' (for sites where species are present)

It seems to only read the abundance part and either turns all my points for all species to the symbol or none at all?


texts: Is it possible to 'unplace' annotations?

texts:
Is it possible to add or return annotations to the Unplaced Annotation list? Is there no way to put an annotation in the list once the annotation is placed?

I have noticed that the Aviation Editing toolbar has the option to hide annotation features, but there does not appear to be a button to add annotations back to the Unplaced Annotation list in the regular editing toolbar.

Using a definition query would be a last resort as there are lots of labels. Creating a layer from selected features and editing those would be a better option, but it would be nice to be able to place and unplace annotations in an editing session.


texts: Best approach for indexing a SQLite database containing geoJSON spatial objects?

texts:
I'm building an iOS application for public release. My architecture currently includes Mapbox, Parse, geoJSON and SQLite for the local store. The app will store spatial objects locally as a string in a SQLite database, and when a user opens the map view the app will query for spatial objects that exist within the map extent. Parsing/querying the raw geoJSON data would be very costly so I want to index the tables that contain any spatial data.

What's the best approach for this?  Is it efficient enough to store the spatial object's bounding box (or envelope) in a separate column and index that column?  Or should I go a step further and split out the envelop into it's component xmin, xmax, ymin, and ymax values, store those in their own columns and index each of those columns?

Or perhaps there's another approach I haven't considered yet? Maybe I shouldn't be storing the geoJSON at all and only the geometries?  I haven't fully committed to any design decisions yet, so I was hoping to hear thoughts from other's with more experience before I start building out the db.


texts: Using ST_AsText() without session object

texts:
I am storing gps locations as feature in a restfull api. The locations are linked by foreign key to "place" objects in a flask/sqlalchemy/geoalchemy framework. I would like to print a user readable string representation of the points that I store using ST_AsText like so:

class GPSPoint(Base):

    __tablename__="gps_points"

    id = Column('id', Integer, primary_key=True)
    geom = Column(Geometry('POINT'),  nullable = False)

    def __init__(self, geom):
        self.geom = geom

    def toJson(self):

        return OrderedDict([
            ("geom", self.geom.AsText())
            ])

However the AsText() method does not actually return a human readable string, but some sort of function object.


texts: Importing Excel into QGIS as layer and not as attribute table

texts:
I'm trying to add Excel sheet as a point layer, that includes coordinates into QGIS.

Instead of a shapefile, it comes as attribute table (layer name is shown and attribute table logo beside it instead of point layer) and the layer doesn't appear in the canvas. What I am doing wrong?


texts:
I am able to run the following polygon filter against my points and it will correctly filter out my points in POSTGIS="2.4.4 r16526"

select propertypost_propertypost.*
From propertypost_propertypost
where st_contains(
  'POLYGON((1616523 -172630,1618559 -170395,1612691 -169037,1610057 -170554,1609857 -175464,1611294 -176382,1616523 -172630))',
  propertypost_propertypost.location)

however when I run it in Linux with POSTGIS="2.4.3 r16312", it first give me  the following error

"ERROR:  Operation on mixed SRID geometries"

so I would fix it as

  st_polygon('POLYGON((1616523 -172630,1618559 -170395,1612691 -169037,1610057 -170554,1609857 -175464,1611294 -176382,1616523 -172630))',4326),
  propertypost_propertypost.location)

then I get this error:

"function st_polygon(unknown, integer) is not unique"

so I fixed it as

where st_contains(
  st_polygon('POLYGON((1616523 -172630,1618559 -170395,1612691 -169037,1610057 -170554,1609857 -175464,1611294 -176382,1616523 -172630))',4326)::geometry,
  propertypost_propertypost.location)

but this error wont go away,

I was wondering how I could get rid of this error


texts: operator does not exist: geometry  geometry

texts:
I'm running PostGIS 2.0 for PostGreSQL 9.2, and I have a customer who wants to know how close a provided geometry is to the nearest feature in a table. So I'm playing around with the "<->" and "<#>" operators.  My problem?  The "<#>" operator isn't recognized.

This works:

SELECT gid,
    the_geom <-> ST_GeomFromEWKT('<my geometry>')
    FROM "<my table>"
    ;

This does not:

SELECT gid,
    the_geom <#> ST_GeomFromEWKT('<my geometry>')
    FROM "<my table>"
    ;

I get an error:

ERROR:  operator does not exist: geometry <#> geometry

Am I missing something obvious?


texts: FeatureLayer creation using ArcGIS JS Api with data hosted on ArcGIS online

texts:
I have a few feature layers hosted on ArcGIS Online and I am trying to add feature layer editing to my web project. The JS API guide says the data must come from an ArcSDE geodatabase (https://developers.arcgis.com/javascript/jshelp/inside_editing.html), but i am hoping they just aren't telling us everything...

I have done a bit of research on this and can't really find any info on it or examples of anyone doing this.

Does anyone know if its possible to do feature creation/editing from a feature service hosted on ArcGIS online or have any experience with this?


texts: r.sun model crashing (GRASS GIS 7.0.4, Windows 7 Enterprise SP1)

texts:
I have been trying to run the r.sun model on a small test area but it crashes after roughly 5-6 minutes. A windows error pops up saying the program stopped working. I used it before to create an incidence angle raster map with no issue and r.horizon runs fine with the same inputs. Does anyone have a suggestion on what I might be doing wrong?

This is the code I use to run the model:

r.sun elevation=DEM@User aspect=Aspect@User slope=Slope@User horizon_basename=horangle horizon_step=11.25 glob_rad=Global day=15 step=1

After I click the windows prompt to close the running script it prompts me that the command is finished but it procuded no output. Below is a screenshot of the issues.


texts: Converting georeferenced .pdf to .tiff in ArcGIS for Desktop?

texts:
I am attempting to convert a USGS georeferenced pdf to a tiff.  I found this help page, which looks pretty simple: https://desktop.arcgis.com/en/desktop/latest/tools/conversion-toolbox/pdf-to-tiff.htm#GUID-E7643B8C-5493-4E07-83F8-6322307415D2

Except my Conversion Tools do not have "From PDF" (I am running Arcmap 10.2.2).

Does anyone have advice on how to troubleshoot or use another method?


texts: Exporting Arcgis Online content as georss or geojson?

texts:
There are number of sites talking about how to add a georss feed to an ArcGIS Online map (see silverlight and flex examples from esri_webapi), but I haven't located anything about exporting a georss feed, or geojson for that matter. Is this possible? how?

Just to be clear, I'm not talking about the basemaps, but the feature layers.

Update: correction again, the silverlight and flex examples linked to aren't for ArcGIS Online. They're for running on your own web server, which can then be imbibed by AGO. However if I had a server for this I wouldn't be using AGO in the first place. I'm beginning to the think the whole question is based on a false premise.


texts: How can I set brightness in openlayers 3.11....?

texts:
Current version of OL 3 is 3.11.2 , I've seen in old example (3.8) there is support for brightness / contrast :
http://openlayers.org/en/v3.8.1/examples/brightness-contrast.html

But newer versions don't support such feature.

Any ways to affect the layers brightness/contrast ?


texts: Getting subtypes from standalone table in ArcObjects

texts:
I have been working in the next code to get all Subtypes from a standalone table in a single string that concatenates it.
I work with ArcMap and arcobjects library in C# (sharp) language.

But I can not put it to works it correctly.

My actual code:

IDocument docRDC = m_application.Document;
IMxDocument mxdocRDC = docRDC as IMxDocument;

string esFeature = "";
string esTable = "";

if (mxdocRDC.SelectedItem is IFeatureLayer2)
    {
    esFeature = "Si";
    MessageBox.Show("Es IFeaturelayer2");
    }

if (mxdocRDC.SelectedItem is ITable)
    {
    esTable = "Si";
    MessageBox.Show("Es Itable");
    }

ITable table = mxdocRDC.SelectedItem as ITable;

System.Text.StringBuilder listadosubtipos = new System.Text.StringBuilder();

if (esFeature != "Si" && esTable=="Si")
            {
                IObjectClass MiObjeClassTabla = table as IObjectClass;
                ISubtypes MisSubtypesTabla = (ISubtypes)MiObjeClassTabla;
                IEnumSubtype enumSubtypeTabla;
                int subtypeCodeTabla;
                string subtypeNameTabla;

            if (MisSubtypesTabla.HasSubtype)
                {
                enumSubtypeTabla = MisSubtypesTabla.Subtypes;
                subtypeNameTabla = enumSubtypeTabla.Next(out subtypeCodeTabla);

                while (subtypeNameTabla != null)
                    {
                    listadosubtipos.Append(subtypeNameTabla);
                    listadosubtipos.Append(", ");
                    subtypeNameTabla = enumSubtypeTabla.Next(out subtypeCodeTabla);
                    }
                }

                textBox3.Text = Convert.ToString(listadosubtipos);
            }


texts: Deployment of Zeno Field 3.11 in windows 8.1

texts: Is it possible to install Zeno Field 3.11 in windows 8.1?


texts: Creating constant slope DEM using ArcGIS Desktop?

texts:
I wanted to create a DEM of constant slope say 1/3000.

How can I proceed in ArcGIS Desktop 9.3?

The area may be 1000m x 1000m or any other as per convenience.


texts:
I'm trying to make a model where I make a copy of 'geodatabase A' to another location as 'geodatabase B', delete all of the datasets within, and create new datasets for geodatabase B. The reason I'm copying a particular geodatabase (geodatabase A) is because it contains all the domains I need for geodatabase B. However, I want to delete the datasets once it's copied so that I can reassign the spacial reference. When I put the copy tool and the iterate datasets together in one model, it iterated the copy function as well, making the model go infinitely. So now I'm thinking I need to put the iteration as a submodel. However, I'm not sure how I can reference the input for the iteration, as the input doesn't exist yet.

Can someone help me or let me know a better way to do this?


texts: Adding a logo to a panel on an app in Google Earth Engine?

texts:
How can I add a logo to a panel on an app  in Google Earth Engine?

I don't even know how to upload a JPG or PNG to an asset and then add that image to the UI. Panel that I have on the right. Similar to the image below.


texts: Text Disappears After Creation using ArcMap?

texts:
I am using ArcMap 10.5.

I'm just trying to add some text, and I have no idea why but for the last maybe 30min, after I create them, they just disappear. I have tried other draw features (lines, polygons), restarting my map, trying different mxd's, restarting ArcMap and even restarting my computer.

I am using the 'Drawing' toolbar.

I am in Data View.


texts:
I´d like to integrate the for loop in my model in ModelBuilder in QGIS3 (similar as the for interator in ArcGIS ModelBuilder).
data:
testdata: point layer with different classes (assignment to different classes)
hexagons: polygon layer (hexagrids)
aim:
For each class of the point layer testdata calculate the number of points per hexagrid and assign the count to the centroid of the grid where the points are situated.
The result should be again 1 point layer (hexagrid centroids) with the number of points per hexagrid for each class as fields in the attribute table. If a hexagrid contains points from different classes, for each class a separate feature (separate row in the attribute table) should be available.
Example:

I suppose the for-loop (iteration) is not integrated as a tool in the QGIS ModelBuilder so a PyQGIS script is necessary. I´ve tried something but it didn´t worked properly.
My idea (without using PyQGIS) was the following one:

select hexagrids with at least 1 point and create a new layer from the selection (extract by location)
create centroids for the hexagrids
join the grid ID to the points (join attributes by location)
join the points to the centroids based on the grid ID (join attributes by field)
split the point layer testdata by attribute (class)
run the tool "Count points in polygons" in the batch mode for each splitted point layer
merge the created polygons together (1 layer)
join the merged hexagrid layer with the number of points per grid to the centroid layer (grid ID as the join field)

(alternatively to 6. and 7., join the hexagrid layers with the point count with the centroid layer in the batch mode and then merge them together into 1 layer).
This workflow functions with a limited number of layers (point classes) but for thousands of classes this is impossible so a PyQGIS script seems as the only appropriate sollution.
Could you help me how the python script should look like or do you have any other suggestions?

Here is a picture for better understanding:

Each centroid (red) of the hexagrid represents the number of points within the same hexagrid. The assignment to the centroid is based on the hexagrid ID (join attributes by location). If a hexagrid contains points from different classes, for each class a new centroid is created (see hexagrids in circles).


texts: Table of Contents(TOC) Widget for ArcGIS JavaScript API in Legacy?

texts:
I want to develop a legend using Table of Contents(TOC) Widget for ArcGIS JavaScript API in Legacy code.

I have tried its development by following the steps in the documentation of TOC Widget for ArcGIS JavaScript API but it seems that they are only for AMD code type. Also I haven't been able to find any example of this widget in Legacy code.

Does somebody know if it is possible? Is there any documentation about it?


texts: I downloaded the GPX to shapefile tool and when I use it I get the following error message: 1 No module named re Failed to execute (gpx2shp).


texts: How to programmatically check if the number of shapes = number of table records?

texts:
I have a handful of approximately 1000 shapefiles that are corrupted (see attached error message).  The shapefiles were generated from eCognition Developer 8.  There is a script tool that seems to repair the shapefile once it is identified as corrupted.

Edit:

I want to create a quick script to loop through all of my shapefiles and check if the number of shapes matches the table records. I can count table records using the following:

# Name: fcCount.py
# Purpose: calculate the number of features in a featureclass

# Import system modules
import arcpy
from arcpy import env

env.workspace = "C:/data"
Sample = "MyShp.shp"
result_dbf = int(arcpy.GetCount_management(Sample).getOutput(0))
print result_dbf

I would ultimately like to create some sort of logic check such as:

if result_dbf = result_shp:
    pass
else:
    print "There is a problem with" + str(Sample)

How can I count shapes directly without accessing the .dbf file?  Or, in other words, what is the best way to programmatically check if the number of shapes matches the number of table records?


texts: Image Classification Toolbar - Training Samples Manager Unselectable - ArcGIS

texts:
I am having an issue with this toolbar.
I have it enabled but many options are greyed out. The training samples manager is unselectable.
I have two raster files I just created from Landsat data that I want to create classes on by using the image classification toolbar.

(I have spatial analyst enabled and I just upgraded from ArcGIS 10.3 to 10.4.)

Here is my Table of Contents:

and the Image Classification Toolbar:


texts: How to reduce a custom bin histogram for a region in Google Earth Engine?

texts:
I'd like to reduce all the conties' data histogram with a custom bin, but what I got is a auto bin histogram, it seems ee.Reducer.fixedHistogram() can help, but I don't know how to pass the parameters(min,max,steps) for all conties.

  // load precip data (mm, daily total): 365 images per year

var precipCollection = ee.ImageCollection('IDAHO_EPSCOR/GRIDMET')
                    .select('pr')   // select  precip band only
                    .filterDate('2017-01-01', '2017-12-31');
print(precipCollection, 'precipCollection');

// load precip data (mm, daily total): 365 images per year

var precipCollection = ee.ImageCollection('IDAHO_EPSCOR/GRIDMET')
                    .select('pr')   // select  precip band only
                    .filterDate('2017-01-01', '2017-12-31');
print(precipCollection);

  // reduce the image collection to one image by summing the 365 daily
  rasters

var annualPrecip = precipCollection.reduce(ee.Reducer.sum());
print(annualPrecip);

// equivalent shortcut syntax

var annualPrecip2 = precipCollection.sum();

// visualize annual precipitation
var precipPal = ['white','blue'] // store palette as variable
Map.addLayer(annualPrecip, {min: 60, max: 3000, palette: precipPal}, 'precip');

  // load regions: counties from a public fusion table, removing
  non-conus states
      // by using a custom filter

var nonCONUS = [2,15,60,66,69,72,78] // state FIPS codes that we don't want
var counties = ee.FeatureCollection('ft:1ZMnPbFshUI3qbk9XE0H7t1N5CjsEGyl8lZfWfVn4')
        .filter(ee.Filter.inList('STATEFP',nonCONUS).not());
print(counties, 'counties');

  // visualize

Map.addLayer(counties,{},'counties');
var steps = 128;

  // get histogram precipitation values by county polygon(how to obtain custom      bin histogram,such as 64/128/256 bins)

var countyPrecip = annualPrecip.reduceRegions({
  collection: counties,
  reducer: ee.Reducer.autoHistogram(),
  scale: 4000 // the resolution of the GRIDMET dataset
});

print(countyPrecip);

// drop .geo column (not needed if goal is tabular data)

var polyOut = countyPrecip.select(['.*'],null,false);

// add a new column for year to each feature in the feature collection

polyOut = polyOut.map(function(feature){
return feature.set('Year',2017);
});

// Table to Drive Export Example

Export.table.toDrive({
  collection: polyOut,
  description: 'GRIDMET_annual_precip_by_county',
  folder: 'GEE_geohackweek',
  fileFormat: 'CSV'
});


texts: IMapServerObjects.get_Map says object is null reference.

texts:
So I am trying to access a map layer, which I can connect to.

try
                {
                    agsconn.Connect();
                    Log.Write("Connected.", "myclass");
                    if (agsconn.IsConnected)
                    {
    ESRI.ArcGIS.Server.IServerObjectManager som = agsconn.ServerObjectManager;
                            IServerContext serverContext = som.CreateServerContext(mapDocName, "MapServer");
                            IMapServer mapServer = serverContext.ServerObject as IMapServer;
                            Log.Write("Server context established MapName=" + mapServer.DefaultMapName, "myclass");
                            ESRI.ArcGIS.Carto.IMapServerObjects mapSO = mapServer as ESRI.ArcGIS.Carto.IMapServerObjects;
                            Log.Write("IMapServerObjects created", "myclass");
                            IMap map = mapSO.get_Map(mapServer.DefaultMapName);
                            Log.Write("Layers map established", "myclass");
                            //int layerIndex = -1;
                            IFeatureLayer fLayer = map.get_Layer(0) as IFeatureLayer;

The map was created with direct connect, I don't know if that makes a difference.

So far, the arcgis connects, and everything goes fine, but it stops here:

IMap map = mapSO.get_Map(mapServer.DefaultMapName);

I checked mapServer.DefaultMapName, and it returns "Layers" so that's all correct.

mapSO also works because mapServer and mapSO are the same object basically.

So how is there a null here?

My log prints out an exception and never reaches "layers map established".

The exception:

Object reference not set to an instance of an object.

I've tried a few other maps too and it doesn't work. But this code is straight from sample code I've seen on ESRI website, and it works for ESRI.

I'm using ArcGIS 9.3.1.


texts: I have three files with RRD extension, which are data files for ASTER and I do not know how to open them in ENVI 5.3?


texts: Rasterio IO Error: Read Failed and Gdal VRT_SHARED_SOURCE

texts:
I trying to follow some code on git which uses sentinel hub images code

I'm stuck with one part of the code which opens and reads a grayscale tif image which for some reason returns a rasterio IOerror.

with rasterio.open(img_path) as src:
    gamma0 = src.read(1)

The error message is as follows:

rasterio.errors.RasterioIOError: Read or write failed. /home/vulcanadmin/landstat-sentinel-fusion-installation/landsat-sentinel-fusion/data/processed/sentinel-1/antananarivo/vv-gamma0.tif, band 1: IReadBlock failed at X offset 0, Y offset 0: TIFFReadEncodedStrip() failed.

I checked the image. It isn't corrupt, no issues opening it with gdal.

Read that it has do with setting a GDAL variable, VRT_SHARED_SOURCE=0. I'm trying to figure out how this variable should be set or if anyone has encountered a similar issue and found a resolution. I'm using a linux system so tried setting it as an environment variable export $VRT_SHARED_SOURCE=0 but that didn't work. Neither did it work with the following,

with rasterio.Env(VRT_SHARED_SOURCE=0):
    with rasterio.open('vv-gamma0.tif') as src:
        gamma0 = src.read(1)

Would appreciate any help. Thanks.


texts: Invalid SRID return from PostGIS ST_Centroid Function

texts:
I'd like to get the centroid of a set of polygons with SRID 4326
I start from this simple query to verify the ST_Centroid function
CREATE TABLE test AS
SELECT ogc_fid, geom, ST_Centroid(geom) centriod_geom
FROM s_table

However, when I check the SRID of two geometry columns in the test table, it appears the SRID of centriod_geom changes to 0.
SELECT Find_SRID('public', 'test', 'geom');
SELECT Find_SRID('public', 'test', 'centriod_geom');

geometry column
find_srid

geom
4326

centriod_geom
0

What should I do if I need the centriod_geom to be 4326?


texts: Update cursor with not populating table with describe/xml metadata

texts:
I am writing a script that will export feature class metadata to a GDB table. The metadata is first obtained through arcpy.Describe and by parsing an exported XML file. The script then creates a GDB table, adds fields, and uses arcpy.da.UpdateCursor to populate the fields. The script runs and entering the various variables in the PythonWin interactive window returns the correct values, but they are not writing to the GDB table, so I end up with a table that has all the correct fields but no information. What is wrong with my UpdateCursor?

Code:

import arcpy
import csv
from xml.etree.ElementTree import ElementTree
from time import strftime
arcpy.env.overwriteOutput = True
# Set Translator folder path
AGSHOME = arcpy.GetInstallInfo("Desktop")["InstallDir"]
xslt = AGSHOME + r"Metadata\Stylesheets\gpTools\exact copy of.xslt"
# Get date and time
dateTime = strftime('%d%b%Y_%H%M')
# Parameters
data = r'T:\FS\Reference\GeoTool\r06\DatabaseConnection\r06_csa2_default_as_myself.sde\S_R06_CSA2.PortOrfordCedar\S_R06_CSA2.PortOrfordCedar_2013'
workspace = r"T:\FS\NFS\R06\Program\ResourceInfoMgmt\GIS\CSA2\Workspace\astratton\R6\SDE_Metadata_Python\Test"
# Get name for new folder in Workspace
descData = arcpy.Describe(data)
outputName = descData.name.replace('.', '_') + dateTime
# Create folder in Workspace
folder = arcpy.CreateFolder_management(workspace, outputName)
desc_fold = arcpy.Describe(folder)
fold_path = desc_fold.path
# Create GDB
table_GDB = arcpy.CreateFileGDB_management(folder, "Tables")
desc = arcpy.Describe(data)
path = str(desc.path)
file = str(desc.file)
print "Working on " + file
dataType = str(desc.dataType)
if desc.dataType != "Table":
    sr = desc.spatialReference
    sr_name = str(sr.name)
else:
    sr = "None"
# Export metadata to xml
outXML = "Temp_Metadata"
outFileXML = os.path.join(fold_path, ("{0}.xml".format(outXML)))
arcpy.XSLTransform_conversion(data, xslt, outFileXML)
tree = ElementTree()
tree.parse(outFileXML)
# Check for presence Metadata 10
spot = tree.find("dataIdInfo/idCitation/date/pubDate")          # Meta 10
if spot != None:
    Meta_10 = "Yes"
else:
    Meta_10 = "No"
# Check for presence Metadata 9
spot = tree.find("idinfo/citation/citeinfo/pubdate")            # Meta 9
if spot != None:
    Meta_9 = "Yes"
else:
    Meta_9 = "No"
# Get title (10 and 9)
spot = tree.find("dataIdInfo/idCitation/resTitle")              # Title
if spot != None:
    title = spot.text
else:
    title = ""
spot = tree.find("idinfo/citation/citeinfo/title")
if spot != None:
    title = spot.text
else:
    pass
# Get publication date (10 and 9)
spot = tree.find("dataIdInfo/idCitation/date/pubDate")          # Publication date
if spot != None:
    pub_date = spot.text
else:
    pub_date = ""
spot = tree.find("idinfo/citation/citeinfo/pubdate")
if spot != None:
    pub_date = spot.text
else:
    pass
table = arcpy.CreateTable_management(table_GDB, file.replace(".", "_"))
fields = ["Path", "FileName", "Data_Type", "sr", "Meta_10", "Meta_9", "Title","Pub_date"]
for field in fields:
    arcpy.AddField_management(table, field, "TEXT")
with arcpy.da.UpdateCursor(table, fields) as cursor:
    for row in cursor:
        row[0] = path
        row[1] = file
        row[2] = dataType
        row[3] = sr_name
        row[4] = Meta_10
        row[5] = Meta_9
        row[6] = title
        row[7] = pub_date
        cursor.updateRow(row)

del cursor
print "Done."


texts: Resolving "ERROR: could not open extension control file" from "create extension postgis;' for PostgreSQL 10?

texts:
I am receiving this error when I create extension postgis;:

  SQL Error [58P01]: ERROR: could not open extension control file
  "C:/POSTGR~1/10.5/pg10/../pg10/share/postgresql/extension/postgis.control":
  No such file or directory

The cause appears to be an incorrect path (an extra postgresql). It should be:
"C:/POSTGR~1/10.5/pg10/../pg10/share/extension/postgis.control"

This PostGIS installer that I used is: http://download.osgeo.org/postgis/windows/pg10/postgis-bundle-pg10x64-setup-2.5.0-1.exe

And I am running PostgreSQL 10.5 on x86_64-pc-mingw64, compiled by gcc.exe (Rev5, Built by MSYS2 project) 4.9.2, 64-bit

While my question is a near duplicate of https://gis.stackexchange.com/a/71476/31, the issue is specific to PostgreSQL 10/PostGIS 2.5, and I don't think the answer is necessarily the same (since there is no spatial_ref_sys.sql file).

Is there way to correct the error? If not, which files in the extension directory have to run (and in which order) to create the PostGIS extension?


texts: Accessing individual raster bands using arcpy Raster class?

texts:
I am working with rasters in arcpy and I am not able to access individual bands using the Raster class. The same code that worked about a year ago no longer works.

I should be able to access a single band of a raster like so:

Band3 = Raster(raster_filename + "/Band_3")

But the resulting Band3 variable has all 6 bands in the raster, not just band 3.

The "/Band_3" syntax appears to work fine in the Describe method in the first half of the following code but doesn't work in the Raster class later on:

    folder = r"C:\Users\elyons\Desktop\Spring_2018\Labs\New_LabData"

    import arcpy
    from arcpy.sa import *
    arcpy.CheckOutExtension("Spatial")

    arcpy.env.workspace = folder
    arcpy.env.overwriteOutput = True

    raster_filename = "Landsat.tif"

    #Create describe objects
    desc = arcpy.Describe(raster_filename)
    desc_band = arcpy.Describe(raster_filename + "/Band_3")

    #These are properties of the describe raster BAND object
    print "Height = " + str(desc_band.height)
    print "Width = " + str(desc_band.width)
    print "Resolution = " + str(desc_band.meanCellWidth)

    #This one is a property of the describe raster dataset object
    print "Number of Bands is = " + str(desc.bandCount)

    #Open bands as raster band objects
    Band3 = Raster(raster_filename + "/Band_3")

    #print number of bands
    print Band3.bandCount

    #Save the raster
    print "saving"
    Band3.save("savebandtest.tif")
    print "saved"

    arcpy.CheckOutExtension("Spatial")

The above code prints out the correct information for the describe objects and then prints 6 for Band3.bandCount, which is not right. It then saves a six band raster to the filename savebandtest.tif that is identical to Landsat.tif. I have tried different multiband rasters, different computers, getting the band name programatically, setting the working directory to the raster, forward slashes, backward slashes, everything I can think of.

Has the syntax changed for accessing single bands from multiband raster files? I can't find any mention of a change in the help.


texts: Layer leaf nodes replicated across all tree nodes using GeoExt?

texts:
This is going on with GeoExt 1.1 and OpenLayers 2.12. I need to build a layer tree dynamically, with layers and nodes that are only known at run time. The snippets I show below correctly create a layer tree with several nodes and adds the layers node leafs. The problem is that though each layer is assigned to a single tree node, it is then shown under each and every tree node created.

During init the tree and its root are set up:

  treeRoot = new Ext.tree.TreeNode
  ({
      text        : 'Project',
      id          : '0',
      leaf        :false,
      expandable  : true,
      expanded    : true
  });

  layerTree = new Ext.tree.TreePanel
  ({
    title: 'Map Layers',
    region: "west",
    width: 200,
    collapsible: true,
    autoScroll: true,
    enableDD: true,
    plugins:
    [{
      ptype: "gx_treenodecomponent"
    }],
    loader:
    {
      applyLoader: false,
      uiProviders:
      {
        "custom_ui": LayerNodeUI
      }
    },
    root: treeRoot,
    rootVisible: false,
    lines: false
  });

Later on new tree nodes can be added according to a tag, and kept in an array:

if(treeNodes[tag] == null)
{
    treeNodes[tag] = new GeoExt.tree.OverlayLayerContainer
    ({
            text: tag,
            leaf: false,
            expanded: true
    });
    treeRoot.appendChild(treeNodes[tag]);
}

Layers are added to the map and to a specific tree node:

var layer = new OpenLayers.Layer.WMS(title, serviceURL, params, options);
leftMap.addLayer(layer);
layer.events.register("visibilitychanged", this, toggleLayer);

var newNode = new GeoExt.tree.LayerNode(
{
    text: title,
    layer: layerName,
    leaf: true,
    checked: false,
    children: [],
    nodeType: "gx_layer"
});
treeNodes[tag].appendChild(newNode);

But then each layer appears several times, once under each tree node.

Why is it happening this way?


texts: I have implemented a new web map application in my organization and  I want to do a survey/questionnaire for the web map users (both internal and external users) to collect information about what they think about the web map and if they have feedback about it. I will be thankful if you could provide me with tips/suggestions on the type of questions that I should ask in my questionnaire. Has anyone done a similar survey before about their GIS systems?


texts: GEE, pixel based SD over time series Sentinel 2 NDVI

texts:
I would like to calculate the SD of a NDVI value by pixel over a time series using Sentinel 2. It is to identify pixels with more or less variations over a specific time period. The expected output would be a collection of images where pixels values are SD,mean, Min MAX, etc.

What I've done so far is :

//Import GEE Feature Collection
    var ROI = ee.FeatureCollection('https://code.earthengine.google.com/?asset=users/floriangirond/multiplearea');

// Create image collection of S-2 imagery for the perdiod 2018-2019
var S2 = ee.ImageCollection('COPERNICUS/S2')

//filter start and end date
.filterDate('2018-01-01', '2018-12-31')

//filter according to drawn boundary
.filterBounds(ROI);

// Function to mask cloud from built-in quality band
// information on cloud
var maskcloud1 = function(image) {
var QA60 = image.select(['QA60']);
return image.updateMask(QA60.lt(1));
};

// Function to calculate and add an NDVI band
var addNDVI = function(image) {
return image.addBands(image.normalizedDifference(['B8', 'B4']));
};

// Add NDVI band to image collection
var S2 = S2.map(addNDVI);
// Extract NDVI band and create NDVI median composite image
var NDVI = S2.select(['nd']);
var NDVI = ee.FeatureCollection(NDVI.map(mapfunc));
var reducer1 = ee.Reducer.mean();

var reducers = reducer1.combine({reducer2: ee.Reducer.median(), sharedInputs: true})
                       .combine({reducer2: ee.Reducer.max(), sharedInputs: true})
                       .combine({reducer2: ee.Reducer.min(), sharedInputs: true})
                       .combine({reducer2: ee.Reducer.stdDev(), sharedInputs: true});

var results = NDVI.reduceRegion({reducer: reducers,
                                geometry: ROI,
                                bestEffort: true});

print(results);


texts: How to print 2 maps in the same PDF page using PyQGIS?

texts:
I am trying to generate a one-page PDF document containing the following items:

1) On the left half, a full state map (say, New York State).

2) On the top right quadrant, a map that zooms in to one part of the map in item (1). For example, the Manhattan area.

3) On the bottom right quadrant, a table with text and some data.

I could generate all items individually but couldn't put both 1) and 2) on the same page at once. Both item 1) and 2) come from the same layer but rendered based on different extent. If I'm not mistaken, QgsComposition only takes in one instance of QgsMapRenderer at a time so only one extent is possible?

Any clue on what I should do? This is for a standalone PyQGIS app. I am using QGIS 2.4.0 on Ubuntu 12.04.


texts:
I want to show a pie presentation of the number of mosquitoes caught in each trap, where the size of the circle is the proportion of mosquitoes in one trap out of the total mosquitoes caught, and on top of that, I want to color each circle according to another parameter, which is how much more or less each point is in relation to the last trappings (in percentage, of proportion). For instance, a large circle for a big proportion of mosquitoes caught, and a red color for increase from last time. small circle for few mosquitoes and green for decrease from last time.

How do I combine the two?


texts: Looking for an open-source Java based GeoTIFF library

texts:
I have written many programs for reading and writing various types of spatial data in the past but I must admit that the GeoTIFF format has often gotten the better of me. It's such a flexible format that writing code to read/write the variations that are commonly found in the wild is extremely challenging. So I am looking for a library that I can use for accomplishing this task.

Here are my requirements:

It must be developed in Java and exist as a single and ideally very small JAR file that can be easily incorporated in existing projects.
It must have a licence compatible with the GNU GPL.
It must be able to handle many different flavours of GeoTIFF including the various types of data compression.

The GeoTIFF web page recommends GeoTIFF-JAI (it's also recommended by the answerer of this GIS.SE question) however the project doesn't appear to have been maintained since 2001 and I am not able to find any documentation to describe its use. I am currently using a modified version of the UCAR GeoTIFF reader but it similarly is not actively maintained nor documented, and does not handle many types of modern GeoTIFFs (no compression formats). I am aware that GDAL can handle GeoTIFF quite well and that there is a JNI interface to link with Java but I am looking for a small-footprint native Java solution. I am also aware that GeoTools has a means of dealing with GeoTIFFs, but as the OP of this question points out, it actually requires a GDAL plugin, which unfortunately I need to avoid.

I realize that I am asking for a lot here but is anyone aware of a GeoTIFF library that meets these conditions? I thought I'd ask before spending another couple of days trying to hack together a solution for dealing with compressed formats.


texts: Exporting a thematic raster in byte format to use with GUIDOS toolbox

texts:
I am trying to export some thematic rasters to use them with GUIDOS Toolbox (specifically the MSPA segmentation functions) - http://forest.jrc.ec.europa.eu/download/software/guidos/

GUIDOS requires the images to be GeoTIFF format, with pixels of data type Byte. It reads 0 bytes as no data (optional), 1 byte as background and 2 bytes as foreground and performs segmentation only on foreground pixels.

In Arc I have rasters with values 2=foreground, 1=background - although changing these is trivial. How do I go about converting these to byte and exporting to GeoTIFF? The lowest pixel depth Arc seems to support for TIFF is 8bit.

Does anyone have a workaround, or some third-party recommendations to do this pre-processing for GUIDOS?