texts: What would cause a geoprocessing service to say "Cannot Open 'FeatureClass' Failed to Execute"?

texts:
I have a really basic geoprocessing script that's published to ArcGIS Server 10.2.2.

It basically just reads in a json object that's passed in as a parameter.
The json just has attribute names and values, as well as a geometry in wkid 4326.

All it does is open an InsertCursor to create features in one specific feature class.

Nevermind asking me "All you want to do is create a bunch of features, why don't you do this in the rest endpoint?" Long story short, I tried several times using the jsapi, and even building my own geometries. In both cases, I got "d is not applied."

Anyway. This Geoprocessing script does work. But sometimes I check the logs and see "SEVERE: 'Cannot open SCHEMA.FEATURECLASS'. Failed to execute (Name of GP service). Failed to execute (name of GP Service).

What would cause this?? I noticed that on times where it works, it was very early in the day (5 AM). And then someone tried to use it around 10 AM and the script fails. From my DBMS knowledge, the only reasons I can think of why someone can't open a feature class for inserting:

Someone has a schema lock (I find this less likely since you usually can't get a schema lock on anything if the ArcGIS Server is running)
A database transaction is in progress (I don't see any evidence of this since the only other system that loads data into this database shows no records of any loading happening around the time of the error log
does arcpy.da.insertcursor fail to open if there's another insertcursor open simultaneously? (I hope not...)

Other info:

This is not a versioned feature class
The arcpy.env.workspace is set to an sde file that connects as the geodatabase admin. This connection is set to use the default version
Here are the properties on the service

#CODE SNIPPET - This is 90% of the business logic here
loadCursor = arcpy.da.InsertCursor(Feature_Name, fields)
oids = []

for feature in FeaturesToAdd:
    geom = feature["geometry"]
    point = (geom["x"], geom["y"])
    attributes = feature["attributes"].values()
    attributes.append(point)
    row = tuple(attributes)
    print len(row)
    oid = loadCursor.insertRow(row)
    oids.append(oid)

del loadCursor

print oids

What else would cause arcpy not to be able to open a feature class for inserting?

Here's the basic code snippet for doing that:

loadCursor = arcpy.da.InsertCursor('SCHEMA.FEATURECLASS', ['field1','field2', 'field3', 'shape@XY']

there's nothing else that opens up any cursors, lists fields, queries the db. etc. This is the only line (not including loadCursor.insertRow())


texts: Help using the US American Community Survey Data

texts:
I want to use the USA American Community Survey Data for spatial analysis, but I'm having difficulty understanding how the survey tables are constructed so that I can relate them to block-group shapefiles.

I've successfully made a shapefile that contains all the US blockgroups, but I can't figure out how the census data is related to it using the STFID code. I have a few thousand text files (named e20095ak0001000.txt for example). When I read them into R they appeared with the following headers:

[1] "ACSSF"    "X2009e5"  "ak"       "X000"     "X0011"    "X0000001"
[7] "X4269"    "X2211"    "X224"     "X157"     "X189"     "X142"
[13] "X77"      "X269"     "X247"     "X169"     "X389"     "X146"  ... up to 200 or so

(I'm assuming the blockgroup code does not match the filename, as I have 12,465 files and 212,083 blockgroups)

I realize I need to figure out what the the ID code is, to link these to the block groups, but I can't find documentation that explains this information (even after reading the technical documentation). Where are these column headers explained?

I am using the 5 year ACS survey.


texts: R: Error when stacking multiple raster files

texts:
I have a large number of rasters in grid format ("asc"). I created a list of files:

files <- list.files("D:/......../EnvironData", pattern='asc',
                    full.names = TRUE)

Then I want to stack all files:

dataStack <- raster::stack(files)

Unfortunately, I get the following error:

  Error in .local(.Object, ...) :

  Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  :
  Cannot create a RasterLayer object from this file.

All raster files have the same spatial dimension and extent.


texts: GDAL TMS minidriver in QGIS, problem or bug?

texts:
Its great that you can now use TMS layers in qgis! Only, my result is broken. At the 3 highest resolution zoom levels it looks perfect, but at lower resolution layers, i get mixtures of tile shift, black tiles, bright noise colors, etc.

Heres examples of lower zoom (nb: vectors are seperate layers).
pic http://zoneblue.org/files/qgis-tms.jpg

Im using qgis 2.4.0. Here is my xml:

<GDAL_WMS>
    <Service name="TMS">
        <ServerUrl>file:///E:/gis/tilesets/ortho/${z}/${x}/${y}.jpg</ServerUrl>
        <SRS>EPSG:2193</SRS>
        <ImageFormat>image/jpeg</ImageFormat>
    </Service>
    <DataWindow>
        <UpperLeftX>1678531.2</UpperLeftX>
        <LowerRightY>5636060.0</LowerRightY>
        <LowerRightX>1703824.0</LowerRightX>
        <UpperLeftY>5669442.4</UpperLeftY>
        <TileLevel>9</TileLevel>
        <YOrigin>bottom</YOrigin>
        <SizeX>63232</SizeX>
        <SizeY>83456</SizeY>
        <TileCountX>247</TileCountX>
        <TileCountY>326</TileCountY>
    </DataWindow>
    <Projection>EPSG:2193</Projection>
    <BlockSizeX>256</BlockSizeX>
    <BlockSizeY>256</BlockSizeY>
    <BandsCount>3</BandsCount>
</GDAL_WMS>

The tileset has 0-9 levels, 0-246 for the folders, 0-325 for the files. 0.4m per px imagery. Everything is transverse mercator.

You just make a xml file with the above in it, add a raster layer with it, and its done. Note that the pixel size and coordinate extents must be for the whole tiles, not just your data extent.

Anyway am i doing something wrong, or is this a bug?

I have read and reread the gdal minidriver spec here, but its not very detailed.
http://www.gdal.org/frmt_wms.html


texts: Is there a source for common GIS terminology in French?

texts: I am trying to find a a resource that provides English/French translations for GIS terminology.  Obviously I could use Google or a French/English dictionary, but there are many cases where there are multiple French words.  For example, the word layer could be translated as couche (which seems to be the most common translation), or strate (which is a geology term and may be more pertinent when talking about GIS)


texts: Rename the zonal statics table same as the raster name in ArcGIS model builder

texts: I need to get zonal statics for multiple raster files. However, I would like the name of the output table to be the same as the raster name but i am getting errors.


texts: GRASS GIS: auto-completion fails in commande line

texts:
I'm running GRASS GIS 6.4.3. on Win 7, SP1, 64 bit.

The path to my grass-DB is: D:\grassdb

Within the GRASS GIS Commande Line, if I press the tab-key to auto-complete only folders, files and documents from "C:\" are listed and none of the maps from GRASS which should appear:

If I type:

r.colors map=D

and press the tab-key after the "D"
this appears:

r.colors map="data"
r.colors map="Documents and Settings"

and so on.

Instead I'd like to see my layers:

r.colors map=DTM_2012
r.colors map=DTM_2013

and so on.

Any hints?


texts: Feasibility of a custom neighborhood-summary calculation (seasonal asynchrony) in GEE?

texts:
The gist of my question is:
Is there an obvious way that I’m not seeing to implement my calculation on GEE, and if so, how?

Details
For a research project I am attempting to calculate a new Image, at global extent and identical resolution to the original input MODIS ImageCollection. In the new Image, each pixel’s value will indicate the rate of drop-off across space in the correlation of that pixel’s time series with time series of neighborhood pixels. I am calling this value the pixel’s ‘asynchrony’ value, such that the resulting Image of asynchrony values would be an ‘asynchrony Image’. (The ImageCollection I am using extends across multiple years, such that the asynchrony will be interpreted as ‘seasonal asynchrony’. Other filtering and analytical steps will be taken to support that interpretation, but those are not pertinent to the present question of this algorithm’s feasibility, so are not herein discussed.)
The pseudocode for the naive algorithm that would calculate the asynchrony Image is:

for each focal pixel pix_foc:

var R2 = []: empty vector for R2 values
var dist = []: empty vector for distance values
var neigh_foc = getNeighborhood(pix_foc, maxDist): set of the pixels within the focal pixel’s neighborhood, defined by a circular kernel of neighborhood radius maxDist
for each neighborhood pixel pix_i in neigh_foc:

var model = linRegression(ts_pix_i ~ ts_pix_foc): regression of the neighborhood pixel’s time series on the focal pixel’s time series
R2.append(model['R2']): store the R2 of that regression
dist.append(getDistance(pix_foc, pix_i): store the distance between the focal and neighborhood cells’ centroids

var overall_model = linRegression(R2 ~ dist): regress the vector of R2 values on the vector of distances
var asynch_foc = abs(overall_model['slope']): assign as the focal cell’s asynchrony value the absolute value of the (likely negative) slope of the overall regression

I understand that GEE does not provide a flexible, general-purpose scripting environment within which that naive algorithm can be implemented as laid out above. Thus I have tried to implement this algorithm within GEE’s client-server, map-reduce paradigm. I have written a function that takes an ImageCollection ('imgColl'), a max neighborhood-cell distance ('maxDist'), and a pixel's centerpoint coordinates ('lon' and 'lat'), and returns an asynchrony value for that pixel. To call that function for all pixels globally, I have then tried the following steps:

get ee.Lists of the lon and lat values for the centerpoints of all pixels in the ImageCollection (vars ‘lons’ and ‘lats’)
use two nested map methods to calculate asynchrony values for all pixels, by calling:

lats.map(function(lat){

  return lons.map(function(lon){

    return calcAsynchrony(imgColl, maxDist, lon, lat)})})

reduce, using ee.Reducer.first(), the ImageCollection of 1-by-n Images returned by that call (where each Image contains the horizontal strip of asynchrony values for a given latitude)
This approach feels really hacky and ugly. But more importantly, it won't scale, because of memory limitations on the List objects created.

I have spent a good amount of time reading through the GEE docs and looking for a set of methods and data structures that might allow me to implement this calculation. At this point, I feel fairly certain that there is no obvious way of hacking this on GEE, and that it is probably worth my time to export my data, run this calculation elsewhere, then import the result to continue analysis on GEE. However, before I make that decision and invest the extra time it will require, I want to make sure that I’ve exhausted my options.
Can anyone think of an approach that I’ve overlooked for running this calculation on GEE? And if so, can you lay out the steps of that approach, so that I can consider them and possibly test-run it?
Thanks for your time and consideration!


texts: How to draw nodes of a graph in Python and link them to a table?

texts: I would like to draw some points with a Python GUI, like a graph, and link those points to a table.


texts: Symbol's size huge in the layer tree legend in QGIS

texts:
When I use Type of symbol -> Generator of Geometry of a vector layer that the symbol has a huge size in the layer tree legend of a QGIS projet.

I use this formula wedge_buffer($geometry, "fyd" , "fov" , 0.00007), it displays well into the canvas, but it is very big into the layer tree legend.

I get the same issue with both QGIS 3.4 and 3.10, the CRS is WGS 84

Is there a way to avoid this display issue or should I forget about using wedge_buffer?


texts: I have my osm server with mapnik, it works fine but now I want to change the background color( black to white), how can I do that? Is it a server or client configuration?


texts: How to extract FileGDB's domains without ArcGIS?

texts:
I'm trying to get an ESRI File Geodatabase domains by using the GDAL\OGR library. By doing so, I wish to edit this datasource using QGIS' Attributes Form with designated value map widgets that point to this domains.

I'm familiar with QGIS\OGR's capabilities of opening, viewing and even altering File Geodatabase (using FileGDB driver). What's missing is the ability to extract the properties of the FGDB itself - such as 'domains'.

Is there any way to extract FileGDB domains only using OSGeo programs without the help of ArcGIS software?


texts: Assign projection to hundreds of shapefiles using OGR2OGR

texts:
I'm trying to assign projection ESPG:32638 to hundreds of shapefiles in batch way like this

for %%N in (d:\lid\i08*.shp) DO ogr2ogr -a_srs epsg:32638  d:\lid\i08\mod\%%N d:\lid\i08\%%N

however, this is not working

Any clues why it's not working?


texts: Accessing Microsoft Azure File Storage from ArcMap?

texts:
I need to store some raster files on Microsoft Azure account.  They are rather large with each file being around 15MB.  It seems the best place on Azure is in File Storage as blobs (unless there is a better place?).

Will I be able to access from inside ArcMap via a Folder connection or a database connection?


texts:
I'm running a for loop with a cursor to disaggregate a large feature class into many smaller ones (of one polygon each). ArcMap crashes after about 1100 iterations, and I am trying to break the job down into more manageable snippets. I just can't get the code right. I have 4909 entries, and want to do it in five sets of 1000.

Below is a copy of my original code.

How do I change it so that it runs for a selection of rows at a time?

fc='C:\\Users\\xxxxx\\popcenters\\miyun_villages_poly.shp'
cursor=arcpy.SearchCursor(fc)
for row in cursor:
    FID = str(row.getValue("FID"))
    where = '"FID" =' + FID
    out_fl='C:\\Users\\xxxx\\miyun_allvillages.gdb\\allvillages_miyun\\fl_{0}'.format(FID)
    arcpy.MakeFeatureLayer_management(fc,out_fl,where)
    SaveToLayerFile_management(out_fl, allvillagesfl)


texts: Is there a python option to "join attributes by location"?

texts: I am trying to perform the function join attributes by location as found on the QGIS menu Vector>Data Management Tools. I am seeking an open source python option for this. I know arcpy has a spatial join function but I'm trying to perform this outside of the ESRI environment.


texts: why sql file generated by postgis was different in notepad and pgAdminIII?

texts:
When i imported a shapefile which contains some polygons into PostGIS by command line on Win 7, a problem puzzled me:
   The sql file generated by shp2pgsql opened with notePad has the geom (Mutipolygon) value ,however when put this sql file in pgAdminIII, the value was empty..why??some one can help me??


texts: Calculation amount or number of certain pixel-values in a raster

texts:
It seems pretty simple but until now I could not find a proper solution:
In a given raster i want to know the number or ratio of e.g. red pixels.

Edit:
to be more precise: I have a number a arial photos and want to estimate vegetation cover, so I tried to estimate/count the number or amount of green pixels. If there is an other solution to estimate vegetation cover, i am open to any other way.

First Attemps were raster statistics and inspecting histogramms, but the results were not clear to me.


texts: QGIS WMS/WMTS layers only renders as "Singleband color data"

texts: I'm trying to visualize some height rasters on QGIS 3.18 (Windows and Mac). The problem is that only the "Singleband color data" band rendering mode is available in the layer's properties/ symbology. I'm missing the modes normally available with raster data (Multiband color, Paletted/Unique values, Singleband gray, Singleband pseudocolor, Hillshade, Contours). I've tried several different WMS/WMTS sources and updated versions of QGIS 3.18, now on 3.18.2. Should all the modes be possible also with WMS/WMTS sources?


texts: Auto-Populating field based on another field in ArcGIS Online?

texts:
I am trying to use the Damage Assessment hosted feature layer and dashboard for local government.

One of the fields that we use a domain (combobox) for is the Inspector who looks over the damage. Is it possible to auto-populate the phone and email fields for the inspector once the object is created in ArcGIS Online?

I know there are ways to do this with the calculate field tool, but that is not automatic. I am honestly not sure if this is even possible.

Do you know if it is?


texts: Detecting and fixing outliers in a GPS trajectory

texts:
I need to find an algorithm or method that can detect outlier latitude longitude points in a trajectory during post-processing, which can then be fixed (brought back into the trajectory's path based on its neighbours).

As an example of the kind of outlier points I would like to detect and fix, I've attached an image demonstrating:

I have tried using an unscented Kalman filter to smooth out the data as best as possible, but this does not seem to work effectively enough for more extreme outliers (raw data in blue, smoothed data in red):

My UKF may not be calibrated properly (but I'm fairly certain that it is).

The trajectories are those of walkers, runners, cyclists - human-powered movement that can start and stop, but not drastically change in speed or position that quickly or suddenly.

A solution that does not rely on timing data (and only on position data) would be extremely useful (as the data being processed may not always contain timing data). However, I'm aware of how unlikely this kind of solution is to exist, so I'm equally as happy to have any solution!

Ideally, the solution would detect the outlier so that it could be fixed, resulting in a corrected trajectory:

Resources I've sifted through:

Smooth GPS data - https://stackoverflow.com/questions/1134579/smooth-gps-data
Common GPS and Geospatial Tracking Challenges and Solutions - http://www.toptal.com/gis/adventures-in-gps-track-analytics-a-geospatial-primer (solution seems to lose precision of data)
What algorithm should I use to remove outliers in trace data?


texts: Why is arcpy ChangePrivileges_management adding a second feature class to map?

texts:
I am new to GIS and pretty new to python, so sorry if this is ugly code.
I am creating a feature class in a postgres database if it does not already exist and the user wants it created.
That works fine, and adds a single layer to my table of contents in my current map.
I added the line to change the privileges of the layer, and now I get a second polyline feature class added to the table of contents.  Named a little different, and it does not stored in the database (the first, correct polyline does get stored).  Where is this layer coming from, and how can I make it not show up?

   is_there = False
   matched_layer = False
   for fc in arcpy.ListFeatureClasses():
       regex = re.compile("skld_geodata\..+\.A%s" % sectional_name)
       if regex.match(fc):
           matched_layer = fc
           is_there = True

   if is_there is False:
       msg = 'No skld layer currently exists for this Arb.  Would you like to create one?'
       create_input = pythonaddins.MessageBox(msg, 'Missing Layer', 4)
       # create a new geodatabase entry with the correct name if it does not already exist
       if create_input == 'Yes':
           arcpy.CreateFeatureclass_management(shp_path, arb_name, 'POLYLINE')
           # change the privileges for everyone to work with
           arcpy.ChangePrivileges_management(shp_path + "\\" + arb_name, "editor", "GRANT", "GRANT")
   else:
       geo_layer = arcpy.mapping.Layer(matched_layer)
       arcpy.mapping.AddLayer(self.df, geo_layer, "TOP")

I can add images of the table of contents if that makes it easier to understand what I am asking.


texts: Missing vector features when import from .sqlite to postgis?

texts:
I imported vector data from a .sqlite file to PostGIS but not all of my vector features imported.

I haven't been able to determine a reason for most of the loss of features between the databases (being potentially erroneous attribution or lack of xy coordinates).

Why am I missing vector features in my database after the import from .sqlite file?

an example, the building layer with 78 features, the 2 features which do not migrate are pointID 221 and 262 being addressed 17896 buckhorn and 56430 buckhorn


texts:
So, in QGIS I'm trying to apply the script "Point layer to route".

(See this website for a general idea about the script: http://anitagraser.com/2013/10/19/a-routing-script-for-the-processing-toolbox/#comments)

The script takes as input: a network layer and point. The goal is to map-match the point to the network and subsequently create a route, from the points, on the network. Below your can see the network in grey, the points in red, and the route in red.

The process, from right to left, goes as it should, until the 2 point in the left lower corner. When the distance between the network and the points becomes less than 10 meters, the process isn;t exectued anymore. Do you guys have any idea on how to solve this?

The scipt is added below:

#Definition of inputs and outputs
#==================================
##Routing tools=group
##points=vector point
##network=vector line
##route=output vector line

#Algorithm body
#==================================
from PyQt4.QtCore import *
from PyQt4.QtGui import *

from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

from processing.tools.vector import VectorWriter

point_layer = processing.getObject(points)
network_layer = processing.getObject(network)
writer = VectorWriter(route, None, [QgsField("order", QVariant.Int)], network_layer.dataProvider().geometryType(), network_layer.crs() )

# prepare graph
vl = network_layer
director = QgsLineVectorLayerDirector( vl, -1, '', '', '', 3 )
properter = QgsDistanceArcProperter()
director.addProperter( properter )
crs = vl.crs()
builder = QgsGraphBuilder( crs )

# prepare points
features = processing.features(point_layer)
point_count = point_layer.featureCount()

points = []

for f in features:
  points.append(f.geometry().asPoint())

tiedPoints = director.makeGraph( builder, points )
graph = builder.graph()

route_vertices = []

for i in range(0,point_count-1):
    progress.setPercentage(int(100 * i/ point_count))

    from_point = tiedPoints[i]
    to_point = tiedPoints[i+1]

    from_id = graph.findVertex(from_point)
    to_id = graph.findVertex(to_point)

    (tree,cost) = QgsGraphAnalyzer.dijkstra(graph,from_id,0)

    if tree[to_id] == -1:
        continue # ignore this point pair
    else:
        # collect all the vertices between the points
        route_points = []
        curPos = to_id
        while (curPos != from_id):
            route_points.append( graph.vertex( graph.arc( tree[ curPos ] ).inVertex() ).point() )
            curPos = graph.arc( tree[ curPos ] ).outVertex()

        route_points.append(from_point)

    # add a feature
    fet = QgsFeature()
    fet.setGeometry(QgsGeometry.fromPolyline(route_points))
    fet.setAttributes([i])
    writer.addFeature(fet)

del writer


texts: Is there a described GIS method for creating catchment areas based on nearest driving distance?

texts: ie. a combination of Voronoi polygons with isochrones, so that the Voronoi polygons are based on driving distance instead of Euclidean distance. Is there a name or a described method for this?


texts: Zonal Statistics as Table output is missing some polygons ONLY when run from Python script

texts:
I have a raster file from which I am attempting to calculate mean values for polygons defined in a shapefile. Using the Zonal Statistics as Table tool from within ArcGIS works great. My shapefile contains 118 shapes, and the output table contains 118 rows.

My problem is that when I try to execute the exact same tool from a python script, not all of the shapes are included in the results. From 118 shapes in the shapefile, my output table only contains 87 rows.

My polygons do not overlap, and they are all much larger than the pixel resolution of the raster.

Here is the code I am running:

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

env.workspace = "E:/Temp"

shapeFile = "test_shape.shp"
shapeFileField = "FID"
inRaster = "test_raster.pix"
outTable = "test_output"

arcpy.CheckOutExtension("Spatial")

outStats = ZonalStatisticsAsTable(shapeFile, shapeFileField, inRaster, outTable, "DATA", "MEAN")

Can anybody give me any hints as to why this is happening?


texts: ArcGIS JS api - how to hide attributes editor when new graphic is added

texts:
I use ArcGIS Javascript api in 3.12 version. I have followed this sample in order to display editing toolbar plus template picker. It works pretty good. However I would like to hide popup window with attributes editor which is displayed after new graphic is added. That's because I already have attached to every FeatureLayer my own InfoTemplate:

    var template = new InfoTemplate();
    template.setTitle("Object");
    template.setContent(getPassportLink);
    var layer = new Dependencies.FeatureLayer(layers[i].serviceUrl, {
        mode : Dependencies.FeatureLayer.MODE_ONDEMAND,
        infoTemplate : template,
        outFields : [ "*" ],
    });

How can I achieve that?


texts: Copying attributes from one point to another in QGIS

texts: I have mapped several points on a shapefile. Some of them share the same attributes. In the attributes table I copied the attributes of the source point, select the target point with empty row, but after paste it, it creates a new row. How a can get my target point "filled" with the same attributes?


texts: How to make QGIS abbreviate labels with a dictionary?

texts:
I'm trying to use qgis "expression dialog" for labeling a layer. What I'm looking for is an expression that permits to apply some abbreviation rules on the layer.
For now i'm using this

CASE
WHEN nome2 is not null THEN  replace("nome1"  || '$' ||  "nome2",  'Wasserfall', 'Wssf.')
WHEN nome2 is not null THEN  replace("nome1"  || '$' ||  "nome2",  'CascatA', 'Casc.')
WHEN nome2 is not null THEN  replace("nome1"  || '$' ||  "nome2",  'Cascate', 'Casc.e')
WHEN nome2 is null THEN    replace("nome1", 'Wasserfall', 'Wssf.')
WHEN nome2 is null THEN    replace("nome1", 'Cascata', 'Casc.')
WHEN nome2 is null THEN    replace("nome1", 'Cascate', 'Casc.e')
END

Obviously it doesn't work because only the first condition (wasserfall) is used for the label.

It happens something like that

Wasserfall$Cascata -> Wssf.$Cascata

instead of

Wasserfall$Cascata -> Wssf.$Casc.

I can solve the problem directly in postgis (with a view or a simple update) but if it is possibile to do it in qgis it would be great.
Most of all it would be great to use an abbreviation dictionary but I'm quite sure that qgis (for now) is not able to doing it.
I'm on W7 qgis 2.4.


texts: How to rotate rasters with GDAL, QGIS or GRASS GIS

texts:
I need to rotate raster files from the TRMM Multisatellite Precipitation monthly product 3B43. For some unknown reason, this data comes rotated 90 degrees clockwise, meaning that the y coordinate is the x and that the x coordinate is the y.

Of the few proposed solutions which can be found online, none worked:

hdp dumpsds -  hdp dumpsds -r 4 -d -o outfile -b 3B43.20000501.7A.HDF.nc
error: SDS with reference number 1: not found
Use ENVI - can't use/don't have ENVI.
CDO specifying a target grid - This approach almost worked, however  while the data extent is correctly transformed (i.e., x and y are swapped), the pixels are placed randomly. My target grid specification are:

gridtype = lonlat

xsize = 1440 # replace by your value

ysize = 720 # replace by your value

xfirst = -179 # replace by your value

xinc = 0.25 # replace by your value

yfirst = -89.5 # replace by your value

yinc = 0.25 # replace by your value

Then run the commands:

cdo setgridtype,curvilinear  3B43.20000501.7A.HDF.nc out.nc
cdo remapbil,target_grid.txt out.nc out2.nc


texts: Confusing results from script that uses arcpy.clip analysis?

texts:
Following on from previous question Python code that uses arcpy.describe skips if statement?, I have now attempted to include clip analysis and this is my result;

I am confused as to what is happening, as I am using .bdyseg.shp which is also in the file geodatabase for comparison and updating. Here is my code with the clip analysis in it;

# Import arcpy module
import os, sys, arcpy

InGDB = "C:\\Users\\mwisniewski\\My Documents\\NRM_Base_Data.gdb"
InShpFolder = "C:\\Users\\mwisniewski\\My Documents\\DPIPWE"

WSdict = {} # empty dictionary

# Local variables
clipRegion = "U:\\Boundary\\NRM_Regions\\DPIPWE_NRMRegionBoundary_Tas_Unk.gdb\\Cradle_Coast_NRM_Region"
outFc = ""

# Count the number of shapefiles clipped
count1 = 0
# Count the number of feature classes being copied
count2 = 0

# put all the names into a dictonary with a key of the name
for (path, dir, datasets) in arcpy.da.Walk(InGDB,datatype="FeatureClass"):
    for thisDS in datasets:
        WSdict[thisDS.lower()] = os.path.join(path,thisDS) # key is upper case to avoid case problems later

arcpy.env.workspace = InShpFolder
try:
 for thisFC in arcpy.ListFeatureClasses():
     thisName,thisExt = os.path.splitext(thisFC) # get the shapefile name without extension (base name)

     if WSdict.has_key(thisName.lower()):
         # the shapefile is found in the database..
        shpFull = os.path.join(InShpFolder,thisFC)# full path of shape file
        gdbFull = WSdict[thisName.lower()]        # get the path from the dictionary
        arcpy.AddMessage(shpFile)
        arcpy.Clip_analysis(InShpFolder, clipRegion, outFc, "", "", "")
        count1 += 1
        # Copy clipped shapefile into file geodatabase
        arcpy.FeatureClassToFeatureClass_conversion(outFc, InGDB + InShpFolder)
        count1 += 1
     else:
        arcpy.AddMessage("Error: Shapefile/s not clipped and copied to file geodatabase")

# Print message for shapefile name clipped
arcpy.AddMessage("\n" "The above {0} shapefiles have been clipped to: {1} and saved in {2}" "\n".format(count1, InShpFolder, InGDB))

# Return any errors
except:
    arcpy.AddMessage(arcpy.GetMessages())


texts:
I'm developing a Qgis plugin to use my own raster resampling method rather than what is provided in Qgis (cubic/bilinear).

And so, I've extended the QgsRasterResampler class. However, I'm having difficulty getting resample method to work.

Based on the C++ api documentation for resample method:
http://qgis.org/api/classQgsRasterResampler.html#aaffb31ed35ed12353e49458be27b06ba

The 'dstImage' arg is passed by ref and I'm not sure how return it to the calling function in python. As such, when the resampling method is applied, I get a blank image instead.


texts: Is it safe to delete files in the Staging > SharingProcesses directory in ArcGIS Pro?

texts:
In ArcGIS Pro, is it safe to delete files in the Staging > SharingProcesses directory after a webmap or web layer has been published online?

I would assume that it is but I have some important web maps that are used by many users on the portal and I don't want to risk experimenting.


texts: Setting projection to a geotiff file when it's shown as Unknown

texts:
I have a geotiff file that has no CRS in the metadata  but I know that the CRS of the raster is WGS 84 UTM 22 N.
How can I modify the geotiff file to have the CRS  WGS 84 UTM 22 N in the metadata.

Best Regards


texts: Find centroids of polyline and create point using Arcobject

texts:
I have one road feature class and one area feature class .What i have to do is i have to create a point at the centroids of the polylines of (road feature class) and find intersection of point with the area boundary.

How can i get the centroids of polylines in arcobject through vb.net code and create point with that (x,y) and then find intersection with area boundary????


texts: How do GPS receivers map out directions?

texts: Although I've never really dug into it, I've always assumed mapping out directions in a GPS unit is done with A* or Fibonacci-heaps at a higher level and Dijkstra's algorithm for smaller graphs. Are there any other algorithms being used that aren't immediately obvious?


texts: ArcGIS Runtime (WPF): Does Layer has FeatureTable property?

texts: According to the API reference, there is a ArcGISSceneLayer#FeatureTable property inherited from Layer. However, I can't find this property in my Visual Studios. Is the API reference wrong?


texts: How to create a geometry generator symbol using PyQGIS

texts:
I am trying to use PyQGIS to convert a Fill Symbol Layer into a geometry generator Symbol Layer in a rule-based renderer with 300 rules.
I am able to iterate over the rules and return each rule's symbol and then the respective (unique) Symbol Layer. I know there's a changeSymbolLayer()
What I am failing to do, is to create a QgsGeometryGeneratorSymbolLayer to overwrite the existing symbol layer. Until now, the code is something like this:
layer = iface.mapCanvas().currentLayer()
rules = layer.renderer().rootRule().children()

for rule in rules:
    label = rule.label()
    symbol_layer = rule.symbol().symbolLayer(0)
    print(symbol_layer)
    color = symbol_layer.color()
    generator = QgsGeometryGeneratorSymbolLayer()
    generator.setFillColor(color)
    generator.setGeometryExpression(expression)
    rule.symbol().changeSymbolLayer(0,generator)

In this case the QgsGeometryGeneratorSymbolLayer() fails by saying that I am not able to instantiate.
I am not understanding how to create the new symbol to apply on each rule.


texts: How to orthorectify pushbroom pixel locations with position and orientation data?

texts:
I'm trying to orthorectify/georectify data obtained from a hyperspectral line scanner via a UAV in Python. My system has both an RTK for highly accurate position data (long, lat), altitude readings (platform height), as well as an INS with the roll, pitch and yaw readings of the platform.
I have looked at multiple sources of mathematical/physical models to do this on a line-by-line basis, where my understanding would be that I can "project" the sensor line down to the ground (using the position and orientation data) and use raytracing with a DEM to determine an updated coordinate for where the pixel truly lies. However, many of these sources all use different notation and different variables, and some require GCPs while others do not.
The image below is an example from this website review on orthorectification, which uses collinearity equations. While it seems straight forward, it doesn't do much explaining on the "ortho-plane" coordinates (X, Y, Z) which makes me think this isn't directly applicable when considering multiple lines versus a single photo.

In terms of Python, I have also looked into existing Python packages like Rasterio to see if they have done similar things, but I have still been unable to find rectification processes that work directly  for line scanners. Building it in Python shouldn't be a huge issue, it's really just defining the mathematical formula for rotating and translating each line (prior to resampling) that I am struggling to wrap my head around.
Any help would be greatly appreciated


texts: Saving a ArcMap layers RGB values into the attribute table

texts:
I am wanting to take a symbolized layer in an mxd, and save it's symbol color (like a polygon fill) into its attribute table as an RGB value. My goal is to use this value in MapInfo to symbolize it's features via Discover.

I can't see any obvious way of using the ArcPy mapping module to read the particular values. Has anyone had a hack at building a simple modern solution?


texts: QGIS geolocate a DBF file with another point layer by matching field

texts:
I have a DBF table and a SHP with points. I want to geolocate the DBF file using the points layer. Both match in a field called "unificationstring".

I know how to do it in postgis but how can this be done in QGIS?


texts:
I am using developer edition of ArcGIS Qt SDK. What I am trying to do is, move camera from one point to another point in a timer. Somehow the camera does not move and is stuck on in its initial position.
// [WriteFile Name=CreateTerrainSurfaceFromLocalRaster, Category=Scenes]
// [Legal]
// Copyright 2019 Esri.

// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// http://www.apache.org/licenses/LICENSE-2.0

// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// [Legal]

import QtQuick 2.6
import Esri.ArcGISRuntime 100.8
import Esri.ArcGISExtras 1.1
import QtQuick.Controls 2.2

Rectangle {
    id: rootRectangle
    clip: true
    width: 800
    height: 600

    readonly property string montereyRasterElevationPath: System.userHomePath + "/ArcGIS/Runtime/Data/raster/MontereyElevation.dt2"
    property double pointX: 29.378955
    property double pointY: 41.159848
    property double pointZ: 500
    property Point test

    SceneView {
        id: sceneView
        anchors.fill: parent

        Scene {
            id: scene
            BasemapImagery {}

            Surface {
                RasterElevationSource{
                    // RasterElevationSource takes a list of raster files. Here we provide a list cotaining just a single file.
                    fileUrls: [montereyRasterElevationPath]
                    // Hook up success/error reporting for the Elevation Source load
                    onLoadStatusChanged: reportLoadStatus(loadStatus, loadError)

                }
                elevationExaggeration: 5.0
            }

        }

        // Once the scene view has loaded, apply the camera.
        GraphicsOverlay {
            id: graphicsOverlay
            AddMarker{
                id:marker
            }
        }
        Component.onCompleted: {
            setViewpointCamera(camera);
            marker.refreshSymbol(pointX,pointY,pointZ)
        }
    }
    function reportLoadStatus(loadStatus, loadError){
        if(loadStatus === Enums.LoadStatusLoaded){
            console.info("Loaded raster elevation source succesfully");
        }
        else if(loadStatus === Enums.LoadStatusFailedToLoad){
            console.warn("Error loading elevation source : ", loadError.message, loadError.additionalMessage);
        }
    }

    function createPoint(x,y){

        var point = ArcGISRuntimeEnvironment.createObject("Point", {
                                                              x: x,
                                                              y: y,
                                                              z: z,
                                                              spatialReference: SpatialReference.createWgs84()
                                                          });

        return point;
    }

    // Create the camera to be used as the scene view's viewpoint, looking at Monterey, California..
    Camera {
        id: camera
        location: Point {
            x: pointX // Longitude
            y: pointY // Latitude
            z: 300.0 // Altiude
            spatialReference: SpatialReference { wkid: 4326 }
        }
        heading: 10.0
        pitch: 70.0
        roll: 0.0
    }
    Timer {
        interval: 100; running: true; repeat: true
        onTriggered: {
            pointX=pointX+0.0001;
            pointY=pointY+0.0001;
            pointZ=pointZ+0.0001;
            marker.refreshSymbol(pointX,pointY,pointZ)
            test= createPoint(pointX,pointY)
            camera.moveTo(test)
        }
    }
}


texts: Projecting Lat Lng onto image coordinate system (reverse georeferencing?)

texts:
We have multiple illustrated maps. They are based on real places but some features are slightly distorted or exaggerated. They are also often slightly rotated.
We've already created several prototypes by georeferencing the illustrated maps in QGIS, exporting XYZ tiles as a directory and displaying them with Leaflet.
This works okay. However, it has the undesirable side effects of distorting and rotating the original illustrated map.
My ideal goal would be to render a Leaflet map in the image's original coordinate system. The Leaflet docs on Non-geographical maps seem interesting for doing this.
However, the difference is that these are maps that represent real-world coordinates (as opposed to video-game coordinates or something).
I also want to bring in other data with real-world lat/lng coordinates and display it on the illustrated map. Some of this data would be created by hand and could be created based on the image's coordinate system. However some of it would come from phone GPS sources (e.g. the user's current location).
The thing I'm looking for is a way to project lat/lng coords back into the image's coordinate system. This kind of feels like the opposite of georeferencing to me.
Before starting to use Leaflet or QGIS, we'd created a prototype that did pretty much this already. It took a bunch of control points, created a Triangulated Irregular Network (using turf) then looked up the triangle you were already in and calculated the new position based on your position in the current triangle. Although that worked, it was hacky and slow and was created before we'd even heard the term "georeferencing".
Does something like this already exist? Can anyone point me towards any examples that do this kind of projection using Leaflet?
Also, while some of the maps are distorted in complex ways, most of them could be georeferenced with only the corners if that makes it simpler.


texts: I have several raster tiles that are adjacent and non-overlapping, so they form a mosaic. I am trying to merge two adjacent rasters using PostGIS but am getting an error that they do not have the same alignment (pixel corner coordinates are not aligned). Can anyone tell me the correct way to combine the tiles into a single raster?


texts:
Is it possible to use MapBox GL JS to render vector-tiles in a Leaflet map, without the access token?

Even if the map object itself, does not come from Leaflet directly, but it's the wrapped MapBox version, is that possible?

My question is, does my data have to reside on MapBox's server? Or it can reside on my own server, being served with TileStache, for example?

All the Mapbox examples use the accessToken, but I just want to use the webgl rendering capabilities, and serve my data, or other base layer data (osm, google, etc).


texts: Georeferencing an NDVI (JPEG) image in QGIS

texts:
I am struggling to georeference a NDVI image.
The image is a JPEG, how would I go about this?
I am using QGIS


texts:
I have a problem. I'm currently trying to get the extent of a vector layer in latitude and longitude coordinates in python. The layer is currently in EPSG:3857, so .fullExtent() is not returning what I want.  I have tried reprojecting the layer using this:

self.layer= QgsVectorLayer(self.input, self.featureName, "ogr")
self.layer.setCoordinateSystem("crs=epsg:4326")

But the second line introduces a type error when loading the vector layer.

So what would be the best way of getting the bounds in latitude and longitude coordinates? It would be nice if I could get the coordinates from the pyqgis API, but I would be happy if there was another way tool that would do it, such as gdal.


texts: Listing available online WMS services (Weather, Land Data, Place Names)?

texts:
I'm Looking for WMS services freely available to the public to overlay on my application, specifically weather and AIS information.

Perhaps this is a good place to list some available online WMS services you have stumbled across?


texts: ArcMap Geoprocessing VBA code/tool - smooth polygon

texts:
I would like to execute the smoothpolygon tool within ArcMap VBA and I'm having some trouble defining the parameters syntax for the tool. I'm getting this error below:

ERROR 000622: Failed to execute (mooth Polygon). Parameters are not valid.
ERROR 000628: Cannot set input into parameter tolerance.

Here is my code:

Private Sub SmoothGeo_Click()

    Dim filePath As String
    filePath = "C:\Projects\OSU_Well_Buffer_Tool\USDE\USDE.mdb"
    Dim inputName As String
    Dim clipName As String
    inputName = "buffer9_smu"

    Dim pGP As Object
    Set pGP = CreateObject("esriGeoprocessing.GPDispatch.1")

    On Error GoTo EH
    pGP.Workspace = filePath
    ' note: the inputs can also be objects; however the output must be a string.
    pGP.smoothpolygon_management filePath + "\" + inputName, "C:\Projects\OSU_Well_Buffer_Tool\USDE\USDE.mdb\smooth", "1", "PAEK", "150"

    Exit Sub
EH:
    MsgBox pGP.GetMessages(), vbOKOnly, "Test"
End Sub


texts: Invalid SQL Error when using InsertCursor on table created using CreateTable

texts:
I am trying to write a list of items to a new table in Arc. When I attempt to add a new row to an Arc table, I get the following error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\__init__.py", line 1170, in InsertCursor
    return gp.insertCursor(dataset, spatial_reference)
  File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 369, in insertCursor
    self._gp.InsertCursor(*gp_fixargs(args, True)))
RuntimeError: ERROR 999999: Error executing function.
An invalid SQL statement was used.
An invalid SQL statement was used. [Table]
An invalid SQL statement was used. [SELECT * FROM Table]

I have narrowed down the error to using arcpy.InsertCursor() on a table created by using the arcpy.CreateTable_management() tool. I have recreated the problem using this code snippet:

tbl = arcpy.CreateTable_management (r'C:\GeoDatabase.gdb',"Table")
arcpy.AddField_management (tbl, "NewField", "FLOAT")
irows = arcpy.InsertCursor (tbl)

I get an error as soon as I try to execute the last line. This seems very simple, so I feel like I must be overlooking something. However, the InsertCursor tool does not use SQL statements, so the error is confusing to me. Can anyone help me understand how to bypass this error. If you know a better way to take a list of items and write them to a new table in Arc, I am all ears.

Thanks for the help.


texts:
I am using QGIS 3.10.2-A Coruña and PyQGIS via the internal Python Console (Plugins -> Python Console). Currently I am trying to import a shape file using the following command:

vlayer = QgsVectorLayer("path/file.shp", "layer_name", "ogr")

It seems to be loaded in the sense that I don't get any error message, I can return the column names (fields) and I can also see the attribute table via:

iface.showAttributeTable(vlayer)

However, the loaded .shp file is neither displayed in the list of layers nor is the corresponding map shown in the main window. I don't have that issue when loading files from a geoPackage or with point vector layers.

What do I have to do differently?


texts:
I produced a NetCDF file by combining two others using NCO (specifically ncbo command). I cannot open the output with QGIS though. Indeed, I get the following error msg:
"Download with "Protocol" source failed, do you want to try the "File" source type?".
Whatever answer I give the file is not opened.
I manage to get info about the file using gdalinfo. Below are the details:
Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: No UNIDATA NC_GLOBAL:Conventions attribute
Driver: HDF5/Hierarchical Data Format Release 5
Files: c_gls_output.nc
Size is 512, 512
Metadata:
  history=Wed Apr 14 16:15:16 2021: ncbo --op_typ=divide a.nc b.nc c_gls_output.nc
  lat_CLASS=DIMENSION_SCALE
  lat_NAME=lat
  lat__FillValue=nan
  lat__Netcdf4Dimid=0
  lon_CLASS=DIMENSION_SCALE
  lon_NAME=lon
  lon__FillValue=nan
  lon__Netcdf4Dimid=1
  NCO=4.7.2
  _NCProperties=version=1|netcdflibversion=4.6.0|hdf5libversion=1.10.0
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

Is there any workaround, may be changing some of the metadata, in order to be able to open the file?


texts: Is there an open source equivalent to the .mpk file format that ArcGIS can still open?


texts: QGIS custom coordinate reference system of a local site

texts:
I am creating a geo database for an archaeological excavation. We want to put data from a local system and from UTM35 in the same QGIS-project. I have Ground Control Points in both systems.
I setted the project to UTM35(EPSG 32635) and now want to create a custom coordinate reference system so that every layer in the local system could be converted on the fly to UTM35.
To calculate transformation parameters I used pyproj. I transformed the coordinates of the ground control points of both systems (local und UTM35) to ECEF WGS84 cartesian. To transform the local points to ECEF, I just assumed the cassini projection (Soldner). Second I calculated the 7 Parameters of the Helmert Transformation (with the code of Gabriel De Luca). I now have the parameters to transform between two ellipsoids.
translation (in m):
X: -195.06498913
Y: 495.68894225
Z:  28.56591918
rotation(counterclockwise):
X:   0.017997038245226384 (rad) ;  3712.1556066732  (arcsec)
Y:   0.00945795348286504  (rad) ;  1950.8429426372  (arcsec)
Z:   0.015533377928717144 (rad) ;  3203.9891888298  (arcsec)
scale factor:  0.9999905702146588  in ppm:  1.0000009999905701
I tried this:
+proj=cass +lat_0=37.53084417 +lon_0=27.27634361 +alpha=0 +x_0=2000 +y_0=2000 +ellps=WGS84 +towgs84=-195.06498913,495.68894225,28.56591918,3712.1556066732,1950.8429426372,3203.9891888298,1.0000009999905701 +units=m +no_defs
and this:
+proj=cass +lat_0=37.53084417 +lon_0=27.27634361 +alpha=0 +x_0=2000 +y_0=2000 +ellps=WGS84 +units=m +no_defs +x=-195.06498913 +y=495.68894225 +z=28.56591918+ s=1.0000009999905701 +rx=3712.1556066732 +ry=1950.8429426372 +rz=3203.9891888298
But the result is not as expected. If I import the local ground control points with it, they don't match with the UTM points. There is a distance of 1 up to 20 meters between them.
Everytime I click enter in the custom coordinate menu, QGIS changes the rotation angles to this:
+rx=0.0000049992 +ry=0.0000026272+ rz=0.0000043148
Why does QGIS do that?
Maybe I should use WKT instead.


texts: How to build SQL Expressions in ArcGIS 10.1 Python script tool dialog?

texts:
I've made a python tool that requires the user to establish which state to select within an input's table (the "Prepared Footprint"). However, once the user drops in the layer and tries to write the SQL statement to select the state which is contained in the state field of the Prepared Footprint, there are no available fields to select from.

Is there some option within the properties of the python script tool that will allow the Query Builder to populate the fields and attributes from the input layer like the native Select by Attributes in ArcMap? If not, how might I code this into my source script to feed fields and their attributes to the Query Builder of the "State" parameter?

(experience level with building tools: Beginner)

EDIT:
Here's a screenshot of the tool's parameters. When I select "State" and try to change the "Obtained from" paramter property nothing happens. There's no drop down, no cursor, nothing but a highlighted row for "Obtained from". Has anyone ever had this problem before?

EDIT:
Found the problem. A "Layer" data type cannot be used as the input for "Obtained from". Feature Layer works instead for this case.


texts: I am trying to join the information from a line to that of a point associated with that point. To explain a bit, it would be adding the information from a bus line to the associated bus stops with that line. The bus stops are not touching the line, so I was thinking about running a buffer to capture the bus stops, but some of them are very close to another and that would potentially add the wrong line information to the wrong point.


texts: Getting started with storing LoST and GML searching

texts:
I am working on a logging and search database application.

The application will receive event records containing XML LoST values (RFC5222 Location-To-Service-Translation). These records will be stored in a PostgreSQL database.

The application will later need to search for those event records using GML Polygon, Circle, ArcBand or Ellipse areas.

Can someone give me a pointer of where to get started with this?


texts: Connecting QGIS to table hosted on ArcGIS Server?

texts:
I have some aspatial tables hosted on ArcGIS Server (10.3.1 with Oracle 11g for the database) that I would like to load into QGIS (2.18).

The Add ArcGisFeatureServer dialog is only displaying the feature layers i.e. the tables are not listed.

Is there any functionality built-in to QGIS or available via a common plugin that allows me to load these tables?

As suggested by Francisco Raga in the comments below I have tried loading the table via the python API.

While this loads the spatial layer at index 0 of the Feature Server:

iface.addVectorLayer("url='https://example.com/.../FeatureServer/0' crs='EPSG:4326'", "title", "arcgisfeatureserver")

Both of these attempts to load the data table at index 1 produce only an "invalid layer" message:

iface.addVectorLayer("url='https://example.com/.../FeatureServer/1' crs='EPSG:4326'", "title", "arcgisfeatureserver")
iface.addVectorLayer("url='https://example.com/.../FeatureServer/1'", "title", "arcgisfeatureserver")


texts: Data for average monthly high/low by zipcode

texts:
I am looking for data that will provide average monthly temperatures by zip codes in the US. I require only a year or two worth of historical data, as I am only looking for general trends.

Wunderground used to provide summaries like this through their API, but they no longer provide free api keys.

Other APIs provide historical data per day, which I can aggregate for monthly averages. However, this would quickly hit the API limits. Please let me know if there is a data dump, or API I could hit that may have this information.


texts:
I saw this script in GEE website for masking clouds of Sentinel 2 imagery:
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

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

  return image.updateMask(mask).divide(10000);
}

However, I see other people changing the second part of the script to:
  var cloudBitMask = Math.pow(2, 10);
  var cirrusBitMask = Math.pow(2, 11);

or:
  var cloudBitMask = ee.Number(2).pow(10).int()
  var cirrusBitMask = ee.Number(2).pow(11).int()

I do not understand why they change the values 1 to 2. Is it something that has to be changed manually? I thought these bits has to remain as they are in the provided script.


texts: Is it possible to add a TMS layer to TileMill to compose a map, if yes, how is it possible to do it?


texts:
I am trying to extract a region from a .osm.pbf file I have and write the output as a .osm.pbf file. I'm not sure on the syntax however to achieve this result. I have continued from an example I found in another question but so far no luck. Here is what I have so far.

osmosis --rb file=germany-latest.osm.pbf -mw file=germany-small.osm.pbf --bb left=8.793491 right=11.5632 bottom=50.5139 top=52.5931

This tells me that the file germany-small.osm.pbf can not be found, which makes sense since it doesn't exist, it is the file I want to create. What would be the correct line to achieve what I want?


texts: Creating rings from existing polygons using QGIS?

texts:
I have been given a shapefile from a client that features pivots, yet they were recorded as 10 circles overlapping. The software my company uses requires that these polygons be rings that don't overlap. I have read about the clip and difference tools, but I think I'm going about it the wrong way.

Can anyone help me clip these polygons so they are rings and not overlapping circles?

I have attached a screenshot of what I mean.


texts: ESRI Javascript API Export to Shapefile and CSV

texts:
I am working on a Javascript Web Application using the ESRI Javascript API v3.3.  My application allows the user to edit and query data and also provides a table for attributes.

I would like the user to be able to export map data to a Shapefile or export the table data to a CSV.

Is this possible with Javascript along?  If not, my web application is running on Amazon's Elastic Beanstalk service with PHP.  I've never written anything in PHP, but can this be done with PHP?  How would I go about doing this in PHP?


texts:
How might I enable Quantum GIS to read ECW imagery?

At the moment I'm mostly concerned about this on Ubuntu 10.4 amd64, but I also use WinXP & Win7.


texts: CQL_FILTER that joins two feature types in WFS

texts:
I want to join two feature types in my WFS in a union. How can I achieve this using CQL? E.g.

http://localhost:8080/example/ows?
service=WFS&
version=2.0.0&
request=GetFeature&
srsName=epsg:4283&
typeNames=type1,type2&

This returns an error:

  Join query must specify a filter

I'm trying to append a cql_filter parameter, but the ones I try don't work.

What filter can I use to grab both types, in any format?


texts: Spatial Analyst Cell Statistics tool in Python gives inaccurate results

texts:
My question is similar to the Cell Statistics gives Wrong Answer question. In Python, I'm using the MINIMUM parameter to execute the cell statistics tool on several reclassified rasters ('0's and '1's) and I'm not getting accurate results.  However, when I use the tool in ArcMap I get what I'm looking for.  It seems as if it's only using two rasters (don't know which ones).  Here's the code I'm using.

arcpy.env.workspace = ("F:\\hdf\\hdf_test\\reclass")
outh10v03 = ("F:\\hdf\\hdf_test\\h10v03\\")

h10v03 = arcpy.ListRasters("*2000*h10v03*", "tif")

for reclass in h10v03:
    print reclass

arcpy.gp.CellStatistics_sa(h10v03, outh10v03 + "h10v03_2000" + ".tif", "MINIMUM", "DATA")

Any ideas as to why it's not using all rasters (I believe that's the case)?


texts: Field calculator 'if/then' statement, with 'and' not working

texts:
I'm trying to do some conditional formatting of an attribute table using Field Calculator. I tried using code below, but am not getting any returns that meet conditions specified (even though the conditions exist). What am I doing wrong? I am new to coding, so please assume I am ignorant to obvious mistakes.

My steps so far:

Created a new field "Raster_code" (this is the field where I am using field calculator)
Open field calculator, set parser to Python, enter the following pre-logic script:

def ifBlock(BGCLMTCLBL,STSRSLBLC1):
    if 'BGCLMTCLBL' == 'AT' and 'STSRSLBLC1' == '00':
        return 1
    else:
        return 0

Note: "BGCLMTCLBL" and "STSRSLBLC1" are two separate fields. I want to locate records containing a certain combination of values from these respective fields

Enter the following in the Raster_code block:

ifBlock(!BGCLMTCLBL!,!STSRSLBLC1!)

The result I get is 0 for all records-- in other words, the conditions were not met.

Here is a screenshot of the above-noted fields in my attribute table.


texts: Converting Longitude range from (0,360) to (-180,180) in Python or by using gdal commands

texts:
I have a raster downloaded.
Visualising it in QGIS.

I want to correctly project it so that its extent varies with longitudes from -180 to 180 instead of 0 to 360.

gdalinfo filename.nc outputs the following as the corner coordinates

Corner Coordinates:
Upper Left  (  -0.1250000,  90.1250000)
Lower Left  (  -0.1250000, -90.1250000)
Upper Right (     359.875,      90.125)
Lower Right (     359.875,     -90.125)
Center      ( 179.8750000,   0.0000000)

I am looking for a gdal command or python function to perform this task of rotating the image such that longitudes vary from -180 to 180.

Similar questions:

Verifying formula that will convert longitude 0-360 to -180 to 180? (Fortran solution)
https://stackoverflow.com/questions/25730625/how-to-convert-longitude-from-0-360-to-180-180 (R function of rotation)


texts: Drawing features with specific dimensions in ArcGIS

texts:
I am having some issues drawing transect lines in ArcGIS. I need to create some lines that have specific lengths (20m and 30m). My Map units are set to 'meters' and coordinate system is NAD UTM 10. When I hit CTRL+L and enter the desired  line length the line is far too long! I can verify this by using the measure/ruler tool.

I also need to draw some square polygons... and am having the same issue!

What is going on?


texts: Best function to smooth data layer in Google Earth Engine?

texts:
As you can see using the code below in Earth Engine, the resolution of the MODIS images are 1km squares that look quite obtuse.

// Create a geometry representing an export region.
var features = ee.Geometry.Rectangle([-85.1417893413635, 45.31413490213395, -86.125065708551, 44.65070625463291]);
Map.addLayer(features);
print("features",features);

// This function clips images to the ROI feature collection
var clipToCol = function(image){
  return image.clip(features);
};

var dataset = ee.ImageCollection('MODIS/006/MYD11A2')
                  .filter(ee.Filter.date('2018-05-01', '2018-06-01'))
                  .map(clipToCol);
print("dataset",dataset);
var landSurfaceTemperature = dataset.select('LST_Day_1km');
print("lST",landSurfaceTemperature);
var landSurfaceTemperatureVis = {
  min: 14375.62,
  max: 15300.38,
  palette: [
    '040274', '040281', '0502a3', '0502b8', '0502ce', '0502e6',
    '0602ff', '235cb1', '307ef3', '269db1', '30c8e2', '32d3ef',
    '3be285', '3ff38f', '86e26f', '3ae237', 'b5e22e', 'd6e21f',
    'fff705', 'ffd611', 'ffb613', 'ff8b13', 'ff6e08', 'ff500d',
    'ff0000', 'de0101', 'c21301', 'a71001', '911003'
  ],
};
Map.setCenter(-85.60371794450282,44.93590436363271, 10);
Map.addLayer(
    landSurfaceTemperature.filterDate('2013-04-10', '2013-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 13');
Map.addLayer(
    landSurfaceTemperature.filterDate('2014-04-10', '2014-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 14');
Map.addLayer(
    landSurfaceTemperature.filterDate('2015-04-10', '2015-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 15');
Map.addLayer(
    landSurfaceTemperature.filterDate('2016-04-10', '2016-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 16');
Map.addLayer(
    landSurfaceTemperature.filterDate('2017-04-10', '2017-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 17');
Map.addLayer(
    landSurfaceTemperature.filterDate('2018-04-10', '2018-05-10'), landSurfaceTemperatureVis,
    'Land Surface Temperature 18');
Map.addLayer(
    landSurfaceTemperature.mean(), landSurfaceTemperatureVis,
    'Land Surface Temperature Mean May');

// Where I am having trouble is coding the export feature that would allow for an overlay to be placed into Google Earth. This is what I have now:
print(landSurfaceTemperature.mean());
// Visualize image
var imageRGB = landSurfaceTemperature.mean().visualize(landSurfaceTemperatureVis);

print(imageRGB);
Map.addLayer(
  imageRGB)

// Define a boxcar or low-pass kernel.
var boxcar = ee.Kernel.square({
  radius: 200, units: 'pixels', magnitude: true
});

// Smooth the image by convolving with the boxcar kernel.
var smooth = imageRGB.convolve(boxcar);
Map.addLayer(smooth, imageRGB, 'smoothed');

// Export the image, specifying scale and region.

Export.image.toDrive({
  image: imageRGB,
  description: 'Overlay',
  scale: 5,
  region: features,
  maxPixels: 400000000,
});

I've been exploring convolution formulas for linear smoothing based on a per-pixel basis given a radius, however all of the examples use landsat images that have multiple bands. In the case of this code, how exactly would you structure it in a way that would achieve a moderate level of smoothing?

This is what I have worked with:

// Define a boxcar or low-pass kernel.
var boxcar = ee.Kernel.square({
  radius: 200, units: 'pixels', magnitude: true
});

// Smooth the image by convolving with the boxcar kernel.
var smooth = imageRGB.convolve(boxcar);
Map.addLayer(smooth, imageRGB, 'smoothed');

That returns the following error:

  Invalid argument specified for ee.Number(): true

What am I doing wrong? Is convolutions the correct method to use to smooth this data set?


texts: Is it possible to retrieve lost versioned data in SQL Server Management Studio?

texts:
One of my staff accidentally deleted a feature class she was creating in her version before I reconciled it.

I have grabbed a back up of the database from yesterday in SQL Management Studio. Is there a way to locate the table with the unreconciled data and restore it?

I'm using ArcGIS 10.1


texts: Convert data type from bytea to polygon psql

texts:
I was trying to convert data type of my column from bytea to polygon,

ALTER TABLE polygon ALTER COLUMN geom type polygon using
geom::polygon;

But I got an error saying cannot cast type bytea to polygon. Is there a solution to change the data type of the entire column or I can only use casting from ST functions.


texts: Plugin not taking acount a newly created field in a called Algorithm

texts:
I'm trying to write a qgis plugin and I'm having trouble when updating fields.

So here's the code:

# Create an Id column to be used as reference to the areas feature after the union
idnumber = 0
LYR_VECTORAREA.dataProvider().addAttributes([QgsField("id_temp", QVariant.Int)])
for feature in LYR_VECTORAREA.getFeatures():
     LYR_VECTORAREA.dataProvider().changeAttributeValues({ feature.id() : {fieldcounter: idnumber}})
     idnumber += 1
LYR_VECTORAREA.updateFields()
# Checking if the field has been created
for field in LYR_VECTORAREA.pendingFields():
     print field.name()

# Run the Union geoprocessing
processing.runalg("saga:union",LYR_VECTORAREA,LYR_VECTORDISTRIBUTION,False,VECTOROUTPUT)

The union processing doesn't take into account the newly added field. Even if the printed fields assess that the field has been created at this stage of the script (so, before processing the union), the new field never gets through to the output of the union.


texts: r.in.lidar Error: 'Raster map not found' in QGIS 2.18.14 GRASS

texts:
First of all I would like to say, that I have checked all the open threads (with almost the same title) about r.in.lidar errors and in general GRASS and it's "raster map not found" errors within QGIS 2.18.x.

The r.in.lidar in QGIS 2.18.14 and in GRASS 7.0.3 in Linux Mint and Windows 10 do the same error when running the tool:

Algorithm r.in.lidar - Creates a raster map from LAS LiDAR points using univariate statistics. starting...
g.proj -c proj4="+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=500000 +y_0=-5000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
g.region n=92000.0 s=91000.0 e=396000.0 w=395000.0 res=1
r.in.lidar input="/home/alen/Desktop/TM_395_91.laz" method=mean type=FCELL zscale="1" percent="100" resolution="1" -e -o output=output412a4072b5e749b8b5eaa46a1e1a183f --overwrite
g.region raster=output412a4072b5e749b8b5eaa46a1e1a183f
r.out.gdal --overwrite -c createopt="TFW=YES,COMPRESS=LZW" input=output412a4072b5e749b8b5eaa46a1e1a183f output="/home/alen/Desktop/obmocje.tif"
Cleaning up temporary files...
Starting GRASS GIS...
Executing '/home/alen/.qgis2//processing/grass7_batch_job.sh' ...
Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default
Projection information updated
ERROR: Unable to open file
ERROR: Raster map not found
ERROR: Raster map or group not found
Execution of '/home/alen/.qgis2//processing/grass7_batch_job.sh' finished.
Cleaning up temporary files...
Cleaning up temporary files...
Starting GRASS GIS...
Executing '/home/alen/.qgis2//processing/grass7_batch_job.sh' ...
ERROR: Raster map not found
ERROR: Raster map or group not found
Execution of '/home/alen/.qgis2//processing/grass7_batch_job.sh' finished.
Cleaning up temporary files..

For Windows 10 I tried for both the OSGeo4W install and the direct install. I tried it also with QGIS 2.16.x and 2.14.x for both the OS. In Linux Mint I apt-get autoremoved/autocleaned always after.
In Windows 10, if I install GRASS 7.2.2 separately and run the tool from it, it works well. I tried to reconnect QGIS with the separate GRASS install, but it does the same error again.

It worked until a few updates ago:
it works well in QGIS 2.18.5 with GRASS 7.2.0.

I noticed that also other r.* GRASS tools don't work anymore.

How can I make r.in.lidar work from QGIS?


texts: Converting LANDSAT-8 TIFF to matrix, data frame, or SpatialGridDataFrame formates

texts:
I need to correct the reflectance values of LANDSAT-8 data for steep terrain and find it too tedious and clunky in ArcGIS. The data I have is:
LANSDAT-8
DEM (ASTER)
Slope and aspect (Derived using DEM in ArcGIS)
Sun Elevation and Azimuth (Obtained from L8 MTL file)

All of which are clipped to a polygon of my study area.

Two methods I wanted to use were the cosine and minneart corrections. I found a package for R called "landsat" - http://cran.r-project.org/web/packages/landsat/index.html thatcontains a function named "topocorr" - http://www.inside-r.org/packages/cran/landsat/docs/topocorr which takes the uncorrected image, slope, aspect, azimuth and elevation then produces a corrected image.

I'm struggling to use it however as the input data type for the input image is "matrix, data frame, or SpatialGridDataFrame" formats.

How can I convert L8 data into this format?


texts: Sunshine hours duration calclulated from r.sun (GRASS Plug in) are same for diffrent date

texts: Currently i am trying to calculate sunshine hours duration from r.sun (GRASS Plug in ) from QGIS 2.12. I want to calculate sunshine duration for 21 st jun and 21st December for same elevation data set. After time consuming process i managed to get sunshine duration. But when i saw output for summer as well as winter both duration hours value are same. Is it related to DEM size or tools issue. because similar process i run with ArcGIS solar radiation tools and got very accurate results. if anybody worked on this please let me know


texts: How to auto-create a schema.ini file for a .csv?

texts:
I have been given some data in a csv file, which I need to import into my ArcSDE geodatabase. ArcGIS is having problems recognizing the field types, and is not importing some data. It is recognizing some fields as Integer, when it should have been Text.

I know that we are supposed to create a schema.ini file, which will contain the appropriate data type, like in this question.

My csv file has over 50 fields, so I was wondering if there was a way to auto-create this schema.ini file, and then edit it to the correct data types, instead of creating this file by hand, from scratch.


texts: Forcing routes to visit closest facility in ArcGIS

texts:
My dataset contains three parts: A) Incident location B) AED locations C) locations of available "first responders"(don't think of these like an ambulance or a fire department, but rather normal civilians)

Calculating the shortest path from C to A or C to B is doable with the Network Analyst tools, but is there a way for me to calculate two routes for each person in C? I.e. I'd like to have both a direct route, C to A as well as the route via the AED, C -> B -> A, where B in this case is the "closest facility"(closest AED) for the individual C.

I tried to work with barriers with negative value, but that didn't help much. Any pointers on how to create the two routes?


texts: Missing parameter INPUT_A in qgis:rastercalculator

texts:
I simply want to sum up two rasters with PyQGIS. To this end, I'm trying to use the qgis:rastercalculator but I encounter an error.
Here's the code:
import processing

output_layer = '/tmp/result.tif'
params = {
    'LAYERS' : [ "Input_DEM@1", "Input_Overlay@1" ],
    'EXPRESSION': '+',
    'OUTPUT' : output_layer
}
processing.runAndLoadResults('gdal:rastercalculator', params)

Which leads to this error:
Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 152, in runAndLoadResults
    context=context)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 153, in runAlgorithm
    raise QgsProcessingException(msg)
_core.QgsProcessingException: Unable to execute algorithm
Could not load source layer for INPUT_A: no value specified for parameter

But I cannot figure out such an INPUT_A field in the doc: https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasteranalysis.html#qgisrastercalculator
Did I miss something? What should INPUT_A be?
Also, it's not 100% clear how to pass the input rasters from the documentation. I may also be wrong in my way of doing it.
Info:
QGIS version 3.16.3-Hannover
QGIS code revision 94ac9f21b8
OS Version Ubuntu 18.04.5 LTS


texts: How do I interpret ESRI Licensing options to correctly build a quote that matches my requirements?

texts:
I don't know why ESRI licensing is so complicated, has so many options, and is so very poorly explained. The cynic in me says this is to force any potential customer to start a conversation with a sales rep. and ultimately prompt more purchases.

Yes I've seen the documentation but it's next to useless in any real world scenario that doesn't involve a county government.

So if anyone can provide assistance with my specific needs I'd be extremely grateful. Unfortunately this is unlikely to help many others as slight variations in requirements could tip the scale in another direction.

My requirements:

ArcGIS Server (10) and ArcSDE installed on separate machines
ArcSDE enterprise geodatabase on Oracle 10gR2
Up to 100 concurrent web users requesting a mix of cached and non-cached data from ArcGIS Server
Single map server machine (at first) with 16 CPU cores - SOCs and SOM both on one machine

From the documentation I've seen restrictions that prevent deployments on multiple machines, but it doesn't specify if this means SOCs and SOM distributed across machines or ArcGIS Server and SDE on separate machines.

There are also limitations on the number of concurrent users, but if I have 100 web users accessing pooled services does this mean 100 users or as many users as pooled services?

I would be really happy if anyone can tell me what I need in this scenario. I currently have access to ArcGIS Server Enterprise Standard and I need to know in advance if it will support the scenario described above. (I think) I've already figured out that the 'Workgroup' option is no good for me.


texts: Is there a way to expose Spherical Mercator Lat Lon in Google Map API?

texts:
The map tiles in the Google Maps are in Spherical Mercator projection (srid 900913), but in the API the methods report coordinates in Lat/Lon WGS84, e.g. map.getCenter().

Is there a way to return the Lat/Lon's in the native Spherical Mercator projection coordinates within the API?

For example, if I centered the map in London and the getCenter() returns 51.5, 0.0, is there a way the API can actually return it in Mercator meters? 6710219.08322, 0.0

(without having write an external algorithm to convert it)


texts: help identifying the right spatial analysis tool

texts:
Let's say I have a set of N towns and a road network that connects the towns. I think some towns might an influence on other towns because they meet some criterion. I have a binary indicator of this criterion (e.g., crit==1 or crit==0) for each town.

I want to construct a new variable that is a measure of how much each town might be influenced by it's neighbors, giving more weight to closer towns with a value of crit==1. I could define neighborhood as all towns within D km driving distance or R km radius.

What is the best tool in this case? Conceptually, I thought a spatial weights matrix might be useful, but I don't want a NxN matrix. I want a Nx1 matrix where the vector of values is a measure of how much each town might be influenced by its neighbors.

I created a toy example on GitHub.


texts: How to use ArcGIS 10.1 to make Python scripts for ArcGIS 9.3?

texts:
Recently I have been making scripts to run geoprocessing tools for a client who is running ArcGIS 9.3. However I only have access to ArcGIS 10.1. As a result I cannot use model builder in 10.1, and have to script everything with arcgisscripting (import arcgisscripting). I cannot find many arcgisscripting examples beyond the 9.3 help files.

My questions are: is there a way to reverse convert 10.1 (arcpy/modelbuilder) scripts to 9.3? Is there a place with good arcgisscripting .py examples? Has anyone else done this?


texts: Is it possible to open rasters as array in NumPy without using another library?

texts:
I know that it is possible to open a raster as an array in NumPy using GDAL, but I want to skip GDAL and use NumPy only, as it is cooler handling rasters with NumPy as matrices.
There is a similar question here: but the answer solutions involve using other libraries.

Can you post the code here if you are aware of a NumPy workaround?


texts: How to get base map ids from basemap gallery ArcGIS JS (i.e. topo)

texts:
Context: I am trying to simply create a 'Save Button' from selection of basemap in basemap gallery...
The only thing missing is I cannot find the 'well known string ids' from the basemaps in the basemap gallery for which I can use to reference later, i.e. streets, streets-vector, topo dark-gray, etc etc.
I am getting the objects responses with the below on any change of selection to the basemap gallery, I am seemingly getting EVERYTHING back but the 'well known string ids' in which I can save and use later to render that basemap. I am getting the title and property ID but the ids to change the map, i.e. topo are no where to be found within the object responses.
basemapGallery.watch('activeBasemap', function(newBasemap, property, object)
      console.log('new active basemap', {newBasemap, property, object});
.......... // also have tried pulling in the basemap obj here as well

i.e. I can programmatically, manually, change the basemap in my application or in the dev console with the below.
myapp.widgets.basemapManager.basemapGallery.activeBasemap = "streets"

The docs lead me to believe I can use portal item id to reference the basemaps but this did not work. I have tried to reference it every which way in accordance to the above line.
Below is a screenshot of my watch function after a selection has been made to the basemap gallery. No named string ids to change the map found, i.e. 'dark-grey' 'topo' 'streets' etc.


texts: User-Scaled Color Relief Map to Google Earth from DTM

texts:
I have a DTM file, and my goal is to load it up onto Google Earth as a color relief map, with the colors gradient at levels I chose. I have successfully done this with RAS MAP using HEC RAS software, but I have been unable to simply export to Google Earth.

Using QGIS I have been able to make a contour map that does the bare minimum of what I need, but the contours are far apart, and the kml file slows down Google Earth a painful amount.

I have scowered the internet, but have not yet found a viable solution.

Looking for original answer, or answers to assist the RAS Mapper or QGIS path I have already taken.


texts: Refreshing sketch geometry in OpenLayers 4?

texts:
I have an Openlayers 4 polygon and line draw interaction. There are also some construction tools to add vertices from a dialog.
These vertices are inserted using the geometryFunction on the draw interaction. The function modifies the coordinates and reapply them to the geometry:

Gui.Draw.Control.prototype._lineDrawGeometryFunction = function (coords,geom){
    if (!geom) {
        geom = new ol.geom.LineString(null);
    }
    if (this._coordinateForAddition != null)
    {
        var popped = coords.pop();
        var popped2 = coords.pop();
        coords.push(popped2);
        if (!(popped2[0] === this._coordinateForAddition[0] && popped2[1] === this._coordinateForAddition[1]))
            coords.push(this._coordinateForAddition);
        coords.push(popped);
        this._coordinateForAddition = null;
    }
    geom.setCoordinates(coords);
    return geom;
};

When drawing using the mouse cursor the sketch geometry is updated on the fly. But when constructed vertices are added the sketch is not updated until the cursor is moved in the map.

How do I force a refresh of the sketch geometry? Or is there a better approach to do this?

Note:

When calling ol.interaction.Draw#removeLastPoint() programatically the sketch is autoupdated. It uses the internal function

this.updateSketchFeatures_()

, but that is not available from the API. Is there a way to trigger this?

What I did try:

I have attempted the following - all without success:

    this.olSketchFeature.set("dummy", createGuid());
    this.olSketchFeature.setStyle(this.olSketchFeature.getStyle());
    this.olSketchFeature.getGeometry().setCoordinates(this.olSketchFeature.getGeometry().getCoordinates());
    this.olSketchFeature.changed();
    this.olSketchFeature.getGeometry().changed();

this.olSketchFeature is cached from the "drawstart" event.

I have also attempted to modify the olSketchFeature geometry, but the changes are removed as soon as the mouse cursor is moved.


texts: MapFish Print – Different background between first and next pages on pagination file

texts:
I’m currently using MapFish to print a paging on a same pdf, like this:

First page: A full map (which contains for example: a polygon as geometry feature + OSM default as background)

A pagination of this same map across next pages

I'm stuck on a problem: printing on the same pdf:

A full map (which contains for example: a polygon + ‘OSM default’ as background) on the first page

Across next pages, a pagination using a different background from first page (for example ‘OSM greyscale’)

I don’t know how to make it works, or if it’s possible... If anyone has a suggestion, or a solution?
I didn’t find any example of this on GitHub (https://github.com/mapfish/mapfish-print/tree/master/examples/src/test/resources/examples)


texts: Making Map so that anyone in ArcGIS Online organization can edit it?

texts: I am using arcgis online and I have created a map with layers and I have made it public with everyone in the organization. The only issue is that though the map is visible and everyone in the organization can see it if they log online. The individuals cannot see or edit the map in the arcgis collector app when they go out into the field. How do I make it so that everyone can edit the map so that it's visible in the arcgis collector.


texts:
I'm new to Mapbox tile-reduce library. I would like to process OSM tiles for Portugal.

I've download the mbtiles from Portugal.
I've installed tile-reduce with:

npm ini
npm install @mapbox/tile-reduce

I've downloaded the tile-reduce/examples/count/ example.
My index.js is:

'use strict';
var tileReduce = require('@mapbox/tile-reduce');
var path = require('path');
var numFeatures = 0;
tileReduce({
    zoom: 12,
    map: path.join(__dirname, '/count.js'),
    sources: [{name: 'portugal', mbtiles: path.join(__dirname, 'portugal.mbtiles'), raw: true}]
})
.on('reduce', function (num) {
    numFeatures += num;
})
.on('end', function () {
    console.log('Features total: %d', numFeatures);
});

My count.js is:

'use strict';

module.exports = function(data, tile, writeData, done) {
  var count = 0;
  if (data.portugal.roads) count += data.portugal.roads.length;
  if (data.portugal.buildings) count += data.portugal.buildings.length;
  done(null, count);
};

I was expecting to have the number of roads and buildings on the output, but the result is 0 features.

jgr@dusseldorf:/media/bonn/mapbox$ node index.js
Starting up 8 workers... Job started.
Processing tile coords from "portugal" source.
1876 tiles processed in 1s.
Features total: 0

What am I doing wrong?


texts:
I am trying to run the following query:

SELECT * FROM fences where ST_Intersects(geom, ST_GeomFromText('POINT(-73.990452 40.994184)', 4326))

When I run Explain Analyze I get the following:

"Index Scan using fencegeomindex1 on fences  (cost=0.41..8.68 rows=1 width=980) (actual time=0.231..0.353 rows=9 loops=1)"
"  Index Cond: (geom && '0101000020E61000008010C990637F52C06D1ADB6B417F4440'::geometry)"
"  Filter: _st_intersects(geom, '0101000020E61000008010C990637F52C06D1ADB6B417F4440'::geometry)"
"Planning time: 0.239 ms"
"Execution time: 0.393 ms"

When I run this query which is using a normal postgresql polygon:

SELECT * FROM fences WHERE poly @> '(-73.990452,40.994184)'

I get:

"Bitmap Heap Scan on public.fences  (cost=464.69..24194.22 rows=6229 width=980) (actual time=0.058..0.087 rows=9 loops=1)"
"  Recheck Cond: (fences.poly @> '((-73.990452,40.994184))'::polygon)"
"  Heap Blocks: exact=7"
"  Buffers: shared hit=14"
"  ->  Bitmap Index Scan on fencepolyindex1  (cost=0.00..463.13 rows=6229 width=0) (actual time=0.042..0.042 rows=9 loops=1)"
"        Index Cond: (fences.poly @> '((-73.990452,40.994184))'::polygon)"
"        Buffers: shared hit=7"
"Planning time: 0.094 ms"
"Execution time: 0.136 ms"

From a cost perspective the query using the postgis ST_Intersects appears to be much cheaper however it is taking much longer to plan and then to execute.  Any idea why this is happening.  I have Vacuumed and Analyzed the table.

Running the two different queries in our real application the one using PostGIS is netting us about a 1/3rd of the throughput as the plain old postgresql query.

One more note, we are currently Unioning multiple queries together to prevent extra network traffic not sure if this would affect how the queries are executed.

example:

SELECT * FROM fences where geom && ST_GeomFromText('POINT(-82.85966597206304 31.44405266605823)') UNION ALL
SELECT * FROM fences where geom && ST_GeomFromText('POINT(-82.85966597206304 31.44405266605823)')


texts: Issue with `-sds` naming convention in gdal_translate

texts:
I'm converting ~100 .hdf files to .tiff using gdal and getting some unfortunate naming conventions when using -sds

for %A in ("...\*.hdf") do gdal_translate -of GTiff -sds "%A" "%A.tiff"

However, it appears that it is first converting the files and then separating the subdatasets, therefore giving me file names like .tiff1, .tiff2, etc (completely useless).

Any suggestions on how to fix this within the same piece of code or with a simple second step?


texts: ESRI JS API - CIMSymbol with text from graphic attribute

texts:
The following sample shows a great use of the CIMSymbol in the ESRI JS API.

https://developers.arcgis.com/javascript/latest/sample-code/cim-symbols/index.html

Is it possible to set the textString property on the markergraphic to an attribute of the graphic?

https://github.com/Esri/cim-spec/blob/master/docs/v2/CIMSymbols.md#cimmarkergraphic

I'd like to place a number (attribute of the graphic) in the centre of the symbol and have tried "{$feature.numfld}" but this just displays the text string, not the expresssion value.


texts: QGIS 3.10.3: label callouts are not appearing

texts:
The 3.10 changelog mentions that label callouts are now available.  However I have not been able to get them to display, no matter what values and options I enter in the Layer Styling > Callouts window.  That is, no lines extend from labels to features.  Instead, I get the previous version labels, which are centered over the features.

I also tried with 3.12; no joy.

What am I missing?


texts: Using where clause to query date ranges in ArcPy

texts:
I would like to query a feature classes date range.  For instance, if there is data within a 2 day range then I would like to print true or do something else.  I thought arcpy.Exist would do the job, but I am missing something here.  I guess my arcpy.Exist is only seeing if the variable exists, how would I check to see if the actual where clause is True for my feature layer?

query = arcpy.MakeFeatureLayer_management(fc, "QueryLayer", """Date BETWEEN '2014-02-16 16:53:25' AND  '2014-02-17 18:53:25'""")

try:
    query

finally:
    if arcpy.Exists(query):
        print "true"
    else:
        print "no data"


texts: Crop and reproject TIFF to match bounds, pixel size of second TIFF

texts:
I have this working but I'm wondering if there's a better/faster way.

Currently, I use the following which creates an intermediate file (rasterb_cropped) from gdal.Warp() then creates the reprojected image with a call to gdal.ReprojectImage().

rastera_handle      = gdal.Open( rastera, gdalconst.GA_ReadOnly )
rastera_proj        = rastera_handle.GetProjection()
rastera_trans       = rastera_handle.GetGeoTransform()
rastera_band        = rastera_handle.GetRasterBand(1).DataType
rastera_xsize       = rastera_handle.RasterXSize
rastera_ysize       = rastera_handle.RasterYSize

# 1) crop rasterb to bounds
gdal.Warp( rasterb_cropped, rasterb_orig, outputBounds=[ lon1, lat1, lon2, lat2 ], cropToCutline=True )

rasterb_handle      = gdal.Open(rasterb_cropped, gdalconst.GA_ReadOnly)
inputProj           = rasterb_handle.GetProjection()
NoDataIn            = rasterb_handle.GetRasterBand(1).GetNoDataValue()

driver              = gdal.GetDriverByName('GTiff')
rasterb_final       = driver.Create( outputfile, rastera_xsize, rastera_ysize, 1, rastera_band )
rasterb_final.SetGeoTransform( rastera_trans )
rasterb_final.SetProjection( rastera_proj )
rasterb_final.GetRasterBand(1).SetNoDataValue( NoDataIn )

# 2) reproject rasterb to match pixel size (etc) of rastera
gdal.ReprojectImage( rasterb_handle, rasterb_final, inputProj, rastera_proj, gdalconst.GRA_Average )


texts: Mapinfo Runtime error when called from PowerBuilder

texts:
When I am trying to call MapInfo from my PowerBuilder application, with map info pro(full) version "Open table" command is working fine but with mapinfo runtime(proviewer) "Open table"  is crashing for the below syntax saying error calling external object.

MIOleObject.application.Do("Open Table 'C:\Jail\DataDir\GcuPerm.tab' as Zones Interactive")

Please suggest me what could be the problem?


texts: How to cut polygons to fit grid ArcGIS 10.3

texts:
I have a map that displays the distribution of different types of karst all over the US. The different types of soil are portrayed through polygons. I've also created a grid that divides the US, and the endgame here is to be able to analyze the fraction of karst within each grid box (i.e. use a command in arcpy or a geoprocessing tool that searches through all the gridboxes and gives back a table with all the labeled gridboxes with their respective fraction of karst). I have two issues here:

When the gridbox is created, I'm pretty sure the polygons aren't "cut" to conform to the lines of the grid. These are very complicated shapes that overlap the gridlines, so is there a way to cut these polygons so that when I go in and calculate the fraction of karst within each gridbox, the program will calculate JUST the karst info that lies within the confines of that gridbox?
What kind of command/tool would I be able to use to search through all the gridboxes and return the fraction of karst within each one?


texts: WFS request with Openlayers; blank request

texts:
I'm having trouble to add a WFS layer over an OSM map. At first I thought that my data had a problem, and I thus decided to adapt one of the example on the Openlayers website. I have no error returned, but the layer is not appearing... Instead, I see in Firebug that there is a WFS request made to the server, but it is empty

http://demo.opengeo.org/geoserver/wfs

No GET request or arguments whatsoever. I can correctly make the request and I get a xml file as a result (in eg Firefox). I thought it was maybe because of the a proxy and I set one up (the one in Openlayers FAQ), but I get the same thing.

Here is the code I use:

<!DOCTYPE html>

    <meta charset="utf-8">


texts: Installing Python or PyQGIS for QGIS-Android

texts:
Is it possible to install python-plugins on QGIS-Android? What are the steps for installing it? I am assuming that just like on Linux the package PyQGIS is needed?

I am using this version of QGIS-Android:

http://qgis.org/downloads/android/qgis-release.apk

System: Google Nexus 7 (2012)


texts: Repair data source for all layers simultaneously using QGIS

texts:
In QGIS 2.12 i opened a map after i changed it location in the hard disk.

Hence, all layers must be repair. I know the option to fix every layer separate, but i didn't find an option to "repair data source" for all layers simultaneously, like in ArcGIS.


texts: Geoserver SLD scale based view for multiple geometry styles

texts:
I have a FeatureStyle that draws line and adda an arraow to end point of the line. And I want to see this line and arrow by scale. So I specified a MaxScaleDenominator property. Now line is not appearing at scale but arrow is appearing. How can I solve this issue? Can I specify multiple MaxScaleDenominator?

   <FeatureTypeStyle>
        <Rule>
            <Name>Lines</Name>
            <ogc:Filter>
              <ogc:PropertyIsEqualTo>
                <ogc:PropertyName>line_type</ogc:PropertyName>
                <ogc:Literal>middle</ogc:Literal>
              </ogc:PropertyIsEqualTo>
            </ogc:Filter>
            <MaxScaleDenominator>6000</MaxScaleDenominator>
            <LineSymbolizer>
                <Stroke>
                    <CssParameter name="stroke">#aa0000</CssParameter>
                    <CssParameter name="stroke-width">1</CssParameter>
                </Stroke>
            </LineSymbolizer>
            <PointSymbolizer>
                <Geometry>
                    <ogc:Function name="endPoint">
                        <ogc:PropertyName>geom</ogc:PropertyName>
                    </ogc:Function>
                </Geometry>
                <Graphic>
                    <Mark>
                        <WellKnownName>shape://oarrow</WellKnownName>
                        <Fill>
                            <CssParameter name="fill">#fff433</CssParameter>
                            <CssParameter name="fill-opacity">0.5</CssParameter>
                        </Fill>
                        <Stroke>
                            <CssParameter name="stroke">#fff433</CssParameter>
                            <CssParameter name="stroke-width">1</CssParameter>
                        </Stroke>
                    </Mark>
                    <Size>15</Size>
                    <Rotation>
                        <ogc:Function name="endAngle">
                            <ogc:PropertyName>geom</ogc:PropertyName>
                        </ogc:Function>
                    </Rotation>
                </Graphic>
            </PointSymbolizer>
        </Rule>
    </FeatureTypeStyle>


texts:
I've been using ArcPy and Data Driven Pages for a while and now I'm trying to do the same thing in ArcGIS Pro but there doesn't seems to be as many resources.
I'm wondering if anyone has any working code to do a batch export of TIFF images in ArcGIS Pro?
This question has previously been asked as Batch exporting to tiff using ArcGIS Pro but unfortunately no useful answers were provided.
I've found some code, but I don't know what relevance it has.
import arcpy
import os
import sys
relpath = os.path.dirname(sys.argv[0])

p = arcpy.mp.ArcGISProject(relpath + "\\BigMapV3.aprx")
l = p.listLayouts()[0]

if not l.mapSeries is None:
    ms = l.mapSeries
    if ms.enabled:
        ms.exportToPDF(relpath + "\\Output\\Ex1_RI.pdf", "CURRENT", resolution=300)

"Traceback (most recent call last):
  File "<string>", line 7, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\_mp.py", line 332, in __init__
    self._arc_object = arcgisscripting._mapping.ArcGISProject(*gp_fixargs((aprx_path,), True))
OSError: \BigMapV3.aprx"

Ok I've made some advancement, though the error now occurs with the writing geotiff tags section, apparently the line of code encountered an unexpected object of "geoTIFF_tags"
import arcpy, os, sys
if not hasattr(sys, 'argv'):
    sys.argv  = ['']
relpath = "P:\\Big Map V3"

p = arcpy.mp.ArcGISProject(relpath + "\\BigMapV3.aprx")
l = p.listLayouts("Field Maps")[0]
print (l.mapSeries)
if not l.mapSeries is None:
  ms = l.mapSeries
  if ms.enabled:
    for pageNum in range(1, ms.pageCount + 1):
      ms.currentPageNumber = pageNum
      print("Exporting {0}".format(ms.pageRow.Id))
      pageName = ms.pageRow.Id
      l.exportToTIFF("P:\\Temp\\" + str(pageName) + "Ayylmao.tif", geoTIFF_tags = True)

Yet another update
I contacted ESRI support and they mentioned that the issue is that I was using ArcGIS Pro version 2.6. The ability to add in GeoTIFF tags was added in 2.7. (see https://pro.arcgis.com/en/pro-app/2.6/arcpy/mapping/layout-class.htm vs https://pro.arcgis.com/en/pro-app/latest/arcpy/mapping/layout-class.htm)
If anyone who is using 2.7 could test whether this works and continue from where I have left the code, that would be fantastic.


texts:
I am working with QGIS Composer 2.0 and I would like to add a second map.

What I would like to do is to have a main map that takes up the majority of the print area.  This map, Main Map, will have four layers active from QGIS Desktop.

The second map, Overview Map, is a small map identifying the overview of the area.  The overview map will be at a different scale than Main Map, and it will also only have one layer active from QGIS Desktop.

The issue I am having is that when I add a new map, the map areas are updated when I change the layers in QGIS Desktop.

QGIS 2.0 has a new function for Overviews, but the documentation is not overly clear on how to use it.

In ArcMap, you can set up different views, which allows for different data, scale, symbology, etc. ensuring that the different maps in the print layout are independent.

Is there a way to do something like this in QGIS Composer?

Two maps with different layers, from the same QGIS Desktop instance?