texts: How to create a PGDUMP from LVBAG format using ogr2ogr?

texts:
I am trying to figure out how to use this GDAL driver to create a SQL dump, using the ogr2ogr commandline. I have tried the following:
ogr2ogr -oo AUTOCORRECT_INVALID_DATA=YES --config PG_USE_COPY=YES -f PGDump /vsizip/BAGNLDL-08122020/9999STA08122020.zip

This runs without errors. However, it  creates a folder called 'PGDump' containing shapefiles, whereas I expected it to create a single SQL file. Upon adding the extra argument 'test.sql' (after the example of this driver) I get:
ogr2ogr -oo AUTOCORRECT_INVALID_DATA=YES --config PG_USE_COPY=YES -f PGDump test.sql /vsizip/BAGNLDL-08122020/9999STA08122020.zip

But then I receive the following error:
FAILURE: Unable to open datasource `test.sql' with the following drivers.

I've found different usage of quotation marks online, but adjusting them (e.g. "PGDump" instead of PGDump) does not change any of the results above. I'm running GDAL version 3.2.1.
Is there anything wrong with the syntax, or am I just misunderstanding how to use ogr2ogr?


texts: Error .ecw is not a supported raster data source

texts:
I have recently installed QGIS 2.0 in my computer which has Ubuntu 13.04 and when I try to use an .ecw file it says:

  Error, .ecw is not a supported raster data source.

I am new in Linux and in QGIS, so could you help me and explain what could I do step by step?


texts: How to reduce Polygon to Raster (Conversion) bit depth?

texts:
I have a polygon featureclass that I am trying to convert to raster format at a low bit depth (i.e. <8 bit unsigned).  The FC contains 50,000 polygons which have five values: 1,3,5,7,9.  I have set the value field to record the previously described values.  I have also tried adding a field populated with 1's and setting the Value Field to that.  The output raster is an 8-bit unsigned raster dataset with a NoData value of 255 every time.  I have tried changing the Environment setting > Raster Storage > NoData to every one of the options from NONE to PROMOTION.  I have even tried entering "0" as a NoData value in the environment setting.  I am pretty sure the tool does not honor that Environment setting.

How can I control the bit depth of the Polygon to Raster output?  Preferably I would like a 1, 2 or 4 bit raster output.  Work-arounds welcome too.

Edit:
Unfortunately, Copy Raster is not an option because the Polygon to Raster conversion is an intermediate step in a larger processing workflow.  Creating a new raster dataset would only add more time to the workflow.


texts: Building orientation and simplification using Opencarto Java library

texts:
I need to find out orientation of buildings and perform simplification of sides of buildings with a tolerance. I came to know from Stack Exchange that OpenCarto has useful libraries for geometrical functions. Basically i have stored data in PostGIS.

Is there any possibility to use Opencarto Java Library to perform this on buildings?

I think it is possible to connect Java with PostGIS through JDBC driver.


texts:
I have a question - where can I find a list of actual GPS position, for example - ships, airplanes etc.

I need it for my project for College - I will create mobile/desktop/web app, which will show on map current positions of this elements, their last road and save info about that elements in database.

Also, is there any website with informations about weather problems with lat/lng info?

Thanks in advance, I've tried to find something about it, but without result..


texts:
Is there any Standard data format for Radio Signal Coverage for GIS?

Something like GPX or KML (or other XML/JSON/etc-based) for routes, tracks, POI, etc. But for signal strength, base radio stations and other.

Or is there any Best Practice for processing, storage and interchange this type of data for Computer Systems?

My case

I want to (1) "impose" tracks from GPX/KML files "on the map" with signal coverage and then (2) select those sections of routes where the signal is lower than the specified value; I want to do it programmatically with popular programming or scripting lang in offline (not with any SaaS).

Observations

I know several websites that allow you to select some Country/Network and view result on top of OSM, Google and other maps. Do they use the same data format for this? Or, in each case, do they adapt different input data to the map elements (lines polygons, etc)? In this case, in what form they receive input data from MNOs? Or not from MNOs?


texts: Finding minimum bounding extent of given pixel value within raster?

texts:
I am wondering if there is a way to find the minimum bounding extent for a raster with a particular value. I clipped a raster from a global image and the extent is set as the global extent with a lot of NoData area. I would like to remove the NoData area from this raster and retain only the area containing the pixels of the particular value. How can I do this?

Here is my example: I would like to extract value = 1 (Blue area) and use the blue area's extent rather than the whole world for further processing.


texts: Getting maximum distance for point features in the same point group

texts:
For each group I want to find the maximum point-to-point distance (and which of the id1 features participate in the max).
For example, if I have the following attribute table:

FID Long    Lat id1 group
0   120.65627   23.649932   15  1
1   120.65677   23.650132   2   1
2   120.65887   23.652732   111 2
3   120.66057   23.654632   17  3
4   120.66167   23.655832   120 3
5   120.67133   23.65794    55  3

I would like to get an output table that looks like this:

Here the group=2 is omitted because there is only one point in the group, so no pt2pt distance to calculate (or it's fine if it is there and is zero). This is maximum distance for group=3 point combinations (others, 55.57 meters and 1,011 meters are less).
The Point Distance tool doesn't have any group-id capability nor any statistics (min, max, etc.) capablities.
I will be doing this in ArcPy, but I am not sure which tools or sequences of tools to use to get started.


texts: Problem with printing detailed navigation chart from PDF to HP Z5400PS

texts:
We have been given a PDF of a navigation chart from the Norwegian Hydrographic Service to print on our HP Z5400 PS wide format printer.  All is well with the printed chart except that the symbol for a submarine pipeline is not printing properly in places, and also the symbol for a restricted area. Someone has said it's to do with the layers in the PDF.

Can anyone help please?


texts: How to handle the "allow multiple selections" in "Value Relation Widget" results? (QGIS Form)

texts: Using the "allow multiple selections" in "Value Relation Widget" results in {value1,value2,...}. But arrays in QGIS are written like ['value1','value2',...]. How do I handle the multiple selection results? For example, if I would get the first entry in my list. With usual arrays I would use this syntax: array_get('array',0).


texts: Remove (or SetNull) areas of overlapping rasters using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.1.

I have one raster layer that is a site suitability layer for a county.  Then I have two feature classes (one for roads and one for hydrology) that I need to remove or just set to Null from the site suitability raster.  This is because the site suitability is for producing crops which obviously can't grow in water or on roads.

I have converted each feature class to a raster, but now I can't figure out how to get rid of the roads and water from the site suitability layer.  I've tried a few things with the Raster Calculator and Set Null tools but I haven't had any luck.  I'm probably doing something wrong because each time I end up with a result that is just whatever feature I am trying to remove from the site suitability layer.


texts: PyQGIS in standalone script, vector loading error

texts:
I am trying to use the QGIS functions in a stand alone Python script. I'm running Mountain Lion OS, 2.7 Python & the QQGISis 2.0.1 Dufour.

My main problem is that I cannot get started - I cannot load a simple shapefile in order to buffer it. After searching the online resources quite a bit, I actually found that many people had a problem with setting the prefix path, so I ran the QgsApplication.showSettings() to identify it.

The code I'm running now is:

 QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents", True)
    QgsApplication.initQgis()
    print QgsApplication.showSettings()
    #Loading a Layer
    layer = QgsVectorLayer('/path/to/shapefile_folder/test.shp', 'test', 'ogr')
    if not layer.isValid():
        print "The layer failed to load."
    #finally I d like to apply a buffer, sth like that maybe?
    buffer_Shp = layer.buffer(10, 5)

Here is the result I get:

 Application state:
    QGIS_PREFIX_PATH env var:
    Prefix:     /Applications/QGIS.app/Contents
    Plugin Path:        /Applications/QGIS.app/Contents/../PlugIns/qgis
    Package Data Path:  /Applications/QGIS.app/Contents/../Resources
    Active Theme Name:
    Active Theme Path:  :/images/themes//
    Default Theme Path: :/images/themes/default/
    SVG Search Paths:   /Applications/QGIS.app/Contents/../Resources/svg/
    User DB Path:   /Applications/QGIS.app/Contents/../Resources/resources/qgis.db

"The Layer Failed to load" (obviously)

bash_profile  (in case that sth is wrong there)

 /# Setting PATH for Python 2.7
    /# The orginal version is saved in .bash_profile.pysave
    /######export PYTHONPATH=/usr/local/lib/python:$PYTHONPATH
    export PATH=/usr/local/bin:$PATH
    export PATH=/usr/local/share/python:$PATH
    /#adding Qgis libs
    export DYLD_LIBRARY_PATH=/Applications/QGIS.app/Contents/MacOS/lib/:/Applications/QGIS.app/Contents/Frameworks/
    export PYTHONPATH=/Applications/QGIS.app/Contents/Resources/python/
    export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH
    /# Added by Canopy installer on 2013-10-06
    /# VIRTUAL_ENV_DISABLE_PROMPT can be set to '' to make bashprompt show that Canopy is active, otherwise 1
    /# VIRTUAL_ENV_DISABLE_PROMPT=1 source /Users/george/Library/Enthought/Canopy_64bit/User/bin/activate

What the heck am I missing?


texts: Renaming filenames to a standard filename format using Python

texts:
I am working within an ESRI ArcGIS environment. I have a folder full of .xyz files that are named to a specific old file format such as 2014516_MS_005.xyz I want to be able to loop through the whole folder and rename each file to a new format such as MS_e1614.xyz. As you can see there is alot of string displacement and new characters introduced.

Can you guys n gals point me to the right direction if there is already a resource similar to this question?

Here is what i have so far.

import os

survey_type = ['MB','SB']

#File path for directory of interest
folder = r'C:\Data Management\Implementation\xyz_data'

#List of all files in set directory
file_list = os.listdir(folder)
print file_list


texts: Why doesn't pgRouting return the best path?

texts:
Let the next portion of a graph:

When I use the function shortest_path between the points A and B, I got the blue path. Why does this happen?


texts:
I'm trying to automate a simple repetitive task on QGIS. I have some shapefiles with varying number of features in a layer and about 10 different attributes for each feature. I need to add the same string value ("negative") to one of these attributes, over and over again. Copy, paste works but since I want to sharpen my Python scripting skills I'd like to achieve this by creating a simple Python script.

I'm able to print values of attributes using something like this;

layer = iface.activeLayer()

for feature in layer.getFeatures():

    print("{acreage} acres in {name}".format(name=feature['fieldname'],acreage=feature['acreage']))

But I'm stuck at this point. How would I write "negative" to the attribute "acreage"? I tried set.Feature() but it just didn't work.


texts:
using ArcGIS, I have a layer of overlapping polygons and a raster grid underneath. The raster is not a continuous grid - the cells are spread out. I simply need to calculate the total area of the raster that falls within each feature of a polygon layer.

I've used Zonal Statistics and it works great for isolated polygons, however the polygons that overlap each other are not resulting in the correct pixel count for each individual polygon.


texts:
I use Proj 6.3 and I try to run Proj 'C' example from here. I tried it using Qt 5.14.2 on Ubuntu 20.04. I created a 'C' project from Qt Creator and tried the example. I get the errors:
g++  -o proj_c main.o
/usr/bin/ld: main.o: in function `main':
/home/edip/workspace/proj_c/main.c:42: undefined reference to `proj_context_create'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:44: undefined reference to `proj_create_crs_to_crs'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:57: undefined reference to `proj_normalize_for_visualization'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:62: undefined reference to `proj_destroy'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:68: undefined reference to `proj_coord'
/usr/bin/ld: /home/user/workspace//proj_c/main.c:71: undefined reference to `proj_trans'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:73: undefined reference to `proj_trans'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:77: undefined reference to `proj_destroy'
/usr/bin/ld: /home/user/workspace/proj_c/main.c:78: undefined reference to `proj_context_destroy'


texts: Populating multiple fields from string using Update Cursor?

texts:
I am trying to split up a comment field from a csv (or a dbf) into three different groups. The comment field is similar to this :

 Pa=68781/Rh=44.31/Ti=-7.56 AWEC Space Agency

I want to split it up so that three new fields (already created but empty) are populated by the Pa, Rh, and Ti values from each row.


texts: Read Pre defined Style in MapinfoTab File in Geotools

texts:
We have open a MapInfoTab File using OGRFeatureSource the file open on map successfully. But one question we have there is a function named GetStyleString in OGRLibrary of Ogr2Ogr but in our OGRSource the same method is not present. How we can override the same method.
And one more thing if we have a Style feature wise means each feature has a different style is there any option to open style feature wise without creating a Rule.


texts: Strange output from QGIS Raster Terrain Analysis plugin

texts:
This is closely related to several other questions on this site: here, here, here etc.
However, my issue is not resolved by following the answers provided.

I am using the Raster Terrain Analysis plugin in QGIS ver 2.14.3. I want to create a Slope layer based on SRTM imagery. I understand that confusion can arise because the lat-lon units are degrees, and the vertical units are in metres. I calculated the z-factor with the formula z = 1/(111320*cos(latitude*pi/180)) provided in this answer. Which gives 0.00000899928 (I'm using imagery which spans the equator so used latitude = 0).

However, I get a layer of either 0 or 90. It is very similar to the question here. Re-saving the layer with WGS84 before using the plugin does not help.

I am able to produce a sensible looking Slope layer if I convert the CRS to World Mollweide, but I would like to know why it isn't working with CRS in degrees? What am I doing wrong?


texts: Postgis Update using calculation from various fields

texts:
Trying to update one table with a calculation using data from another table.
I construct a subquery that relates acreagetot and acreage_wetland fields and use it to perform a calculation on UPDATE THEN clause.
Planner says 'Subquery has too many fields' SQL state: 42601.

WITH
type_2 AS (
SELECT v.parcel_id,v.acreagetot,w.acreage_wetland,v.geom
FROM mi.wl_int_sum AS w JOIN mi.bay_mi0 AS v ON w.parcel_id=v.parcel_id WHERE acreage_wetland>0)

UPDATE mi.bay_mi0
SET acreage_dry = (
    CASE
      WHEN parcel_id IN (TABLE type_2)
      THEN acreagetot - acreage_wetland
      ELSE NULL
    END
);

Both geometries uses same srid, 4326

And Postgresql versions are: -PostgreSQL 10.0, compiled by Visual C++ build 1800, 64-bit
-POSTGIS="2.4.0 r15853" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.2, released 2017/09/15" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER


texts: Failed to publish web layer - ArcGIS Pro and ArcGIS Portal

texts:
I'm attempting to create a new hosted Service Layer, from a FileGDB, on ArcGIS Portal, using the publication tools inside of ArcGIS Pro.
I have done this successfully before with other data.
This morning however, when attempting to create a new Hosted layer, it fails.
Using the standard "Share as Web Layer" tools in ArcGIS Pro, I configure a Feature layer type and share with everyone in organisation. None of this is new.
However this time, I am getting a generic error message "Failed to Publish Web layer" inside of ArcGIS Pro. When attempting to create the service layer from the ArcGIS Portal admin page, the process just hangs.
"Failed to Publish Web Layer" message is not very helpful at all and there are no further details in the application, so I'm just looking for some starting points to help narrow down the investigation. At this stage I've only found this article, which doesn't help unfortunately.
https://community.esri.com/thread/210811-arcgis-pro-failed-to-publish-web-layer-to-portal-for-arcgis


texts: Clip Raster by Mask Layer - Nodata output

texts:
I am trying to clip a raster around a polygon shapefile. I have managed to do it for a series of polygons with the same raster, however one polygon is returning a blank clipping with 'nodata' in each of the raster's bands.

I am working in EPSG:32755 coordinates and have checked previous posts but don't think it applies to my problem. The output raster is blank and has geometry outside the bounds of the above polygon.

Is the problem that the shapefile is too large to process? Each pixel contains 365 bands containing rainfall data from each day of the year for that area.

Edit: I have split the polygon into smaller pieces and the raster seems to be clipping to the mask fine, so how do I handle a larger polygon with more data?


texts: How to sum multi values for each pixel in PyQGIS?

texts:
I have a raster layer and a vector point layer. I calculated distance between center each pixel and each point based on @YoLecomte and @mgri, this link and xunilk, this link. For each pixel there are multi distance,I want to sum this multi distance for each pixel in raster layer and write them into a new raster file. In the last portion of this code to sum distances for each pixel, what is the error?

def pixel2coord(x, y):
    xp = (pixelWidth * x) + originX + (pixelWidth/2)
    yp = (pixelHeight * y) + originY + (pixelHeight /2)
    return QgsPoint(xp, yp)

pntLayer = QgsVectorLayer("/Data/points.shp","pointLayer",'ogr')

feats = [ feat for feat in pntLayer.getFeatures() ]

# Open tif file
ds = QgsRasterLayer("/Data/study.tif","Study")

pixelWidth = ds.rasterUnitsPerPixelX()
pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the layer
ext = ds.extent()

originX ,originY = (ext.xMinimum(),ext.yMinimum())
src_cols = ds.width()
src_rows = ds.height()
drive = gdal.GetDriverByName( 'GTiff' )
src_ds = gdal.Open("/Data/study.tif")
outBand = src_ds.GetRasterBand(1)

outData = numpy.zeros((src_cols, src_rows), numpy.float32)

pntRstList = []

for i in range(0, src_cols):
    for j in range(0, src_rows):
        rspnt = pixel2coord(i,j)
        pntRstList.append(rspnt)

sumP = []
for rpoint in pntRstList:
    for ft in feats:
        vgeometry = ft.geometry()
        rgeometry = QgsGeometry.fromPoint(rpoint)
        dist=vgeometry.distance(rgeometry) # get distance between cell center to point
    if dist < 200:
        sumP.append(dist)
        print sumP


texts: Putting Labels at the center or inside a polygon layer in QGIS

texts: How can I put labels at the center or at least inside a polygon layer using QGIS Lisboa (old version). I tried options provided on Layer Properties under placements but still my label just falls outside desired layer.


texts: Does ArcGIS Server 10.7.x always provide/support both Python 2.7 and 3.6, or sometimes the 2.7 one only?

texts:
Reading this ESRI support information from the Internet (select English language at the bottom of the page, other translations aren't always updated), it's not clear to me if a 10.7.x version of ArcGIS Server always comes with both 2.7 and 3.6 versions of Python language, or may sometimes come with the 2.7version only.

Is an ArcGIS Server component the same whether it is linked to an ArcGIS Desktop installation or an ArcGIS Pro one?


texts: Python error message with OpenLayers Plugin

texts:
I have just recently downloaded QGIS version 1.9.90-Alpha, Alpha, 5155176. I then went to install the OpenLayers plugin.  That successfully downloaded and installed, or so I thought, but when I went to use it I get an error message.

It says the TypeError is: sip.methoddescriptor object is not callable

Does anyone know what this means or how I can fix it?  I am using Python version 2.5.2

This is the error:

An error has occured while executing Python code:
Traceback (most recent call last):
  File "C:/Users/czuck/.qgis/python/plugins\openlayers_plugin\openlayers_plugin.py", line 48, in addLayer
    self.__plugin.addLayer(self)
  File "C:/Users/czuck/.qgis/python/plugins\openlayers_plugin\openlayers_plugin.py", line 197, in addLayer
    QgsMapLayerRegistry.instance().addMapLayers([layer])
TypeError: 'sip.methoddescriptor' object is not callable

Python version:
2.5.2 (r252:60911, Feb 21 2008, 13:11:45) [MSC v.1310 32 bit (Intel)]

QGIS version:
1.9.90-Alpha Alpha, 5155176

Python path:

['C:/OSGeo4W/apps/qgis-dev/./python', 'C:/Users/czuck/.qgis/python', 'C:/Users/czuck/.qgis/python/plugins', 'C:/OSGeo4W/apps/qgis-dev/./python/plugins', 'C:\\OSGeo4W\\bin\\python25.zip', 'C:\\OSGeo4W\\apps\\Python25\\DLLs', 'C:\\OSGeo4W\\apps\\Python25\\lib', 'C:\\OSGeo4W\\apps\\Python25\\lib\\plat-win', 'C:\\OSGeo4W\\apps\\Python25\\lib\\lib-tk', 'C:\\OSGeo4W\\apps\\qgis-dev\\bin', 'C:\\OSGeo4W\\apps\\Python25', 'C:\\OSGeo4W\\apps\\Python25\\lib\\site-packages', 'C:\\OSGeo4W\\apps\\qgis-dev\\python\\plugins\\fTools\\tools']


texts: Function argument for all values in R (raster extract, and velox raster extract)

texts:
I'm using the package velox to speed up raster extraction.

Simple question (and probably an easy one): what is the argument in R for fun=all values?

Context:
In the raster package the function extract fun=NULL is to return all pixel values. In velox the extract requires an argument "fun=". What is the function argument that will return all the cell values like the fun=NULL of the extract in the raster package?

Bonus points if you have a function that will take the mean while ignoring NA values (which is what I'm trying to do anyhow). Something along the lines of

meanvalue <- raster::extract(testraster, testpoly, fun=mean if x !=NA)

Id then use it on a velox raster (should work the same) like this rather than the above code:

meanvalue <- velox.testraster$extract(testpoly, fun=mean if x !=NA)

Fun=mean doesn't work because some of my polygons have NA values for pixels since they are not covered by the raster.


texts:
I need to use import ogr in my python window provided by ArcGIS 10.2, but i can't configure ogr and gdal on my system.

How can I fix this issue and import ogr in my python window provided by ArcGIS 10.2.


texts:
I am attempting to set up an Identify button that will only allow the identify functionality to be performed when the button is pressed.  I have gotten to the point that the identify only works when the button has been pushed the first time but I can't seem to turn the button off.  I have pasted a few bits of code below or you can follow the link to see the entire thing.

http://geoville.org/viewers/PopUp/

Funtion to activate Identify

function activateIdentify(){
    if (dijit.byId("tool_identify").checked) {
        dojo.connect(map, "onClick", executeIdentifyTask);
    }
  }

Button to Activate function

Identify

Function to toggle buttons

function toggleButtonIcon(tool) {

    //only the tools in the toolbar are dijit togglebuttons so can iterate thru them
    dijit.registry.byClass("dijit.form.ToggleButton").forEach(function(togbtn) {
      if (togbtn == tool) {
        togbtn.attr("checked", true);
      }
      else {
        togbtn.attr("checked", false);
      }
    });
  }  //end toggleButtonIcon

Thanks for your help.

-Mike


texts:
Is there any convenient way to change the opacity of a ol.style.Fill or ol.style.Stroke Object?

As far as I understand it, there are 2 ways to define a color with an opacity.

fill.color = [255,0,0,myOpcitay];

or

fill.color = 'rgba(255, 0,0,myOpacity)';

But what If I get a color attribute in hex format, like

var color = myStyle.getFill().getColor() // color #FF0000;

Do I really have to decompose the hex color definition in order to add the opacity? There must be a more user-friendly way to achieve this...


texts:
I am currently working on ECMWF Data for a curricular activity and for starters I'm working with aermr04 data (param: 4.210).
There are several things that bother me about the documentation in general, but this poses a relevant question: by typing dataset.variables the following result is presented:
int16 aermr04(time, level, latitude, longitude)
scale_factor: 2.1468113041240392e-13
add_offset: 7.034241919092827e-09
_FillValue: -32767
missing_value: -32767
units: kg kg**-1
long_name: Dust Aerosol (0.03 - 0.55 um) Mixing Ratio
unlimited dimensions: time
current shape = (2, 44, 201, 237)
filling on}
How do I work with those scale_factor and add_offset?
The data that I have studied presents a aermr04 of max value: 0.00103 and min value: 1.3157193486989465e-11


texts: I tried to connect uDig with Sql server 2008 but could not do it the following error poped that it missed imageio ext and jdbc driver for sql server. I added it to jre/ext/lib (both of them). Now the error is not present but when I click next and click finish there is no change.What is the Error in here?


texts:
I have a large amount of .TIF which come in EPSG 31467. I want to transform them into EPSG 3044. My gdalwarp command throws me various errors, no matter what I try:

gdalwarp - s_srs EPSG:31467 -t_srs EPSG:3044 infile.tif outfile.tif

throws:
ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for .\infile.tif. There is no affine transformation and no GCPs. Specify transformation option SRC_METHOD=NO_GEOTRANSFORM to bypass this check.

If I try then:

 gdalwarp - s_srs EPSG:31467 -t_srs EPSG:3044 -to SRC_METHOD=NO_GEOTRANSFORM

I get completely unusable files. Any idea how to transform those .tif's ?


texts:
I have a connection with basicAuthentication, the getCapabilities and getFeatureStore return me correct data but when i try to commit a transaction I get the next exception: "HTTP response code:401 for url ...."

Here's my code:

public boolean anadirIncidencia(String getCapabilities, String capa,
        SimpleFeature incidencia, String username, String password) throws Exception {
    boolean b = false;
    Map connectionParameters = new HashMap();
connectionParameters.put("WFSDataStoreFactory:GET_CAPABILITIES_URL",
                getCapabilities);
connectionParameters.put("WFSDataStoreFactory:USERNAME",username);
connectionParameters.put("WFSDataStoreFactory:PASSWORD",password);

    System.out.println("Connexion Parameters: " + connectionParameters);
    DataStore data = DataStoreFinder.getDataStore(connectionParameters);

    SimpleFeatureSource source = data.getFeatureSource(capa);
    List<SimpleFeature> features = new ArrayList<SimpleFeature>();
    features.add(incidencia);

    Transaction transaction = new DefaultTransaction("add");

    SimpleFeatureStore featureStore = (SimpleFeatureStore) source;
    featureStore.setTransaction(transaction);
    System.out.println("FeatureStore:" + featureStore.getDataStore().getNames());

    try {

        System.out.println("ADD: addFeatures");
        List<FeatureId> added = featureStore.addFeatures(DataUtilities
                .collection(features));
        System.out.println("added: " + added.toString());
        transaction.commit();
        b = true;

    } catch (Exception problem) {
        problem.printStackTrace();
        transaction.rollback();
        b = false;
        error = problem.getLocalizedMessage();
        System.out.println("ADD: rollback ----->" + error);
    } finally {
        transaction.close();
    }
    return b;
}

And the error happens in transaction.commit(); and then after commit fails, transaction.rollback() throws the same error.

Someone can help me?


texts: Using an IF statement in Field Calculator

texts:
I'm trying to use an if statement in the field calculator:
This is what the code looks like in Python:

a= 2
FacilityID= !FACILITYID!
if a%2==0:
   MH_ID_Calc=FacilityID[0:6]
else:
   MH_ID_Calc=FacilityID[7:14]

print MH_ID_Calc

  06M244

so I tried to run that in ArcGIS field calculator:

def myCalc(FID,FACID):
if FID%2==0:
return FACID[0:6]
else:
return FACID[7:14]

myCalc( !FID! !FACILITYID!)

But it fails each time.


texts:
I'm trying to create a curve between two points in PostGIS that meets some target distance. So my input is point_a, point_b, target_distance, and my desired output is a curve between those two points of approximately target_distance. I'd want to find curves both in the 'left' and 'right' direction.

I have hacked together the below using another StackOverflow post as a source. This produces a curve like I want (and I can negate the st_length to invert the curve), and I can modify the figure I divide ST_Length by to increase/decrease the size of the curve.

However, I can't immediately figure out the maths that will allow me to get from a target distance to this divisor. Any suggestions?

Of course, if there is a simpler way to achieve this I am all ears!

\pset pager off
\set point_a ST_SetSRID(ST_Point(0.0098956, 51.4831476),4326)
\set point_b ST_SetSRID(ST_Point(0.0000000, 51.52156),4326)
\set line_ab ST_MakeLine(:point_a, :point_b)

-- Direct distance: 4328.567m
SELECT ST_Length(:line_ab::geography);

-- Produces a slight curve of distance 4377.397m
SELECT ST_Length(
        ST_CurveToLine('CIRCULARSTRING(' ||
                       st_x(st_startpoint(:line_ab)) || ' ' || st_y(st_startpoint(:line_ab)) || ', ' ||
                       st_x(st_centroid(ST_OffsetCurve(:line_ab, st_length(:line_ab)/10, 'quad_segs=4 join=bevel'))) || ' ' ||
                       st_y(st_centroid(ST_OffsetCurve(:line_ab, st_length(:line_ab)/10, 'quad_segs=4 join=bevel'))) || ', ' ||
                       st_x(st_endpoint(:line_ab)) || ' ' ||  st_y(st_endpoint(:line_ab)) || ')')::geography
);

-- Produces a large curve of distance 9615.056m
SELECT ST_Length(
        ST_CurveToLine('CIRCULARSTRING(' ||
                       st_x(st_startpoint(:line_ab)) || ' ' || st_y(st_startpoint(:line_ab)) || ', ' ||
                       st_x(st_centroid(ST_OffsetCurve(:line_ab, st_length(:line_ab)/1, 'quad_segs=4 join=bevel'))) || ' ' ||
                       st_y(st_centroid(ST_OffsetCurve(:line_ab, st_length(:line_ab)/1, 'quad_segs=4 join=bevel'))) || ', ' ||
                       st_x(st_endpoint(:line_ab)) || ' ' ||  st_y(st_endpoint(:line_ab)) || ')')::geography
);


texts: User-friendly maintenance of the PostGIS geodatabase

texts:
We have a PostGIS geodatabase that is populated by several people. The base is gradually growing, and it becomes unclear to users what data we have and how to find it. In fact, only the person who loaded the data into the database knows about them. This can lead to data duplication.
The only thing that comes to mind is maintaining an Excel table describing the data.
Are there better ways to find and catalog information in a PostGIS geodatabase?


texts: Is there a .NET ArcGIS Server REST Wrapper library that supports 10.6 or later?

texts:
I am looking for a .NET wrapper for the ArcGIS Server REST API.  I have found several on nuget and github but they haven't been updated in several years.
I didn't see a a REST wrapper on https://github.com/Esri?q=&type=&language=c%23 but maybe they rolled that in to something else?
What are all the cool kids using in 2021?
The ArcGIS Server version I am consuming is 10.6


texts: Maxent error 'java.lang.NumberFormatException: For input string:"Z^,,,"'

texts:
I converted a shp file into asc file in QGIS. Then I tried running the Maxent with that enviromental layer and came across that error message:

java.lang.NumberFormatException: For input string:“Z^,,,”

Can anyone help me?


texts: Openlayers-6 .fit(extent, {options}) not working with options

texts:
I have upgrade openlayers in IONIC V3 app from 4.6.5 to 6.1.0. Now map.getView().fit(extent, {duration:1000}) stop working, when i remove the options object it works fine.

Am i doing something wrong or anything else is needed.?

var extent;
extent = s_layer.getSource().getFeatures()[0].getGeometry().getExtent();
map.getView().fit(extent);
//this works fine.

map.getView().fit(extent, {duration: 1000});
//Doesn't fit the extent,just a slight movement on map.

I want to achieve a slight animation on map, when fit the extent of feature. Also tried animate() method of view class but not able to get the desired results.


texts: PostGIS- slow spatial UPDATE queries on 50mil rows

texts:
PostgreSQL 11.8, AWS hosted RDS instance, SSD, 8GB RAM, 3.0 GHz Intel Scalable Processor. Plenty of storage space.
I'm finding running updates on a large table using spatial queries is taking a very long time (>24 hours). How can I speed it up? I looked into HOT updates, but I can't use that as I have an indexed column (the geometry column, which is essential for fast spatial processing)- as I understand it, HOT updates don't work if any column is indexed.
Here is my process:
I create a new table of 50 million points as a copy of another table (I don't want to alter the original).
CREATE TABLE schema.table AS
SELECT id, geom
FROM schema.orig_table;

I add about 40 new columns which will store '1's where a spatial condition is true. This doesn't have to be a 1 though, I may need the number to be a 2,3,4 etc.
ALTER TABLE schema.table
ADD COLUMN col1 integer DEFAULT 0,
ADD COLUMN col2 integer DEFAULT 0,
ADD COLUMN col3 integer DEFAULT 0, etc etc

I set the id column as the primary key.
ALTER TABLE schema.table ADD PRIMARY KEY (id);

I assign a geometry type to the geometry column.
ALTER TABLE schema.table
ALTER COLUMN geom TYPE geometry(POINT,27700) USING ST_SetSRID(geom,27700);

I add a spatial index on the geometry column.
CREATE INDEX table_gix ON schema.table USING GIST(geom);

From here on, I run multiple UPDATE queries on the table, setting one of the 40 columns with the value 1 (or 2/3/4... etc I decide on a query by query basis) where the point intersects another table (I have a load of geospatial tables representing various polygons, points, line features. All spatially indexed and working correctly).
UPDATE schema.table a
SET col1 = 1
FROM schema.polygons b
WHERE ST_INTERSECTS(a.geom, b.geom);

...

UPDATE schema.table a
SET col29 = 1
FROM schema.another_polygon_table b
WHERE ST_INTERSECTS(a.geom, b.geom);

etc etc.... x200

Another 200 of these UPDATE queries occur. All the other geospatial tables used in the query are clean, fast, and without geometry errors. And most of the time, these ST_INTERSECTS queries don't find anything, so the value in the table stays at 0. The problem is, some of these queries can take up to 1 hour each. So 200 of these? Causes a headache.
The time taken for the whole script to run can take 24 hours. But putting a LIMIT 1000000 as a test, shows that it can run quickly on smaller sample sizes. It also shows that the overall script time grows exponentially with the amount of rows updated in the table.
How can I speed up this whole process? What's the most efficient way of running 100s of spatial UPDATE queries on a table of 50 million rows?


texts:
I'm having a problem that is universal to any shapefile that I'm trying to convert into a DXF with QGIS.

I've used the "Export DXF" option in the toolbar as well as the "Save As" option when right clicking the layer and neither will work.  If I use the "Save As" the data will successfully export into a DXF file.  The problem is when I try to open that DXF in AutoCAD.  All of the information is there, but it's at a microscopic scale (1100 feet turned into something like .003).  If I measure the same distance in the original shapefile in QGIS, it's spot on. I also just tried opening the DXF back into QGIS and that lines up perfectly as well.

I also tried scaling the drawing in CAD (not perfectly precise but just to get something workable) and it wouldn't scale correctly.  If I choose two points and then specify a length it doesn't come anywhere near matching that length. It's almost like it's a difference between the "Length" measurement in the properties menu and the units being used to perform the actions (although I don't know of any units that scale to 1:5,000.

Does anyone know why this happens or how to work around it?


texts: Remove sliver gaps between polygons with R

texts:
Is there a way to eliminate the small "slivers" between polygons using R? The desired solution would create a new SpatialPolygonsDataFrame where the shared boundaries between polygons are coincident. I'm specifically interested in a solution that uses R, rather than ArcMap or QGIS.

I would also be interested to hear an explanation as to why these gaps exist in the first place.

Here's a reproducible example of the spatial data I'm working with:

library(rgdal)
library(sp)
library(tigris)
library(magrittr)
library(leaflet)
library(gplots)

# This project will use WGS 84 projected coordinate system
crs_proj <- CRS("+init=epsg:4326")

# These are the FIPS codes of the specific block groups in my study area
sel <- c("530330079005", "530330079001", "530330079004",
         "530330085002", "530330085003", "530330086003",
         "530330087003", "530330085001", "530330090001",
         "530330091001", "530330091002", "530330092001",
         "530330092002", "530330086001", "530330090002",
         "530330086002", "530330079003", "530330079002",
         "530330087002", "530330087001")

# Create polygons
polygons <- tigris::block_groups(state = "WA",county = "King") %>%
        .[.@data$GEOID %in% sel,] %>%
        spTransform(CRSobj = crs_proj)

# Map the result
leaflet() %>%
        addProviderTiles("CartoDB.Positron") %>%
        addPolygons(data = polygons,
                    stroke = F,
                    fillColor = col2hex("red"), fillOpacity = 1)

As you can see in the above screenshot, there are small gaps between the census block group polygons. These location of these gaps is rendered differently depending on the zoom level, but there are always some gaps visible.

Can anyone recommend an R function (or a combination of functions) to programmatically adjust the polygons to eliminate these gaps?


texts: What does the -d flag 'Dither' in GRASS GIS mean? For example in r.composite


texts: determine whether a circle intersects a county boundary in postGIS

texts:
I'm using PostGIS 2.0.1 and need to determine whether a circle (center + radius) intersects a county boundary (polygon).

The query I'm using is this:

SELECT DISTINCT(u.id) FROM siteareas_poly sa, userlocations ul,
            clienttypes ct INNER JOIN users_clienttypes uct on ct.id = uct.clienttype_id
            INNER JOIN users u on uct.user_id = u.id
            WHERE sa.id = 'concord' and ct.id = 2
            AND ul.geog_lnglat is not null
            AND  ST_Intersects(st_exteriorring(sa.boundary), ST_BUFFER(ul.geog_lnglat,ul.radius_km * 1000)::geometry)

where the ul.geog_lnglat is a

geography(Point,4326)

and the sa.boundary is

geometry

I've tried a few variations on this and always gotten no results, when I expect to get at least one, as I have a row in userlocations that lies in that county (by zip lookup).

Am I making any glaring errors in my approach?


texts: Can we hide/show a Torque layer and/or set SQL?

texts:
Is it possible to select a Torque layer that is part of a visualization?

I am using CreateVis and creating button and layer selector so that the user can choose what appears on the map. However, the Torque layer do not act like the others sublayer and I cannot select it.

Is there a way to set SQL so that only part of the data are presented when a button is clicked for a Torque Layer, while having other "normal" layers on the map?


texts: Editing and manipulating tables with ArcPy?

texts:
I’m new to ArcMap and was wondering if anybody could help me out. I have a shapefile with counties and wells in that county. In the table I’m trying edit it so that for each county calculate the total numbers of wells that fall into that county, and update the county table with the total count. I also want to calculate how many wells there are per area unit, and update the table with that value as well.

import arcpy
featureClass = "C:/Lesson5_Data/Colorado_cities.shp" fieldList = "*"

with arcpy.da.SearchCursor(featureClass, fieldList) as cursor:
for row in cursor:
print row[4

population = 10000

whereClause = '"POP2007" >= ' + str(population


texts: TileStache:TileStache.Core.KnownUnknown: Couldn't get a layer from data source!

texts:
I want to use tilestache to get geojson files. I hava run tilestache-server.py -c tilestache.cfg ,and it works well.

However,I want to get tiles from data stored in postgis. Maybe it will need to run tilestache-seed.py. So I run`

tilestache-seed.py -c ~/TileStache-1.49.8/cfg1 -b 32.000 114.000 31.000 115.000 -l postgis 12  -e geojson

then I got

Traceback (most recent call last):

File "/usr/local/bin/tilestache-seed.py", line 331, in
    mimetype, content = getTile(layer, coord, extension, options.ignore_cached)
  File "/usr/local/lib/python2.7/dist-packages/TileStache/init.py", line 64, in getTile
    status_code, headers, body = layer.getTileResponse(coord, extension, ignore_cached)
  File "/usr/local/lib/python2.7/dist-packages/TileStache/Core.py", line 414, in getTileResponse
    tile = self.render(coord, format)
  File "/usr/local/lib/python2.7/dist-packages/TileStache/Core.py", line 500, in render
    tile = provider.renderTile(width, height, srs, coord)
  File "/usr/local/lib/python2.7/dist-packages/TileStache/Vector/init.py", line 583, in renderTile
    layer, ds = _open_layer(self.driver, self.parameters, self.layer.config.dirpath)
  File "/usr/local/lib/python2.7/dist-packages/TileStache/Vector/init.py", line 453, in _open_layer
    raise KnownUnknown('Couldn\'t get a layer from data source %s' % source_name)
TileStache.Core.KnownUnknown: Couldn't get a layer from data source PG:dbname='osmgis' user='postgres' host='127.0.0.1' password='xxxxxx'port='5432'

here is my cfg1

{

"cache":
  {
    "name": "Disk",
    "path": "~/TileStache-1.49.8/stache",
    "umask": "0000"
  },
  "layers":
  {
    "postgis":
     {
       "provider": {
                "name": "vector",
                 "driver": "PostgreSQL",
                "parameters":
                         {
                                "dbname": "osmgis",
                                 "user": "postgres",
                                "password":"xxxxxx",
                                "host": "127.0.0.1",
                                "port": 5432,
                                "query": "SELECT ST_AsGeojson(geom) from wuhan_osm_roads"}}
                        }

I am a beginner ,and I don't know how to solve the problem.

by the way, if i delete "12" behind "postgis" in code,it works but output nothing.


texts: Adding rows in attribute table together using ArcGIS Pro

texts:
I am using a parcel map and I have a layer with a ton of polygons and those polygons have different names. I combined all the names into categories in the symbology tab but when I go to the attribute table for that layer the ones I combined are no longer combined like in the symbology. I don't want to merge these polygons because I want to be able to count how many polygons are in the categories but I want the rows in the attribute table combined into the groups I put them in in the symbology.
I am using ArcGIS Pro


texts: GeoServer Raster Data Sources: satellite data's contrast issue

texts:
I am trying to use my own satellite image as the basemap for a webapp, instead of OSM, Google etc imageries. In order to do that, I made a Raster Data Sources store inside my GeoServer, which is pointing towards a .tif file (having .prj and .tfw files along-with). I am on Ubuntu 14.04. Everything is working fine except the following issue.

As can be seen in the attached images, visualization of satellite image in QGIS is giving a good contrast and brightness. But importing the very same image inside GeoServer is producing low contrast/bright image, plus seams, which can be seen at the bottom section. Although I am using Tiling:Tiled in the GeoServer, it shouldn't create seams. Zoomin and out are also changing the contrast and brightness of the GeoServer image.

Kindly explain the reason and remedy of this unexpected effect.


texts:
I have a problem with the installation of QGIS 1.8.0.

After finishing the installation two icons show on the desktop: Qgis desktop and Qgis browser.

When I click on these, it shows the same message, that there is no file msvcr70.dll in the computer and that it can help to reinstall the program (it didn't).

Thanks for any ideas.


texts: Example of QGIS for landscape management and maintenance plans?

texts:
I'm writing about landscape management and maintenance plans, and that most of them are in a Word file with perhaps a table of tasks, as opposed to a smart document such as a QGIS file with a linked database of tasks.

A plan would include a list of all tasks that might need to be undertaken, and then a table showing when each task should take place.

Writing it in Word means there is no way to filter tasks, such as only showing tasks that are required in week 10, or all tasks for a specific location.

My suspicion is that this would be far better in QGIS, showing areas in plan and a filterable range of tasks associated with each area that can be exported as worksheets for site staff.

Does anyone know of an example of something similar? It could be habitat management or landscape management.


texts: I have a project with a large number of raster files, but it is becoming very slow. I wonder which is the best way to handle all this. I'm simply georeferencing these files creating .wld file (geoTiff is better?) and than I open them as raster. Am I too naive?


texts:
I have multiple locations saved in WGS84 - it looks like e.g. latitude: 12345678 longitude: 87654321. I tried to search it on the map this way: 12.345678 87.654321, but it doesn't match the place it should be - it is several meters away. I just need to get proper GPS coordinates. Am I doing it the right way? Is this actually WGS84? Or do I just need to convert it somehow to use it? Thank you.

Update:

Real example:

latitude: 50075482
logitude: 14430068

It should be Legerova - Rumunská street in Prague, Czech republic. I can get this place by dividing the numbers with a dot, but it doesn't look like WGS84.


texts: Is it possible to put Geoserver behind a proxy and access the WMS/WFS using OpenLayers

texts:
I just start to work with Geoserver and Openlayers and one of the main concerns of the project now is do not expose the geoserver web interface outside of the internal network. No access to www.mydomain.com/geoserver/

However the application that is available in the Internet needs to get data (wms/wfs) from the GeoServer. This all through OpenLayers.

Is this possible?


texts: Nearest Neighbor calculation in PostGIS?

texts:
I'm learning how to use PostGIS and spatial databases for analysis. What I am trying to do is perform a calculation to get the distance for the nearest polygon in a file, using the edge calculation, rather than vertices.

Using this answer from Paul Ramsey to Finding minimum edge to edge distance of polygons using ArcGIS Desktop? which is a similar question:

  CREATE TABLE mytable_distances AS
  SELECT a.id, b.id,
  ST_Distance(a.geom::geography,
  b.geom::geography) as distance FROM
  mytable a, mytable b;

I am attempting to apply it to my spatial database. I don't understand the structure of this query though. I think CREATE TABLE mytable_distances AS creates a table to store the result but after this part I'm lost. Are a and b column names? If so, why would I specify two columns to calculate this?

My table is called TestArea and I have experimented with some basic queries successfully:

SELECT
  "TestArea".hgt
FROM
  public."TestArea"
WHERE
  "TestArea".area > 100

The structure of the database in PGAdmin III is as follows, with my table called TestArea. I'm not sure what the nearest neighbor calculation should look like using my column headers (all of these objects are polygons).


texts: CartoDB json to HTML Table using Javascript or Jquery

texts:
I am trying to display the dataset inside the cartoDB as HTML table. When i give request to the postgre database  by SQL API, it is returning as the format (json) below in the snapshot. I want the returning dataset to be displayed as HTML table in the browser. I had tried with the below code but it is displaying table for var myList=[{"name" : "abc", "age" : 50},
                       {"age" : "25", "hobby" : "swimming"},
                       {"name" : "xyz", "hobby" : "programming"}];

But its not displaying for the below format.

And code i had tried is

<script type="text/javascript">
        data1 = "transit_point_region_1";
        json_link = 'http://development.localhost.lan:8080/api/v2/sql?q=SELECT * FROM ' + data1;
        var data = $.getJSON(json_link);
        // Builds the HTML Table out of data.
        function buildHtmlTable() {
            var columns = addAllColumnHeaders(data);
            for (var i = 0 ; i < data.length ; i++) {
                var row$ = $('<tr/>');
                for (var colIndex = 0 ; colIndex < columns.length ; colIndex++) {
                    var cellValue = data[i][columns[colIndex]];
                    if (cellValue === null) { cellValue = ""; }
                    row$.append($('<td/>').html(cellValue));
                }
                $("#excelDataTable").append(row$);
            }
        }

    function addAllColumnHeaders(data). {
        var columnSet = [];
        var headerTr$ = $('<tr/>');
        for (var i = 0 ; i < data.length ; i++) {
            var rowHash = data[i];
            for (var key in rowHash) {
                if ($.inArray(key, columnSet) === -1){
                    columnSet.push(key);
                    headerTr$.append($('<th/>').html(key));
                }
            }
        }
        $("#excelDataTable").append(headerTr$);
        return columnSet;
    }
</script>
</head>
<body onLoad="buildHtmlTable()">
    <table id="excelDataTable" border="1">
        </table>


texts: I import map from file osm.pbf using osm2po. Why the low speed function call pgr_dijkstra in database Postgresql?

texts:
Good day!
I import map from file russia.osm.pbf using osm2po. Why the low speed function call pgr_dijkstra in database Postgresql (>2 minute)?

My steps:

java -Xmx6g -jar osm2po-core-5.0.0-signed.jar prefix=pgr tileSize=45x45,1,c  RU.osm.pbf

set PSQL="c:\Program Files\PostgreSQL\9.4\bin\psql.exe"
set PGPORT=5432
set PGHOST=localhost
set PGPASSWORD=1
cd pgr
%PSQL% -U postgres -d pgrouter -q -f "pgr_2po_4pgr.sql"

The database created a table with an index on the field source and target:

CREATE TABLE pgr_2po_4pgr
(
id integer NOT NULL,
osm_id bigint,
osm_name character varying,
osm_meta character varying,
osm_source_id bigint,
osm_target_id bigint,
clazz integer,
flags integer,
source integer,
target integer,
km double precision,
kmh integer,
cost double precision,
reverse_cost double precision,
x1 double precision,
y1 double precision,
x2 double precision,
y2 double precision,
geom_way geometry(LineString,4326),
CONSTRAINT pkey_pgr_2po_4pgr PRIMARY KEY (id)
)

My query:

SELECT * FROM pgr_dijkstra(
    'SELECT id, source, target, cost FROM pgr_2po_4pgr',
    30837, 682499,false, false
)

Query execution time is more than 2 minutes. Why?

This query is very fast (~1 second):

http://localhost:8888/Osm2poService?cmd=fr&source=30837&target=682499&findShortestPath=true&ignoreRestrictions=false&ignoreOneWays=false&routerClassId=0&heuristicFactor=0.0&maxCost=0.0&format=geojson


texts: Creating child entries in QGIS form with Spatialite many-to-many relationship

texts:
I am trying to reproduce the example from Introducing many-to-many (N-M) relations using Spatialite. Creating the tables and the relations works fine, but when I try to create a child element in a form in QGIS, the junction table cannot be saved because of a failed foreign key constraint. What am I doing wrong?
This is what I did:

Create a new Spatialite database using spatialite_gui (SQLite version 3.22.0, SpatiaLite version 4.3.0a)

Create three tables, airports, airlines, and a junction table airports_airlines using the following SQL statements:
CREATE TABLE airports (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    airport_name text NOT NULL
);

CREATE TABLE airlines (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    airline_name
    text NOT NULL
);

CREATE TABLE airports_airlines (
    id INTEGER PRIMARY KEY AUTOINCREMENT,
    airport_fk INTEGER
        REFERENCES airports (id)
        ON DELETE CASCADE
        ON UPDATE CASCADE
        DEFERRABLE INITIALLY DEFERRED,
    airline_fk INTEGER
        REFERENCES airlines (id)
        ON DELETE CASCADE
        ON UPDATE CASCADE
        DEFERRABLE INITIALLY DEFERRED
);

Load the three tables in QGIS (version 3.12.0).

In Project Properties | Data Sources activate the option "Automatically create transaction groups where possible"

In Project Properties | Relations, set the relations with the junction table as shown below:

In the layer properties, set the cardinality of the relation to airports (id) and airlines (id), respectively, add the relations widget to the drag-and-drop form and remove the id text widget, as shown below:

Try to create a new airline and an associated airport at the same time. This seems to work fine at first, but trying to save the layers fails with the following error messages
CRITICAL    Layer airports_airlines : SQLite error: FOREIGN KEY constraint failed
          SQL: INSERT INTO "airports_airlines"( "airport_fk","airline_fk") VALUES ( ?,?)
CRITICAL    Layer airports_airlines : SQLite error: FOREIGN KEY constraint failed
          SQL: INSERT INTO "airports_airlines"( "airport_fk","airline_fk") VALUES ( ?,?)


texts: Preserving uniqueness for overlapping points when converting from points to raster?

texts: I'm trying to optimize outflow points based on lowest elevation within some tolerance. Using the SnapPourPoint function which does exactly this, but outputs a raster file. I've been converting this raster back to points using RasterToPoint however it creates an issue when multiple points that are within close proximity and snap to the same raster cell. The resulting RasterToPoint output will only preserve one of the points. Points sharing the same minimum elevation are to be expected, but I would like to preserve a unique identifier for each point, is there any way to achieve this?


texts: Is it possible to create a tooltip with a map service?

texts:
I am trying to modify the ESRI Parks Finder application to include a tooltip that shows the name of a park when you mouse hover over the icon. It is a map service. Most of the documentation I have seen shows doing this with either a feature service or graphic. In fact, there are other areas in the application that utilize tooltips but they are based on graphics not the actual symbology of the data.

Any ideas or help on how I can go about this?


texts:
I have a rectangular polygon I want to create another polygon that is 200 metres bigger in each direction from the original.

Is there a way to create a square buffer around the original polygon or another method to do this?


texts:
Does the specification for DTED 0 - 2 also apply to 3 - 5? I know the finest resolution supported by spec MIL-PRF-89020B is 0.1 arc-seconds which I think would barely cover up through DTED level 4, right?  I'm guessing then that my code for reading DTED 0 - 2 using just the header information should still work for the later DTED levels but I have no way of proving this.

Anyone know where to download sample DTED Levels 3 and up? I could use that as well.

Thanks!
Matthew Hanna


texts: Spatial Join a point layer to polygon layer using st_join in order to create a choropleth map using R

texts:
I have three spatial layers, 1 USA boundaries and 2 point layers for Arby's and Hardee's locations across USA. Goal of the assignment is to create 2 choropleth maps for density of these fast food joints shown by state.
What I am able to achieve so far;

#importing layers

usa <- read_sf("/Users/username/Desktop/r-intro/data/seminarClass/usa.shp")
arbys <- sf::st_read("/Users/username/Desktop/r-intro/data/seminarClass/restaurants.gdb", layer="arbys")
hardees <- sf::st_read("/Users/username/Desktop/r-intro/data/seminarClass/restaurants.gdb", layer="hardees")

#crs check

usacrs <- st_crs(usa) # "ESRI",102003
usacrs

arbyscrs <- st_crs(arbys) # "EPSG",4269
hardeescrs <- st_crs(hardees) # "EPSG",4269
usageom <- (usa$geometry) #usa geometry

#crs matching

arbys_usa <- arbys %>% st_set_crs(st_crs(usa))
st_crs(arbys_usa)
hardees_usa <- hardees %>% st_set_crs(st_crs(usa))
st_crs(hardees_usa)

# spatial join - THIS IS WHERE I AM HAVING TROUBLE

arbyJoin <- usa %>% st_join(arbys_usa) %>% group_by(STATE_FIPS) %>%
  summarise(...)
hardeeJoin <- usa %>% st_join(hardees_usa) %>% group_by(STATE_FIPS) %>%
  summarise(...)

All coding except the part in spatial join is done by me and I understand what's going on in the code. I emailed my instructor in order to ask for help and he send me this line shown in the spatial join section of my code;
arbyJoin <- usa %>% st_join(arbys_usa) %>% group_by(STATE_FIPS) %>%
  summarise(...)

I don't understand exactly what's going on in this part therefore I can't figure out what to do with the summarise part.
I have tried
summarise(arbyJoin, sum=(STATE_NAME))

and
summarise(arbyJoin, sum=(geometry))

I don't necessarily get an error but I don't understand if they work or not as well, since I don't understand the code.
Can anyone shed some light on this line of the code and point me the right direction on how to proceed with this summarise function?
Ultimate result I'm trying to achieve is to create a choropleth map using these joint layers.


texts: Projection mismatch when plotting precipitation with ggplot

texts:
I'm trying to plot accumulated precipitation in Brazil for a given period, see code below:

Load packages

library(raincpc)
library(SDMTools)
library(raster)
library(ggplot2)
library(rgdal)
library(maptools)
library(rgeos)
library(mapproj)

Use cpcrain package to obtain data

cpc_get_rawdata(2014, 1, 01, 2014, 1, 08)

rain1 <- cpc_read_rawdata(2014, 1, 01)
rain2 <- cpc_read_rawdata(2014, 1, 02)
rain3 <- cpc_read_rawdata(2014, 1, 03)
rain4 <- cpc_read_rawdata(2014, 1, 04)
rain5 <- cpc_read_rawdata(2014, 1, 05)

rain_tot <- rain1 + rain2 + rain3 + rain4 + rain5

Zoom in Brazil

lon_vals <- seq(275, 330     , 0.5)
lat_vals <- seq(-60, 20      , 0.5)
reg_box <- expand.grid(lon = lon_vals, lat = lat_vals,
stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)
reg_box$rain <- extract.data(reg_box, rain_tot)
reg_box$rain_chunks <- cut(reg_box$rain, breaks = c(0, 80, 160, 240, 320), include.lowest = TRUE)

Use ggplot for visualization

gfx_gg <-
  ggplot(data = reg_box) +
              geom_raster(aes(lon, lat, fill = rain_chunks)) +
              scale_fill_manual(values = c("lightgrey", "red", "green", "blue")) +
              theme(axis.text = element_blank()) + # , axis.ticks = element_blank()
              labs( x = NULL, y = NULL, fill = " Precipitation (mm)") +
              ggtitle("Rainfall over South America for selected period")

which, as expected, yields:

Add state boundaries

brazil1 <- getData("GADM", country = "BRA", level = 1)

gfx_gg +
  geom_polygon(data = brazil1, aes(x = long, y = lat, group = group),
               color = "black", show_guide = FALSE, fill=NA) +
  labs(x = "x", y = "y") +
  theme_bw()

But I'm getting an unexpected error:

How can I fix this (projection?) error?


texts:
I'm going to show KML file using Leaflet's plugins named Omnivore KML (https://www.mapbox.com/mapbox.js/example/v1.0.0/omnivore-kml/)

I've read the tutorial, implement the source code, and it generate a map like this :

But now, what i want to do is changing it's layer/terrain style to be like this :

How can i achive that?

I've read that there's another Leaflet Plugin named Leaflet Bing Layer (https://github.com/digidem/leaflet-bing-layer) to show map's layer/terrain like the image above.

Can i combine Omnivore and Leaflet Bing Layer to achieve this?

Or i can just add some code into my leaflet javascript code to change it's layer/terrain style?


texts: CartoCSS : Style a road using a road_label rule

texts:
The parisian ring road (and probably others) usually looks a bit off (especially at low zoom levels), because in OSM it has been created as two independant 'main' roads, instead of one 'highway'.

What I could do is apply my highway style to all the OSM ways of the ring road, using a [osm_id=xxx] selector, but it's going to be a bit tedious.

All the OSM ways of the ring road have a 'road_label' tag with the same 'Boulevard Périphérique Intérieur' value. In this case would it be possible to apply a style on the road, using a road_label property? If so, how?

Thanks !


texts: How to perform a Cost Path analysis in ArcGIS 10.1?

texts:
Is it possible to model the easiest pathway through a landscape between many points (rather than just two)?

I want to carry out a cost-path analysis,  I would like to display the easiest pathways through the entire study region (not just one route but all)- by taking into account slope and proximity of sites to rivers or valley floors.

I have a DEM raster (obtained from UK ordnance survey - Landform-panorama). River systems were in the form of a polygon shapefile, which I have converted in to a raster layer. The lat/lon coordinates of the 10 sites where added to the DEM raster as a feature class.

I was wondering what would be the best way to achieve this? and whether anyone could explain in layman's terms how to do it?

Any advice would be massively appreciated!


texts: Getting the coordinates for tripoints where three counties intersect

texts:
I'm trying to make a simplified triangulation of a US state (say Ohio) based on the counties. The way I accomplish this is to treat all boundaries, no matter how crooked, as straight lines whose vertices are all the points where three counties meet. How would I find the requisite coordinates from a shapefile for further processing?

RN I plan on getting the coordinates for use in Matlab for a personal algorithm that is based on Delaunay Triangulation, and I've attempted to use shapefiles describing said counties on Google Earth Engine and ArcGIS without much luck, unless I'm taking things in the wrong direction.

EDIT: The file is from http://ogrip-geohio.opendata.arcgis.com/datasets/7212b60cc5fc4949869995813628182c_12

and I am trying to use QGIS 2.18.14 to get the individual points. Is there an easy way to highlight and pinpoint all the places where at least 3 counties touch?


texts: Creating Sub-Blocks Within Existing Polygon using ArcGIS Desktop

texts:
I am new to ArcGIS Desktop.

I'm trying to create some sub-blocks within a block. The idea is to divide the block to sub-blocks which follow natural features like river/stream, block's boundary, topography etc. And, the neighbouring sub-blocks should share the boundary. I don't how to make it right, snapping is not the right choice for me when tracing the stream and border especially on location with lot of curve.

Could someone suggest the best way to make it done?

I posted some photos just to share what I try to do.


texts: QGIS field calculator expression: selecting min/max values within field groups

texts:
I am new to using QGIS (previously used Arc) and trying to get my head round the expressions and Python.

I have series of polylines representing routes and require the linear distance from the midpoint (along the polyline line) to the start (or end) point of the polyline.

So far I have been able to crack creating the points: I created the points using Qchainage pulgin and now have polylines at the start point, mid point and end point for each polyline. - see below the attributes of these points where "fid" relates back the polyline id and "cng_(meter" is the distance along the polyline created by Qchainage.

I have created a column to easily select only the start points ("START" column), by selecting only "0" values and now need to populate the "MID" column.

I now require selecting only the midpoint values. (A task similar to this thread in Arc Query to select MIN and MAX values by groups)

Obviously this is easily done in excel using the dbf, but I wish to challenge my skills in writing expressions for the calculator or in python.

Is there an easy way to do this?

I would then use the MMQGIS "Hub distance" to create linear polylines from the start point to the mid point. Does anyone have any other thoughts on this process too?


texts: Understanding why hotspots from Getis-Ord Gi* are perfect circles?

texts:
I've been using optimized hotspot analysis and the hotspot maps I am getting back frequently have rings of perfectly circular hotspots.

Is this because of the way Getis-Ord Gi* tests calculates neighborhood sums?

Natural events don't often happen with perfect symmetry and when I look at the underlying data, it doesn't appear to be set out in circles. I feel like it must be due to the statistical processes underlying hot spot analysis.


texts:
I am developing a standalone app using qgis librarys and run the python inside qgis. I import the librarys with no errors but I have one last error: ImportError: no module named site. I define the paths like this:

 'PYTHONHOME': 'C:\OSGeo4W64\apps\Python27',
 'PYTHONPATH': 'C:\OSGeo4W64\apps\Python27\Lib',
 'PATH': 'C:\OSGeo4W64\apps\qgis\bin;C:\OSGeo4W64\bin;C:\OSGeo4W64\apps\Python27'

I see that site.py is in 'C:\OSGeo4W64\apps\Python27\Lib', so I still have the error. The app doesn't read the correct site.py?


texts: Identify which map services are locking feature class

texts:
I have a feature class in one of my enterprise geodatabases (10.3.1) that I would like to modify but it is locked by a number of ArcGIS Map Services.

Is there a way to identify specifically which Map Services are locking my feature class?  In the GeoDatabase Administrator window (from ArcCatalog) I can see the server and user, however that server and user has hundreds of map services running, and I can't just shut them all off to modify just the one feature class.

Is there a tool/function/command/URL I should be looking at that could tell me specifically which services I need to stop to release the lock on my feature class?


texts: How to get the pixels from a geotiff file in gdal python for a given bbox?

texts:
I have a GeoTIFF file in WSG 84 Projection reference (it's a GeoTIFF of SRTM). I can use python & gdal to read the whole lot of data into an array (with ReadAsArray()), but I want to only have a subset of the data, everything inside one bounding box (defined by a top latitude, bottom latitude, right longitude, left longitude). But I'm very new to gdal, and I can't figure out how to do it.

Here is some python code to confirm it is in WSG 84, etc.

>>> from osgeo import gdal
>>> gg = gdal.Open('srtm_34_02/srtm_34_02.tif')
>>> gg.GetGeoTransform()
(-15.000416327675339, 0.0008333333333333334, 0.0, 55.00041712672993, 0.0, -0.0008333333333333334)
>>> gg.GetProjectionRef()
'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]'

It doesn't matter if the resultant array of data includes points that are just outside the bbox, due to rounding. Just so long as it doesn't exclude any points that are in the bbox, and it's as close enough as is practical.


texts: Calculation of a text value, The value of which is separated by a special character

texts:
How can I calculate a text value, the value of which is separated by a special character?

For example: dog | duck, I would like to keep the first value: dog
or
elephant | cat,  I want to keep elephant.

The length of characters are different.


texts: Deleting control points in ArcGIS Georeferencing?

texts:
I'm trying to georeference a raster in ArcMap 10.2 using a rivers vector (polyline) feature class. I started by shifting and rotating the raster. However when deleting a control point to reduce my RMS, both shift and rotation are "being automatically" canceled and I need to shift and rotate the raster back again.

It only occurs in "Auto Adjust" mode.;

Does somebody have a clue about the origin of the problem?


texts: How to find all points with certain elevation difference between them in QGIS

texts:
I have a point layer (xyz) and want to find all points that have elevation difference of 200m or more. I would like to then draw a line between each pair with the largest elevation difference.

The goal is to identify areas that have large elevation difference and (eventually) nearest distance between them in QGIS.


texts: Where can I find more information regarding using the Python interface in QGIS with the Time Manager plugin. I made an animation in QGIS and am trying to automate some processes with Python, however I cannot find much info on the internet.


texts: Atlas generation: get attributes of all features from the current view

texts:
I would like to show the attributes of every feature in the current view of the atlas.

The background: On every page of the atlas is a map with some markers. The markers have got numbers and a text description. I would like to show the text description of the markers (which are in the current view) in a text field on top of the map.


texts: ArcMap - ArcObjects Text in Layout mode not showing

texts:
Refactoring VBA to VB.NET and most parts of the code is working expect for the part where I want to show the layername with today's date in the layout mode (pTextElement.Text)

Any ideas why the text is not showing up?

Public Sub PlotTag()

        Dim x As Long
        Dim y As Long

        Dim pMxDoc As IMxDocument
        Dim pActiveView As IActiveView
        Dim pGraphicsContainer As IGraphicsContainer
        Dim pLayer As ILayer
        Dim pMaps As IMaps
        Dim pMap As IMap

        pMxDoc = My.Document

        If Not pMxDoc.ActiveView Is pMxDoc.PageLayout Then Exit Sub

        pActiveView = pMxDoc.ActiveView
        pLayer = pMxDoc.SelectedLayer

        If pLayer Is Nothing Then
        MsgBox("Please choose a layer!")
            Exit Sub
        End If

        pMaps = pMxDoc.Maps
        pMap = pActiveView.FocusMap
        Dim strLayerName As String
        Dim strMapName As String
        strLayerName = pLayer.Name
        strMapName = pMap.Name

        MsgBox("Data Frame is: " & pMap.Name & " " & pLayer.Name)

        Dim pTextElement As ITextElement
        pTextElement = New TextElement
        Dim dteDate As Date
        dteDate = Today
        pTextElement.Text = "Layout" & strLayerName & ". Plotted " & dteDate & "."

        Dim pElement As IElement
        pElement = pTextElement
        pGraphicsContainer = pMxDoc.PageLayout
        Dim pGraphicsContainerSelect As IGraphicsContainerSelect
        Dim pGraphicsContainerProp As IGraphicsContainerProperty

        pGraphicsContainerSelect = pGraphicsContainer
        pGraphicsContainerSelect.UnselectAllElements()

        pElement.Geometry = pActiveView.ScreenDisplay.DisplayTransformation.ToMapPoint(x, y)
        pGraphicsContainer.AddElement(pTextElement, 0)
        pGraphicsContainerSelect.SelectElement(pTextElement)

        Dim pPageLayoutGraphicsSelect As IGraphicsContainerSelect
        pPageLayoutGraphicsSelect = pMxDoc.PageLayout
        pPageLayoutGraphicsSelect.UnselectAllElements()
        pMxDoc.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, Nothing, Nothing)
    End Sub


texts:
I have a mbtiles file and want to publish it as wms. My config:

services:
  demo:
  wms:
    srs: ['EPSG:3857']
    image_formats: ['image/png']

layers:
   - name: my_layer
     title: WMS layer from tiles
     sources: [mbtiles_cache]

caches:
  mbtiles_cache:
    sources: []
    grids: [GLOBAL_WEBMERCATOR]
    cache:
      type: mbtiles
      filename: D:/dev/mbtiles/export/vtagil.mbtiles

But I don't see anything in demo or QGIS.

What have I missed?


texts: QGIS : strange tiedpoints after makeGraph

texts:
I try to use the Dijkstra shortest length path calculation.

On the following program, when I use the Network.shp file (found here : https://geospatialpython.googlecode.com/svn/network.zip) with  pStart = QgsPoint(-0.8095638694, -0.1578175511) and pStop = QgsPoint(0.8907435677, 0.4430834924), the value of tiedPoints is [(-0.814696,-0.121406), (0.880931,0.394498)], so OK. (CRS = WGS 84)

When I just change the Shapefile  to "routeToulouse.shp" and the pStart and pStop coordinates to pStart = QgsPoint(572777,6277751) and pStop = QgsPoint(571644,6279403)  (CRS = RGF93 / Lambert-93), I got this for tiedPoints :
[(9999,-7.23497e-310), (9999,-7.23497e-310)]

Why do I have this strange values?

I use QGIS Desktop to draw the "routeToulouse.shp" file and the 2 points and you can see that the 2 points are on the extent of the map, not outside.

Here's the source code :

# -*- coding: utf-8 -*-
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *
import os,sys

class MapViewer(QMainWindow):
    def __init__(self):
        QMainWindow.__init__(self)
        canvas = QgsMapCanvas()
        canvas.useImageToRender(False)
        canvas.setCanvasColor(Qt.white)
        layer = QgsVectorLayer("routeToulouse.shp", "layer1", "ogr")
        #layer = QgsVectorLayer("Network.shp", "layer1", "ogr")
        if not layer.isValid():
            raise IOError("Invalid shapefile")
        director = QgsLineVectorLayerDirector(layer, -1, '', '', '', 3)
        properter = QgsDistanceArcProperter()
        director.addProperter(properter)
        crs = layer.crs()
        print crs.description()
        builder = QgsGraphBuilder(crs)
        pStart = QgsPoint(572777,6277751)    #1.42510,43.58740,572777.082,6277751.358
        pStop = QgsPoint(571644,6279403)   #1.41067,43.60206,571644.110,6279403.609
        #pStart = QgsPoint(-0.8095638694, -0.1578175511)
        #pStop = QgsPoint(0.8907435677, 0.4430834924)
        tiedPoints = director.makeGraph(builder, [pStart, pStop])
        print tiedPoints
        QgsMapLayerRegistry.instance().addMapLayer(layer)
        layerExtend = layer.extent()
        canvas.setExtent(layerExtend)
        canvas.setLayerSet([QgsMapCanvasLayer(layer)])
        layout = QVBoxLayout()
        layout.addWidget(canvas)
        contents = QWidget()
        contents.setLayout(layout)
        self.setCentralWidget(contents)
        canvas.show()

if __name__ == '__main__':
    QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'], True)
    QgsApplication.initQgis()
    app = QApplication(sys.argv)
    app.processEvents()
    frame = MapViewer()
    frame.show()
    app.exec_()


texts: How to get long, lat and altitude info from geotiff file using GDAL and C++

texts:
I need to get longitude, latitude and altitude info for each pixel in GeoTiff image. I have installed GDAL in my system. Since I am new to GDAL I don't know how to proceed further.

How can I get the information that I need using GDAL and C++?


texts: How to associate function to custom controls in google maps api?

texts:
In the following code I've added a custom control on a map using Google Maps API V3.

    <title>Simple Map</title>
    <meta name="viewport" content="initial-scale=1.0">
    <meta charset="utf-8">

Now I'd like to associate a function at the custom control tha simply return the click coordinates when a user click on the map.

I can do it (without associate the function to the customcontrol ... ), in this manner

  google.maps.event.addListener(map, 'click', function (e) {
                  alert("Latitude: " + e.latLng.lat() + "\r\nLongitude: " + e.latLng.lng());
              });

but I'd like to activate / deactivate this function clicking on my custom control.

Suggestions / examples?


texts: I've developped a plugin in QGIS which can generate some maps automatically. But, for certains legend's entry the text is too large for the block. It it possible in QGIS (python) to wrap legend's text? Not with a caracter, but automatically?


texts: How to determine area of pixel in meters with GDAL and Python?

texts:
I have a raster loaded in Python using GDAL bindings.  I'd like to be able to determine the area of a pixel in meters although I'm not confident that the projection is in meters.  On this page I see a reference to OGRSpatialReference::GetLinearUnits(), but no corresponding Python bindings.

Any suggestions on how I can determine the are of a pixel in meters?

EDIT: Here's an example of how to get at the linear units:

srs = osr.SpatialReference()
srs.SetProjection(dataset.GetProjection())
print srs.GetLinearUnitsName(), srs.GetLinearUnits()


texts: Monteverdi Optical Calibration: Spectral Response File Error in input format

texts:
I have been trying to upload a spectral response file for optical calibration in the Monteverdi (standalone) program. I have followed the required format, but continue to get the persistent error: ERROR: itk::ERRORERROR: itk::ERROR: SpectralSensitivityReader(0x7f8358890340): Invalid file format: SpectralSensitivityReader(0x7f835881a3b0): Invalid file format

My input file is a simple text file for a Pleiades image with 4 bands (min, max, stepsize and then the sequential spectral response values at that stepwise), please see below. I haven't had success in finding an example input spectral response file and only have the guidelines in the software:

0.428 0.552 0.0025
0.007
0.069
0.093
0.115
0.125
...
0.496 0.621 0.0025
0.008
0.015
0.023
...
0.598 0.746 0.0025
0.061
0.119
0.226
...
0.752 0.948 0.0025
0.066
0.111
0.157
...


texts:
I am using ArcGIS pro for the first time. I have a csv excel sheet with X, Y and ID.

I have tried to make an XY event layer but I get the 000055 error

How do I add xy data in this system?


texts: Walkthrough: Add-In Command Bars root menu

texts: I am wanting to create a new ArcMap Add-In project with just a custom root menu. According to the walkthrough, I can access the Add-In Command Bars (which has the Menu element) by adding a new item to a project. What kind of project do I create in the first place given the selection of Add-In Types?


texts: GeoServer Layer Preview Opens OpenLayers instead of GeoExplorer

texts: I am using OpenGeo 3.0 and trying to get GeoExplorer to launch through the Layer Preview. I imported PostGIS tables into GeoServer and created a layer.  However, when I go to view the layer and select "GeoExplorer-->Go", it opens OpenLayers.  Is there something I need to configure?.


texts:
Having a bit of trouble with the basics of getting the (awesomely discreet) new user prompt/warning to properly display.  Here is the (troubled) code I am using:

self.iface
...
iface.messageBar().pushMessage("Error", "Please select a Topographic Map for your simulation!", level=QgsMessageBar.CRITICAL, duration=2)

...and a bit more complete code for help in resolving this issue:

class GraphicGroundwater:
def __init__(self, iface):
    # Save reference to the QGIS interface
    self.iface = iface
    # initialize plugin directory
    self.plugin_dir = os.path.dirname(__file__)
    # initialize locale
    locale = QSettings().value("locale/userLocale")[0:2]
    localePath = os.path.join(self.plugin_dir, 'i18n', 'graphicgroundwater_{}.qm'.format(locale))
    # affirm window orientation
    if os.path.exists(localePath):
        self.translator = QTranslator()
        self.translator.load(localePath)
        if qVersion() > '4.3.3':
            QCoreApplication.installTranslator(self.translator)
# run add base
def run3(self, iface):
    if FindDirectory._destination is 0:
        iface.messageBar().pushMessage("Error", "Please select a central folder location to save your base layer(s)!", level=QgsMessageBar.CRITICAL, duration=2)
    elif FindDirectory._destination is not 0:
        if GetMaps._TOPOcrs is 0:
            iface.messageBar().pushMessage("Error", "Please select a Topographic Map for your simulation!", level=QgsMessageBar.CRITICAL, duration=2)
        elif GetMaps._TOPOcrs is not 0:
            crs = str(GetMaps._TOPOcrs)
            Base =  QgsVectorLayer('Polygon?crs=' + GetMaps._TOPOcrs + '&field=LAYER:string(8)&field=BOTTOM:string(8)&field=TOP:string(8)&field=K:string(10)&field=R:string(10)&index=yes', "Base", "memory")
            QgsMapLayerRegistry.instance().addMapLayers([Base])
            pr = Base.dataProvider()
            BaseFeat = FindDirectory._destination+"/Base.shp"
            writer = QgsVectorFileWriter(BaseFeat, "CP1250", pr.fields(), QGis.WKBPolygon, pr.crs(), "ESRI Shapefile" )
            del writer

And here is the resultant conflict:

  Traceback (most recent call last):
    File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\GraphicGroundwater\graphicgroundwater.py", line 132, in run3
      iface.messageBar().pushMessage("Error", "Please select a central folder location for your simulation!", level=QgsMessageBar.CRITICAL, duration=2)
  AttributeError: 'bool' object has no attribute 'messageBar'

Any clues for a qgis_python_novi?


texts: ArcGIS Server License Update - Additional steps for Enterprise GeoDatabase?

texts:
I've already posted more or less the same question at esri's GeoNet Forum, but until now I have got any response. As my license will expire shortly, the is my second try to get some information before license-expiration.
Last year, I have struggled with the some problems with my enterprise geodatabase after my license expired.

My environment is the following:

Windows Server 2012
ArcGIS Server 10.3.1
GeoEvent Processor 10.3.1
Enterprise Geodatabase (postgreSQL) on the same machine as ArcGIS
Server/GeoEvent Extension

My current licenses will expire after the 29th of September (so in 1,5 hours) but I have already a new license-file for the next year.
For the license-update process, I would perform the following steps:
Using the license-file and the ArcGIS Administrator I'll authorize my ArcGIS Server instance including the GeoEvent Processor.

To avoid problems with my enterprise geodatabase: Are there any additional steps, which I have to perform (like inject the license key directly in my database)?!
And what's about my GeoDatabase-Connection from my personal computer to the enterprise geodatabase? Do I have to establish a new connection using a "new generated" keycodes file?


texts: Atmospheric correction of larger collection of Sentinel 2 L1C images (in Google Earth Engine) for NDVI time series

texts:
System: Windows 10

Software: Google Earth Engine browser, (QGIS, SNAP)

Data: Sentinel 2 L1C & L2A, NDVI calculations

I'm currently working on generating a NDVI timeseries in Google Earth Engine for 135 points in an area which is about 750m². I want to do this with Sentinel 2 imagery, starting in 2015.
Since L2A imagery only starts at dec 2018, I want to atmospherically correct all L1C images of my ROI before that timeperiod, but so far I have only found methods to correct the images one-by-one in programs such as FMask and SNAP with Sen2Cor.
I appreciate the fact that these programs produce images of similar quality as the L2A images, but is there a program or script that can correct more images at once (also with a method like Sen2Cor), preferably in Google Earth Engine?


texts: Polygon shapefile doesn't render certain polygons at full extent

texts: This is a snapshot of a polygon in arc catalog, in preview mode. Even when brought into ArcMap, it won't display the polygons that are in the white space until you zoom in to an extent of 50,000 or less. I looked and there is nothing in the data frame properties to cause this and the scale range for this polygon is set to show layer at all scales. There aren't any definition queries either that would only display geologic units of certain types. Any clue on how to bring these polygons back from the void?


texts: SQL View- WMS Viewparams - parameter for table name

texts:
I've made a SQL view using parameteres.
So, I did something like this:
SELECT * FROM '%tablename_parameter%'
The view is requesting data from a table by a parameter inserted from web page.
The geometry layer in table is Polygon or MultiPolygon.
The view is working, it shows data but it shows points not polygon.
I've saved the sql view to show me Polygon or MultiPolygon but still shows me points.
I don't know why or how to resolve this. Can someone help me?


texts: pgrouting: function driving_distance does not exist

texts:
I tried to run driving distance function but got error.

   function driving_distance(unknown, integer, integer, boolean, boolean) does not exist
    LINE 1: SELECT * FROM driving_distance('
                          ^
    HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

SELECT * FROM driving_distance('
          SELECT gid AS id,
              start_id::int4 AS source,
              end_id::int4 AS target,
              shape_leng::float8 AS cost
          FROM roads',
          15137,
          100000,
          false,
          false)


texts: Setting band names when writing multiple layer rasters using GDAL with Python?

texts:
I want to write a raster file (using GDAL+Python) with multiply bands. However, I do not know how to named bands. Their names are always "Band 1", "Band 2", etc. I search and found a function that should do the trick but failed (SetRasterCategoryNames). I hope that someone could help me how to name when writing multiple-layer rasters in python.

Following code is just an example that writing multiple bands to a raster.

        from osgeo import gdal

        #Read projection information from original raster
        path = '/home/12101008/MCD43A1.A2007001.h26v06.006.2016111212159.hdf'
        GRID_NAME = 'MOD_Grid_BRDF'
        DATAFIELD_NAME = 'BRDF_Albedo_Band_Mandatory_Quality_Band1'
        gname = 'HDF4_EOS:EOS_GRID:"{0}":{1}:{2}'.format(path, GRID_NAME, DATAFIELD_NAME)

        gdset = gdal.Open(gname)
        cols = gdset.RasterXSize
        rows = gdset.RasterYSize
        geoTransform = gdset.GetGeoTransform()
        proj = gdset.GetProjection()

        #Read example layer to test
        result1 = gdset.ReadAsArray()
        #print geoTransform, proj, cols, rows

        #Setting driver
        driver = gdal.GetDriverByName('NetCDF')
        driver.Register()

        #Create new file and set parameters from original raster
        path1 = '/home/12101008/Processed/test.nc'
        outDataset = driver.Create(path1, cols, rows, 6, gdal.GDT_Int16, ['COMPRESS=PACKBITS'])
        outDataset.SetGeoTransform(geoTransform)
        outDataset.SetProjection(proj)

        #Write 6 raster from "result" variable
        for i in range(0,6):
            outBand = outDataset.GetRasterBand(i+1)
            #Attempted to set band name like this but failed
            outBand.SetRasterCategoryNames(['NDVI'])
            outBand.SetNoDataValue(-10000)
            outBand.WriteArray(result1, 0, 0)

        gdset = None
        outBand = None
        outDataset = None


texts: Display in legend only features present on a map extent

texts:
I have lots of layers with lots of rule-based styling in my QGIS project, although I would like to display in legend for my print composer only those features that are actually present on a current map extent.

In ArcGIS I would use an option "Only show classes that are visible in the current map extent", so I wonder if there's something similar in QGIS.

Of course there is a manual solution but it is not an option in many cases.


texts:
I would like to compare vector vs raster data. So i try to fix a problem using both ways.

I made it in vector first :

coastline<-readShapeSpatial("shp")
wets<-readShapeSpatial("shp")
wetb<-gBuffer(wets,width=500,byid=T) #doing a buffer of 500m of the wetlands
wetb2<-gBuffer(wets,width=1000,byid=T) #doing abuffer of 1000
wetb1$vat=0 # assigning a new column of 0 in the buffer i created
wetb2$vat=5 # assigning a new column of 5 in the second buffer
coastline$vat=10 #assigning a new column of 10 in the rest of the land

Now how would i write the same thing in raster?

My thoughts are that i should turn the coastline and the wetlands into raster first off.How would i do that and what after that?


texts:
i am running arcgis server 10, installed visual web developer express 2010 edition, now wanted to run samples from arcgis.com for .net web api, for example the following one

http://help.arcgis.com/en/sdk/10.0/serveradf_net/conceptualhelp/index.html#/ADF_Tutorials/000200000025000000/

I downloaded the vb.net zipped folder shown for all files on above link. Unzipped it. Opened the "ADFTutorials_VBNet2010.sln" file. The visual web developer 2010 express says that

 One or more projects in the solution were not loaded correctly.
 Please see the Output Window for details.

The output shows the following

 http://localhost/ADFTutorials_VBNet/ScriptableControls/MapCoordinateDisplayWebSite :
 error  : Unable to open the Web site 'http://localhost/ADFTutorials_VBNet
 /ScriptableControls/MapCoordinateDisplayWebSite'.  The Web site 'http://localhost
 /ADFTutorials_VBNet/ScriptableControls/MapCoordinateDisplayWebSite' does not exist.

 http://localhost/ADFTutorials_VBNet/AccessArcGISServer : error  : Unable to open the
 Web site 'http://localhost/ADFTutorials_VBNet/AccessArcGISServer'.  The Web site
'http://localhost/ADFTutorials_VBNet/AccessArcGISServer' does not exist.

 http://localhost/ADFTutorials_VBNet/AccessArcIMS : error  : Unable to open the Web
 site 'http://localhost/ADFTutorials_VBNet/AccessArcIMS'.  The Web site
 'http://localhost/ADFTutorials_VBNet/AccessArcIMS' does not exist.

 http://localhost/ADFTutorials_VBNet/ADFWebControls : error  : Unable to open the Web
 site 'http://localhost/ADFTutorials_VBNet/ADFWebControls'.  The Web site
 'http://localhost/ADFTutorials_VBNet/ADFWebControls' does not exist.

  http://localhost/ADFTutorials_VBNet/UsingCommonAPI : error  : Unable to open the
  Web site 'http://localhost/ADFTutorials_VBNet/UsingCommonAPI'.  The Web site
  'http://localhost/ADFTutorials_VBNet/UsingCommonAPI' does not exist.

Please can anyone explain how to run the samples provided for .net? thanks a lot in advance!!!


texts: Creating a spatial reference in arcpy from projection name

texts:
I am using ArcGIS 10.2.2, Python 2.7. I want to create a spatial reference from the name of a projection, e.g. NAD_1983_BC_Environment_Albers, using the arcpy tool arcpy.CreateSpatialReference_management().
Consider this:

projName = "NAD_1983_BC_Environment_Albers"
#Replace underscores with spaces, works fine.
sr = arcpy.CreateSpatialReference_management(projName.replace('_',' '))

No issues here. Also, I tried the example in the arcpy documentation here:

projName = "NAD 1983 StatePlane California V FIPS 0405 (US Feet)"
sr = arcpy.CreateSpatialReference_management(projName)

This also works fine. So, why does the following generate an error?

projName = "NAD_1983_UTM_Zone_9N"
#Tried both with and without converting underscores to spaces
sr = arcpy.CreateSpatialReference_management(projName.replace('_',' '))

ERROR 000622: Failed to execute (Create Spatial Reference). Parameters are not valid.
ERROR 000628: Cannot set input into parameter spatial_reference.

Any ideas??


texts:
I am quite new on spatial R. I would like to apply the land mask to monthly SST data from https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.highres.html.
I tried to use "mask" from the raster package but the output raster (NoLand) is missing the Z (Date).

nc = raster("N:/Papier DLM/sst_extraction.nc")
mask = raster("N:/Papier DLM/lsmask.nc")
mask2 <- crop(mask,extent(nc))
NoLand<-mask(nc, mask2)
writeRaster(NoLand, "N:/Papier DLM/NoLand.nc")

What should I use to apply the mask over time (for each month)?


texts: Where is the work in the GIS/RS consulting industry?

texts:
I'm trying to get an understanding of what type of GIS/RS work organizations are most often bidding out to consultants, and why?  For example, is most of the freelance work in web mapping these days because in-house programmers have not caught up with this rapidly expanding field?  Maybe all the work involves niche specialty jobs, such as GIS/RS research support for temporary environmental projects because non-profit organizations do not have the expertise.  Perhaps, oil and gas jobs are dominating the freelance market because of the fracking boom.

What type of GIS/RS work are private/corporate/governmental organizations most often bidding out to consultants these days, and why?


texts: Fast Union with Known Geographic Locations

texts:
I am trying to do a union across thousands of different geometries. The geometries can be very complex so a naive union can take many hours. I then want to dump these geometries so each polygon becomes its own shape rather than one very large multipolygon.

However, I have extra information about the geometries. The total area covered by the geometries is split up into squares. Each geometry is contained within a square and thus can only touch the geometries of the 8 adjacent squares. Is there some way to use this extra information about which polygons could possible intersect in order to reduce the time of the union?

I can do this either in PostGIS or in Python.


texts:
When I project (lon,lat) coordinates I get a different result depending on the version of OGR/OSR.

osgeo.version_info

sys.version_info(major=3, minor=6, micro=9, releaselevel='final', serial=0)

import ogr, osr
lon,lat = 6.16800,44.88358
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(lon, lat)
inSpatialRef = osr.SpatialReference()
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(4326)
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(32631)
coordTransform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
point.Transform(coordTransform)
point.GetX(),point.GetY()

(750192.6531073075, 4974901.243746961)

When I run the same code with another distribution of python, osgeo, I need to swap lon and lat in AddPoint to get the correct result.

point.AddPoint(lon, lat)

osgeo.version_info

sys.version_info(major=3, minor=6, micro=7, releaselevel='final', serial=0)

Is it documented somewhere? I would like to handle this in my code but I'm not sure from which version it has changed..


texts: How to undo/redo the splitting/union of a feature

texts:
I have looked everywhere for examples on this however all I can find are the normal operations, add delete edit and move. Since split and rotate involve creating two features instead of one or combining two features into one, I am not sure how the process would work out.

For example the process for moving a vertex would be as follows:

//evt1 = dojo.connect(dojoEditToolbar, 'onVertexFirstMove', function (evt) {
                    //    oldGraphic = evt;
                    //});

                    //evt2 = dojo.connect(dojoEditToolbar, 'onVertexMoveStop', function (evt) {
                    //    if (oldGraphic != null && oldGraphic != undefined) {
                    //        var operation = new esri.dijit.editing.Update({
                    //            featureLayer: evt._graphicsLayer,
                    //            preUpdatedGraphics: [oldGraphic],
                    //            postUpdatedGraphics: [evt]
                    //        });

                    //        myUndoManager.add(operation);
                    //    }
                    //});


texts: WMS GetCapabilities format does not return a wms content

texts:
My geoserver WMS GetCapabilities request returns an XML file download popup. But I want to see xml result on my browser. And I want to request via ajax request. But could not.

WMS Get capabilities like this:

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE WMT_MS_Capabilities SYSTEM "http://schemas.opengis.net/wms/1.1.1/capabilities_1_1_1.dtd" >

<WMT_MS_Capabilities version="1.1.1">
  <Service>
    <Name>OGC:WMS</Name>
    <Title>Web Map Service</Title>
    <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://.."/>
  </Service>
  <Capability>
    <Request>
      <GetCapabilities>
        <Format>application/vnd.ogc.wms_xml</Format>
        <DCPType>
          <HTTP>
            <Get>
              <OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="http://..."/>
            </Get>
          </HTTP>
        </DCPType>
      </GetCapabilities>


texts: How create a field of sequential number with autoincrement() that start from zero in arcpy?

texts:
I would like to create sequential numbers in a field of a  point shapefile  that start from 0 with an interval of 0.050.
The main code works whenever the start number is not zero.

The original code to put in the codeblock from (https://support.esri.com/en/technical-article/000011137) looks like this:

rec=0
def autoIncrement():
   global rec
   pStart = 1
   pInterval = 1
   if (rec == 0):
      rec = pStart
   else:
      rec += pInterval
   return rec

How I can make it in arcpy with a “codeblock” and “Expression”?


texts: How to serve both WFS and WFS-T from one server?

texts: In GeoServer, how can I serve two types of WFS from one server? Setting the type of the service is done once in global settings. Is it possible to have WFS from workspace A and WFS-T from workspace B?