texts: GeoServer with viewparams returns ORA-00907: missing right parenthesis

texts:
I am trying to call a JSON request from geoserver which has an Oracle NG Layer.  I have one parameter in the call:

viewparams=customer:One Customer

The parameter value has to have a space in it.  I cannot get it to work, and am getting the following error:

<ServiceExceptionReport version="1.2.0" xsi:schemaLocation="http://www.opengis.net/ogc http://schemas.opengis.net/wfs/1.0.0/OGC-exception.xsd"><ServiceException>
      java.lang.RuntimeException: java.io.IOException
java.io.IOException
null
ORA-00907: missing right parenthesis

</ServiceException></ServiceExceptionReport>

This is the url:

http://pc121212:8085/geoserver/sf/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=sf:view_ae&viewparams=customer:One%20Customer&outputFormat=JSON

When I go into "Edit SQL View" for the layer in GeoServer, this is the query:

select * from view_ae
where customer = %customer%

The parameter in GeoServer is added with the following values:

Name = customer
Default value = 'One Customer'
Valid Regular Expression = ^[\w\d\s]+$

Does anyone know what is wrong?


texts: Trying to define a range of values as NoData in QGIS

texts: I have a raster from GEBCO with pixel values ranging from -5888 to 5485 m. I only want to analyze the Bathymetry so I am looking for a way to assign noData values to all values >0. I also have a polygon of the coastline so if there is some way to do a clip outside of a polygon that would work also. I am working in QGIS. Thanks


texts:
I am using ArcGIS 10.6 ArcMap License.

I am looking for a way to calculate the distance between polygons and the centerline of the street network. The aim to use this vale as the buffer value to identify the polygons that face the street (line shape).
For instance I tried the near feature, then taking near distance + tolerance as the buffer value. But still some polygons are not covered by the buffer and others totally covered depending on street width. When I increase the buffer value, that cause problem for the next step modeling. How Can I add the value to the buffer_value dynamically as the function of street width? Or to increase the buffer_value up to the limit that intersect with the near features polygon? I tried iterate but didn’t gave desired result.


texts: What is best way to convert DWG to KML using ArcGIS or alternatives?

texts:
What's the best way to convert Dwg file to Kml ?

i have tried to convert Dwg file to shape file then convert it to KML, but i hadn't very good result , i had some shift from where should be seen in google earth, i heard about some software performs better than others.

i would like to know also the reasons of the shift, the DWG file used contains projected coordinates.


texts: GDAL tools for QGIS on OSX - what is the exact directory for gdal binaries?

texts:
When I try to use the gdal raster tools, eg clipper, I get the message "The process failed to start. Either the invoked program is missing, or you may have insufficient permissions to invoke the program."  I have read the forums discussing this issue, and they seem to have solved it by using multiple paths for python and pymods, separated by "/usr/bin:/Library/Frameworks/GDAL.framework/Versions/1.7/Python/"

I am not sure of what this means, but I do know that setting the path to gdal binaries in the gdal tools settings to the path above does not work.

Is there a dumbed-down version of how to correctly set the paths in qgis?  Any help would be greatly appreciated.  The OS I am using is OSX 10.5 (Leopard), and I am using the 'GDAL complete' package.  I tried versions 1.7 and 1.8, not sure if that makes a difference.


texts: Setting properties of geoJSON in OpenLayers

texts:
I have a map with layers that have similar attributes. I want to add a function that when you select a feature in a layer, it searches for a feature with similar properties in another feature (in this case NAME) and if a similar feature property is found, then the styling of the matching feature is altered.

I have a JSfiddle that tries to achieve this here: http://jsfiddle.net/TimLucas/nmd2vsww/6/

The function is called by using the second button. As you can see the get() of features works from the function called by the first button.

Here I mainly have problems with how you can set properties (instead of getting them) in the VectorSource GeoJSON file. I use .attr() but I am not sure if that actually works.

Does anyone know how to do this?


texts: Mirror of Michigan State University's Python geoprocessing course materials

texts:
Been learning geoprocessing with Python using MSU's GEOG 825 ( python geoprocessing ) course. Unfortunately the link is dead now. I know this is a long shot but does anybody have a copy of the materials from that course? I would really like to finish it.

Already checked out the other online Python geoprocessing courses. Found one in Penn State and on Utah State. Those focus on using tools though, whether they be ogr/gdal or ArcGIS' arcpy/arcgisscripting. Only Michigan State U's course discusses the algorithms and data structures as far as I know.


texts: Why will my field calculator python function always return Null if any inputs are Null?

texts:
I've written a simple python function that outputs a list based on yes/no fields. It appends a word to a list if it has a "yes" in another column then writes the list, formatted according to the number of values in it.

It looks like this:

def proctype(prestm, presdm, postm, posdm, beam, wem):

    mylist = []
    if prestm == 'YES':
        mylist.append("PreSTM")
    if presdm == 'YES':
        mylist.append("PreSDM")
    if postm == 'YES':
        mylist.append("PoSTM")
    if posdm == 'YES':
        mylist.append("PoSDM")
    if beam == 'YES':
        mylist.append("Beam")
    if wem == 'YES':
        mylist.append("WEM")

    if len(mylist) == 1:
        return mylist[0]
    if len(mylist) == 2:
        return mylist [0] + " and " + mylist[1]
    if len(mylist) == 3:
        return mylist[0] + ", " + mylist[1] + " and " + mylist[2]
    if len(mylist) == 4:
        return mylist[0] + ", " + mylist[1] + ", " + mylist[2] + " and " + mylist[3]
    if len(mylist) == 5:
        return mylist[0] + ", " + mylist[1] + ", " + mylist[2] + ", " + mylist[3] + " and " + mylist[4]
    if len(mylist) == 6:
        return mylist[0] + ", " + mylist[1] + ", " + mylist[2] + ", " + mylist[3] + ", " + mylist[4] + " and " + mylist[5]

This works exactly as I want it to if the fields its looking at are 'yes' and 'no', but if there are any Null values in the six input fields the output will always be null.

Does anyone know why this is? Can the field calculator not handle Null inputs, or does the list.append have a problem with Null values?

I've tried setting the values to a number if theyre null like this:

if prestm is None:
    prestm = 1

But that doesn't help.

EDIT: This problem only occurs when the inputs are from other fields in the table. inputing "" works but

  proctype( !PRESTM!, !PRESDM!, !POSTM! , !POSDM! , !BEAM_MIGRATION! ,
  !WEM_RTM! )

doesnt.

EDIT 2: The problem is not with my code it is with passing in other attributes as parameters where any of them are Null

def testing(prestm, presdm, postm, posdm, beam, wem):
    return "test ok"

Returns "test ok" when all of the fields are populated but null

if any of the input fields are Null


texts: Preventing Google Earth Engine exporting image with black border?

texts:
When I import my image into ArcGIS Desktop I get a black border around the region that didn't exist before in Google Earth Engine (GEE).

How can I prevent GEE from doing this when I export an image?


texts: Google earth engine - set null band values as zero

texts:
I have a function which assigns a score depend on the value of a band, which I then map to an image collection.
However, in some cases, I have null values, where no value is given. In these cases I'd like to convert the null values to 0. I tried to set null values as zero in the code, but I'm getting an image error. My new code, with the extra line is below.
function currentHS_function(image){
  var currentHS = image.expression(
 "(b('current') < 0.13) ? 1.0" +
 ": (b('current') > 0.15) ? 0.0" +
 ": (b('current') = null) ? 0.0" + // new line - set null values as 0
 ": -50.0*current + 15.0/2.0 ",
 {"current" : image});
   return currentHS;
 }

// apply heath score function to image collection
var currentHS_collection = currentSpeed.map(currentHS_function);


texts: problem of using TIMESAT to processing MODIS NDVI

texts:
I want to use TIMESAT software to analyse time series MODIS NDVI data. Firstly, I use MRT tool to mosaic,subset and project MODIS NDVI data and output data as GEOTIFF format, and save as ENVI standard using ENVI software. secondly, using TSM_imageview to display the image file, and the result is ok. then i use TSM_GUI to open list of image files, the problem is arising--"Error while evaluating uimenu Callback" .

so, if anyone knows the answer, please help me, thanks a lot!


texts: Cut line label on tile border in geoserver tile cache

texts:
I have a problem with line label clipped on tile boundary. For point layer labels, it was solved by using "buffer" parameter in tile caching option and using "partials" in SLD. But in my opinion, labelling in line is dynamic, not fixed label. So I can't solve this problem using buffer and partials. How can I solve this problem? I'm beginner in GeoServer.
Here are my tile cache option.

meta tile size : 8*8

gutter size : 100

image format : png8

tile image : 512*512

labeling vendor option
<sld:VendorOption name="allowOverruns">false</sld:VendorOption>
<sld:VendorOption name="followLine">true</sld:VendorOption>
<sld:VendorOption name="group">true</sld:VendorOption>
<sld:VendorOption name="maxAngleDelta">25.0</sld:VendorOption>
<sld:VendorOption name="charSpacing">5</sld:VendorOption>
<sld:VendorOption name="partials">true</sld:VendorOption>

I think there are two alternatives(avoid label or label). But I don't know how to that.

How can I avoid line label overlap with tile border? (I have tile border layer.)
How can I solve clipped line label?


texts: Print and write the WKB to another shapefile using OGR_G_ExportToWKB in C?

texts:
I want to convert the geometry to well known binary (WKB) format of my shapefile. I want to print and want to write it into another shapefile.

My code is:

#include <stdio.h>
#include <stdlib.h>
#include <ogr_api.h>
#include <cpl_conv.h>
#include <cpl_port.h>
#include <gdal.h>
#include <string.h>
#include <gdal_utils.h>
#include <cpl_progress.h>
#include <cpl_minixml.h>
#include <ogr_core.h>
int main()
{
    char path1[200],path2[200],filename1[]="/1.shp",filename2[]="/g_copy.shp";
    unsigned char *buf;
    char *sz;
    GDALDatasetH inShapeDS;
    OGRLayerH inShapeLayer;
    OGRFeatureH inShapeFeature;
    OGRFeatureDefnH inFeatureDefn;
    OGRFieldDefnH inFieldDefn;

    int inFieldCount;
    OGRGeometryH inShapeGeometry;

    GDALAllRegister();

    getcwd(path1, sizeof(path1));
    strcpy(path2,path1);
    strcat(path1,filename1);
    strcat(path2,filename2);

    inShapeDS = GDALOpenEx(path1, GDAL_OF_VECTOR, NULL, NULL, NULL );
    if( inShapeDS == NULL )
    {
        printf( "Open failed.\n" );
        exit( 1 );
    }

    printf("number of layers:%d\n",GDALDatasetGetLayerCount(inShapeDS));

    inShapeLayer = GDALDatasetGetLayerByName(inShapeDS, "1");

    inFeatureDefn = OGR_L_GetLayerDefn(inShapeLayer);
    printf("number of features in this layer:%d\n",OGR_L_GetFeatureCount(inShapeLayer,1));

    inShapeFeature = OGR_L_GetNextFeature(inShapeLayer);

    inShapeGeometry = OGR_F_GetGeometryRef(inShapeFeature);

    printf("size of binary geom:%d\n",OGR_G_WkbSize(inShapeGeometry));

    buf = (unsigned char *)CPLMalloc(OGR_G_WkbSize(inShapeGeometry));

    OGR_G_ExportToWkb(inShapeGeometry,wkbXDR,buf);

    //sz=buf;
    sz = (char *)buf;

    printf("binary geom:%s",sz);

    const char *pszDriverName = "ESRI Shapefile";
    GDALDriverH outDriver;
    GDALDatasetH outDS;
    OGRLayerH outLayer;
    OGRFieldDefnH outFieldDefn;

    outDriver = GDALGetDriverByName( pszDriverName );
    if( outDriver == NULL )
    {
        printf( "%s driver not available.\n", pszDriverName );
        exit( 1 );
    }
    outDS = GDALCreate( outDriver, path2, 0, 0, 0, GDT_Unknown, NULL );
    if( outDS == NULL )
    {
        printf( "Creation of output file failed.\n" );
        exit( 1 );
    }
    outLayer = GDALDatasetCreateLayer( outDS, "g_copy", NULL, wkbUnknown, NULL );
    if( outLayer == NULL )
    {
        printf( "Layer creation failed.\n" );
        exit( 1 );
    }

    outFieldDefn = OGR_Fld_Create( "GeomB", OFTString );
    OGR_Fld_SetWidth( outFieldDefn, (OGR_G_WkbSize(inShapeGeometry)));
    if( OGR_L_CreateField( outLayer, outFieldDefn, TRUE ) != OGRERR_NONE )
    {
        printf( "Creating Name field failed.\n" );
        exit( 1 );
    }
    OGR_Fld_Destroy(outFieldDefn);

    OGRFeatureH outFeature;
    outFeature = OGR_F_Create( OGR_L_GetLayerDefn( outLayer ) );
    OGR_F_SetFieldString( outFeature, OGR_F_GetFieldIndex(outFeature, "GeomB"), sz );

    if( OGR_L_CreateFeature( outLayer, outFeature ) != OGRERR_NONE )
    {
        printf( "Failed to create feature in shapefile.\n" );
        exit( 1 );
    }
    OGR_F_Destroy( outFeature );

    CPLFree(buf);
    OGR_F_Destroy(inShapeFeature);

    GDALClose(inShapeDS);

    return 0;
}

Now the problem is in print function, it is giving some symbol. And in the newly created shapefile, it is not printing anything. Another thing that in the new shapefile the width of the GeomB field is truncated to 254, but the width of binary version of the geom is almost 2000 bits.

It is getting very much difficult in C. How to fix it?


texts: Why is ST_X and ST_Y of recently inserted GEOGRAPHY(Point) for certain coordinates different than value inserted?

texts:
I have a table of cities that includes the name and location of that city.

CREATE TABLE cities (
    name        varchar (255),
    location    GEOGRAPHY(Point)
)

For certain cities, I'll set the location:

UPDATE cities SET location='POINT(32.4543 -99.7384)' WHERE slug='abilene';

But when I retrieve the longitude and latitude for that location, it is different than the value I inserted:

select ST_X(location::geometry),ST_Y(location::geometry) from cities where slug='abilene'

The ST_X is correct at 32.4543 but the ST_Y is incorrect at -80.2616

Why is the latitude different than the value I inserted?

I am using PostgreSQL 9.6.9 and Postgis 2.5


texts: ArcMap 10: How do I escape certain characters when labeling?

texts:
I am labeling certain Oil & Gas areas by Operator in ArcMap 10, using this expression:

"" & [OPERATOR] & "" & vbNewline & "" &  [UNIT_NAME] & ""

The only problem is that when the operator includes an ampersand (&), the label messes up because it's reading it as part of the expression, and ends up looking like it does in the screenshot below.

I know that in PHP (from my experience programming unrelated to GIS), there is a way to escape certain characters, but I haven't been able to find it in VBScript.

Screenshot: http://i.stack.imgur.com/gdKxB.png


texts: Reading .osm-file: OGR and XML-Parser return different number of ways and nodes in ways

texts:
I try to read in an .osm-file with python and GDAL/OGR, and then I want to save way IDs and associated node coordinates in a dictionary. For some strange reason, OGR seems to skip either entire ways or some nodes of a way.

I get different results when using an ordinary XML-Parser to read in the file - no ways or nodes of ways are skipped.

Code used with OGR:

count = 0
for wayFeature in layer2:
    processWay(wayFeature) # do smthg with the feature
    count=count+1

Code used with ElementTree to parse xml:

for way in root.findall('way'):
    id = way.get('id')
    processWay(way) #do smthg with the way
    count = count+1

With ogr,I get 253 ways, with the XMLParser, I get 794 - a huge difference. 794 is what I expect,as I did a find&count on my .osm file.

How can that be?

I already set the osmconf.ini file to not ignore any nodes or ways,no matter if they have relevant tags or not.

Does it maybe have to do with the way's IDs in the file? I already played around with a dummy file, OGR sometimes seems skip ways where ID's are not ascending...but not always.

I could use the XML parser, but OGR has some nice features like filters and intersection tests, which I would like to use.


texts: How to change Python version for GRASS in Windows

texts:
I'm using GRASS GIS 7.0.3, the default python version is python2.7. I need python3.5 for my own module. After I edit the 'etc/env.bat' file to change GRASS_PYTHON and PYTHONHOME to the location of my new python .exe file and the new python directory. Then I run GRASS GIS, but I got the following error:

  File "C:\Program Files\GRASS GIS 7.0.3\etc\grass70.py", line 251
  os.makedir(tmpdir, 0700) SyntaxError: invalid token

Am I do it in the right way?


texts: QGIS Python script with raster calculations

texts:
I have a problem with this code, that works, because it gives me the output rasters that I want, but the problem is that these outputs show wrong values because they are all 1.79769e+308.
I think the formula is right, because when I print it is what I would expect, and also the output path is good because all the output rasters are saved in the right folder.
Maybe I make some mistakes in importing the right input rasters at every loop?
Because if I make calculations (by commenting the others) only with ras@1 or only with ras@2 or only with ras@3 the code works perfectly, so I guess the problem could be when I import the three rasters together at every loop.
Here is the code:
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
Rg=[6.4, 8.6, 11.4, 14.3, 16.4, 17.3, 16.7, 15.2, 12.5, 9.6, 7.0, 5.7]
d=[31,28,31,30,31,30,31,31,30,31,30,31]
suffix_output=['Harg_jan', 'Harg_feb', 'Harg_mar', 'Harg_apr', 'Harg_may', 'Harg_jun', 'Harg_jul', 'Harg_aug', 'Harg_sep', 'Harg_oct', 'Harg_nov' ,'Harg_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']
suffix_input2=['Tmaxjan', 'Tmaxfeb', 'Tmaxmar', 'Tmaxapr', 'Tmaxmay', 'Tmaxjun', 'Tmaxjul', 'Tmaxaug', 'Tmaxsep', 'Tmaxoct', 'Tmaxnov', 'Tmaxdec']
suffix_input3=['Tminjan', 'Tminfeb', 'Tminmar', 'Tminapr', 'Tminmay', 'Tminjun', 'Tminjul', 'Tminaug', 'Tminsep', 'Tminoct', 'Tminnov', 'Tmindec']

outputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/PET_Hargreaves/'
inputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_average/'
inputpath2 = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_max/'
inputpath3 = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_min/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
for i in list:
  inputrasterfile = QgsRasterLayer(inputpath + suffix_input[i] + ".tif")

  entries = []
  ras = QgsRasterCalculatorEntry()
  ras.ref = 'ras@1'
  ras.raster = inputrasterfile
  ras.bandNumber = 1
  entries.append( ras )

  inputrasterfile2 = QgsRasterLayer(inputpath2 + suffix_input2[i] + ".tif")

  entries = []
  ras2 = QgsRasterCalculatorEntry()
  ras2.ref = 'ras@2'
  ras2.raster = inputrasterfile2
  ras2.bandNumber = 1
  entries.append( ras2 )

  inputrasterfile3 = QgsRasterLayer(inputpath3 + suffix_input3[i] + ".tif")

  entries = []
  ras3 = QgsRasterCalculatorEntry()
  ras3.ref = 'ras@3'
  ras3.raster = inputrasterfile3
  ras3.bandNumber = 1
  entries.append( ras3 )

  new_path = outputpath + suffix_output[i]+ '.tif'

  hargreaves = str(d[i]) + '* (0.0023 *' + str(Rg[i]) +'* (17.8 +  ras@1) * (ras@2 - ras@3)^0.5)'

  print(hargreaves)

  Pet = QgsRasterCalculator(hargreaves, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )

  Pet.processCalculation()
  iface.addRasterLayer(new_path)


texts: How to manage GEE assets with service account (Permission issue)

texts:
I need to export some images to my GEE Assets using the Python API and I am getting this error: "Insufficient permissions to create asset 'users/ramizsami95/kyrgyz/ndvi_anomaly".
This is my code:

task = ee.batch.Export.image.toAsset(
        image=anomaly,
        description='NDVI-anomaly',
        assetId='users/ramizsami95/kyrgyz/ndvi_anomaly')
task.start()

In the Python code, I have authenticated with a service account that is created from my Gmail account and it seems like it can't access my GEE assets. Now I have two questions:
1) Is there a way that I can give owner access to this kyrgyz folder to my service account? (I have given writer permission already)
2) Do we have separate assets for the service account that I am missing?


texts: How to Find Local Lowest Point(s) on a Landform?

texts: I am trying to write a code to find the local lowest points using Python. This approach will be applied to raster landform data such as digital elevation models. What approach can I take to find the lowest point(s) in a digital elevation model?


texts:
Can spatial views be edited via ArcMap and, if it is possible, then how do I do it?

I do have my Enterprise Geodatabase on Oracle.


texts: Curl request to get times from layers stored as an image mosaic in GeoServer

texts:
I have an image mosaic in GeoServer with some .tif which its name are the date of that raster:

20190929.tif
20191004.tif
...

So I can load every raster by passing the date in OpenLayers:

date="2019-09-29";

var myTileLayer = new TileLayer({
  source: new TileWMS({
    url: 'http://localhost:8080/geoserver/myworkspace/wms',
    params: {'LAYERS': 'myworkspace:mylayer','TIME':date,'TILED': true},
    serverType: 'geoserver',
  }),
  extent: extent,
});

The problem is that I want that the user could choose the date for example in a calendar (in the example I put the date manually) but for that I need to know which dates are avalaible in the image mosaic so I was trying to search if there was avalaible any get request in GeoServer to retrieve this information like:

curl -u admin:geoserver -XGET localhost:8080/geoserver/rest/layers/myworkspace:mylayer/time.json
//Of course this doesn't work

Expected return:

{
"time":[20190929,20191004]
}

The other approach I think it would be to save this information in a database but that would more complicated.


texts: Creating buffer using ArcGIS API for JavaScript?

texts:
How can I create a buffer that is based on the shape of the geometry that has been drawn?

My map application creates a default square buffer extent and I want to create a buffer round a circle based on the shape not the min and max x y coordinates:

if (geometry.type === "polygon" || geometry.type === "multipoint" || geometry.type === "polyline") {
    var searchExtent = geometry.getExtent();
    var biggerExtent = new esri.geometry.Extent(searchExtent.xmin - 500, searchExtent.ymin - 500, searchExtent.xmax + 500, searchExtent.ymax + 500, new esri.SpatialReference({
                wkid : 27700
            }));

Can I give the esri.geometry.Extent a different searchExtent?


texts: Geocoding accuracy of Google Maps vs Tiger

texts:
I recently compared result of US Census Geocoder(http://geocoding.geo.census.gov/geocoder/) with that of Google Maps. For any given address both of them are giving different set of lat, lon. Why there is a difference between both of these?

My assumption is that Google Maps are based on US Census Geography data so both should be same.


texts: Installing QGIS3 on a Mac with Python in Anaconda

texts:
I am trying to install QGIS3 via Kyngchaos installer. But as I install python 3.6 in anaconda, the installer cannot find it. I know I can install a new environment or use symlinks to connected to Python that installed by homebrew as suggested in Installing QGIS3 on Mac? But as I have already installed a lot of packages in anaconda environment, I still hope to connect QGIS libs to it. Is it possible to do that?

Python 3.6 is installed in /Users/.../anaconda3/bin


texts: Viewing ESRI File GDB features without installing software

texts:
I'm looking for a way (Open source or otherwise) to view the contents and the features in a File GDB without installing any software locally.

Are the any current websites or web/map services or online viewrs that enable the uploading and viewing of a .gdb?

Again, I can obviously use QGIS or ArcGIS Explorer, but I need a solution for a local computer network that doesn't allow external installation, but does have internet access.


texts: Android call WFS GetFeature Service returns empty

texts:
I am trying to call a WFS GetFeature service to "implement" onClick get feature attribute popup (info window) but I am getting empty JSON. Here is WFS link that I am calling:

http://geo-polis.gistemp.com/geoserver/wfs?request=GetFeature&version=2.0.0&typeName=geopolis:zwnh_ymhttou&BBOX=37.896877%2C23.79923%2C37.89999%2C23.799334,EPSG:4326&OUTPUTFORMAT=json&service=WFS&

The response that I am getting is:

{"type":"FeatureCollection","totalFeatures":0,"features":[],"crs":null}

Update:

I want to implement the onClick GetFeatureInfo feature on Android, is this right service to use? or is there any other method (WMS) or some thing that I can use?


texts:
I have a Sentinel-2A and I can't fully classify it beacause it's cloudy.
I watched the Luca congedo YouTube tutorial but it is applied to the Level-1C data. And the cloudmasking plugin in used on Landsat data only.

How can I remove/ filter those clouds?


texts: How to import a digital elevation model in Germany Gauss Krüger from .dat file?

texts:
i am very new to qgis.
I want to create a hillshade file and also calculate the inclination.

The structure of the digital elevation model data looks like this:

5390.0125 5659.9875 251.8
5390.0125 5659.9625 250.6
5390.0125 5659.9375 247.5
5390.0125 5659.9125 240.6
5390.0125 5659.8875 236.1
5390.0125 5659.8625 234.8
5390.0125 5659.8375 235.6

The file format is *.dat. The coordinate system is Gauss Krüger, Bessel, 5. meridian.

I have no idea how to import the data in qgis.


texts: How to link data to a polygon ArcMap 10.3

texts:
I have a map with different polygons (118). Each of this polygons represent diferent vegetation types.

For each of those polygons I have data of which species there is and how abundant they are (the abundance is measured with a scale from 1 to 3).

I would like to "link" (I am sorry if its not the correct way of name it) the data of which species and abundance there is in each of those polygons.

I would like afterwards to be able of selecting the different polygons containing the same abundance and species.

In the polygons I created there is no atributes associated other than (FID, shape and Id. the 3 of them created alone by the program when I drew them) and the tables I have in excel with the data I want to "link" to the polygons contain a column called species and another abundance.

As you can see, I'm pretty new using the program...

As I have tried to do what you all suggest in the comments down here, and I couldn't achieve what I was wanting to do (provably because of my bad explanation before...)

I leave a picture:

As you can see in the attribute table, that is how every attribute table of each polygon (which represents one area) looks like. On the right there is how I have my data in excel.
Could you please suggest something I could do in order to add to the polygons attribute table the data belonging of species and its abundances (polygon number 1 has to have all the species (Acer mons, Acer pseu, Acer repens….) and its abundance (which are the numbers of the row Area1).


texts:
I work on Windows; I use Christoph Gohlke wheels; Python 3.7

After installing GDAL-3.0.1 and Fiona-1.8.6

an error occurs on import of fiona
from fiona.ogrext import Iterator, ItemsIterator, KeysIterator
ImportError: DLL load failed:

by downgrading to GDAL-2.4.1 I could surpass the error, yet rasterio starts rasing an error :
from rasterio._base import gdal_version
ImportError: DLL load failed:


texts: Why meridian time zones are not uniform?

texts:
Looking at this image:

Why are time zones are not uniform? (I mean not really parallel !!).

For example, time zone -11 and -10 are pretty overlapping near Alaska and time zone 0 changes the form to grab Island and so on ....


texts: How to mosaicking image that doesnt align each other

texts:
I need to mosaic a couple of images (satellite, 0.5*0.5 m high resolution). After mosaicking, I will do rectification using my own GCP from the field.

For the example, I have image 1 and image 2, I want to put image 2 to covering cloud area in image 1 (logically image 1 has a bigger area than 2).

My problem is, the image doesn't align with each other. Therefore,

My Question 1:
What is the best solution, should I choose option #1 = doing "image to image" correction from image 2 to image 1, then do mosaicking and after that geometric correction or option #2 = doing geometric correction first for each image using my own GCP then do the mosaic process. Which one has the best result?

My Question 2:
If I choose option #1, a new problem occurs. When I want to align visually both images I get high RSME on the table (which is my client want low RSME), but if I stay with low RSME I get a bad visual image, here is my screenshot when doing "image to image".

I have been trying use less or many points from image 1 act like GCP for image 2 to do "image to image" correction, been trying all transformation (1st, sec, third polynomial, spline and others) but still can't get good visual and low RSME, does anyone know how to fix this?

Moreover, if I choose option #2, which is a bad solution considering each image have a different distortion, is not it?


texts: Update PostgreSQL table; insert data from subquery, on conflict do update duplicate id

texts:
Context
I run into (some good old) troubles when trying to update a PostgreSQL (10.12) table called nodes containing nodes with some few attributes from OSM + some few others.
I only want to update the OSM attribute (lat, lon, osmid, highway, geom) but not the few others I also have on the table nodes.
In order to update this table nodes without breaking it, I took a security margin and created an other clean table called osmnodes where I first fetch the OSM basic data I need (see above) using ogr2ogr (...do some pre-processing...) and, only after that, I upgrade my main table nodes as follows:
INSERT INTO schema.nodes
   (id, lat, lon, osmid, highway, geom) -- I only want to upgrade these basic OSM fields, not the others
OVERRIDING SYSTEM VALUE -- this is mandatory, see the note hereunder
SELECT
     b.fid, b.x, b.y,
     b.osmid::bigint, b.highway, b.geom
FROM osmschema.osmnodes b;

Note that id of table schema.nodes was initially defined as: id INTEGER PRIMARY KEY GENERATED ALWAYS AS IDENTITY.
Error 1
This obviously sends me a duplicate key error:
ERROR:  duplicate key value violates unique constraint "nodes_pkey"
DETAIL:  Key (id)=(1) already exists.

So I naturally went for a coupled update such as:
INSERT INTO schema.nodes
   (id, lat, lon, osmid, highway, geom)
OVERRIDING SYSTEM VALUE
SELECT
     b.fid, b.x, b.y,
     b.osmid::bigint, b.highway, b.geom
FROM osmschema.osmnodes b

ON CONFLICT (osmid) DO
UPDATE SET
(lat, lon, highway, geom) -- removed the id here as I want to keep its already present value
=
b.x, b.y, b.highway, b.geom
-- FROM osmschema.osmnodes b -- the line to be freed to get to the next step hereunder
WHERE b.osmid = nodes.osmid; -- update the matching id

Error 2
Then I get:
ERROR:  missing FROM-clause entry for table "b"
LINE 13: (b.x, b.y, b.highway, b.geom)
          ^

Error 3
And, If I add the FROM clause which is asked for (un-comment the line in the previous snippet):
ERROR:  syntax error at or near "FROM"
LINE 14: FROM osmschema.osmnodes b

Question
I'm missing something there (...), I thought I would be simple...
The OVERRIDING SYSTEM VALUE seems not to work after it's been called.
What's wrong with my statement?


texts: I just convert hdf file to tiff with MRT but my result is grayscale map. How to get MODIS data with 3 band and step to change hdf file to colored map?


texts: Are coordinates taken with an handheld GPS device inherently in ESPG:4326 WGS84 and reprojected when added to a Google basemap?

texts:
I am just getting started using QGIS 3.6.3. I added a Google satelite imagery using QuickMapServices, which is projected using EPSG:3857 - WGS84/Pseudo-Mercator.  Then, while importing some points in a .csv file, I thought that I should choose the same CRS so it lines up with the basemap.  But that puts my data on the other side of the world.  When I select ESPG:4326 - WGS 84, it lines up whith the basemap.

Are coordinates taken with an handheld GPS device inherently saved in ESPG:4326 and then reprojected when added to a Google EPSG:3857 - WGS84/Pseudo-Mercator basemap?


texts: 'ST_Subdivide' function (3.5.0+ required) GEOS

texts:
Im running a query with ST_SubDivide() and receiving the error bellow.

ERROR:  The GEOS version this PostGIS binary was compiled against (34) doesn't support 'ST_Subdivide' function (3.5.0+ required)

Im on centOS7 with:

$ rpm -qa postg*
postgresql95-9.5.4-2PGDG.rhel7.x86_64
postgis2_95-client-2.2.2-3.rhel7.x86_64
postgresql-libs-9.2.15-1.el7_2.x86_64
postgresql95-libs-9.5.4-2PGDG.rhel7.x86_64
postgis2_95-2.2.2-3.rhel7.x86_64

$ rpm -qa geo*
geos-3.5.0-1.rhel7.x86_64
geos-devel-3.5.0-1.rhel7.x86_64

Question : How should i interpret this error of 'The GEOS version this PostGIS binary was compiled against (34)' when on my view i have GEOS 3.5.0 installed , what are my options ?


texts: Using dictionary and cursors to calculate a field based on values from another table - KeyError [0]

texts:
I'm new to dictionaries and am still a Python beginner. I'm trying to update field(s) based off of field values in another table using the OID field as a key. The search cursor dictionary has been filtered down (in an effort to conditionally update the values in the Line_Count_13_8 field) with a where clause for the portion of the script where I am receiving an error. I am getting:

  KeyError: [0] error

on this line:

if search_feats_138[0] == upd_row[0]:

Here is the whole block:

import arcpy
from arcpy import env
import time
start = time.time()
print start
print "Start"
arcpy.env.OverwriteOutput = True

subLayer='C:\Topo_Check_Tess_V5.gdb\Sub'
transLayer='C:\Topo_Check_Tess_V5.gdb\TL'
ppLayer='C:\Topo_Check_Tess_V5.gdb\PP'

pp_trans_GUJoin ='C:\Topo_Check_Tess_V5.gdb\pp_trans_GU_SpatialJoin'

SubLayer_pp_trans_GUJoin = 'C:\Topo_Check_Tess_V5.gdb\SubLayer_pp_trans_GUSpatialJoin'

arcpy.SpatialJoin_analysis ("SubLayer_Layer", "pp_trans_GUJoin_Layer", SubLayer_pp_trans_GUJoin, 'JOIN_ONE_TO_ONE', 'KEEP_ALL', fieldMappings, 'WITHIN_A_DISTANCE', .000002)

THIS SECTION WORKS

#add in Line_Count field
arcpy.AddField_management(subLayer, "Line_Count", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

search_feats = {f[0]:f[1] for f in arcpy.da.SearchCursor(SubLayer_pp_trans_GUJoin,["OBJECTID","Join_Count"])}

with arcpy.da.UpdateCursor(subLayer,["OBJECTID","Line_Count"]) as upd_cur:
    for upd_row in upd_cur:
        upd_row[1] = search_feats[upd_row[0]]
        upd_cur.updateRow(upd_row)
    print "Line_Counts Updated"
del upd_cur

THIS SECTION DOES NOT WORK - KeyError: [0]

#add in Line_Count_13_8 field
arcpy.AddField_management(subLayer, "Line_Count_13_8", "SHORT", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

#where clause
where_13_8 = ' "VT" = 13.8 '

#arcpy.MakeFeatureLayer_management(SubLayer_pp_trans_GUJoin, "Sub_PP_Trans_Join_13_8", where_13_8)

search_feats_13_8 = {f[0]:f[1] for f in arcpy.da.SearchCursor(SubLayer_pp_trans_GUJoin,["OBJECTID","Join_Count"], where_13_8)}

with arcpy.da.UpdateCursor(subLayer,["OBJECTID","Line_Count_13_8"]) as upd_cur:
    for upd_row in upd_cur:
        if search_feats_13_8[0] == upd_row[0]:
            upd_row[1] = search_feats_13_8[upd_row[0]]
        else:
            upd_row[1] = 0
        upd_cur.updateRow(upd_row)
    print "Line_Count_13.8 Updated"
del upd_curl

end = time.time()
print (end - start)/60

What should I change?


texts:
I need a tool for converting norweigan SOSI files to ESRI Shape using a Mac. It could also be a command line tool for Windows.

I tried to build GDAL with SOSI support on my Mac following instructions on this page but when I try compiling fyba I get this error message:

INQTID.cpp:18:11: fatal error: 'sys/vfs.h' file not found

Any ideas?


texts: Table join erases previous field calculations

texts:
I have issues with a normal join of two layers. The first set consists of a large list (n=259102) and I need to pick out 56 objects from this list. For that I have a second list with the 56 object IDs which are consistent with an ID field in the first list, as is pointed out here.

My first issue was that the two fields containing ID was in different formats, so i have tried changing the type in my smaller list using the field calculator. The way I did this is just creating a new field with the right type and then copying my ID in to that new field, as is suggested here.

But now to the issue. When i try to join these to layers I only get null fields in the large list and my objects are not picked out. When i look in the smaller list, my new IDs are all null as well. This was not the case when I first looked at the attribute table and it seemed to have changed during the join. I have saved the edit before i did the join, so that should not be the issue.

Is this an issue that anyone else has encountered and has a solution to?

I'm using QGIS 2.18.13.


texts: Removing duplicate highway shield symbol labels from ArcMap document

texts:
I have highway shield symbols in my map document that are showing duplicate labels in the current map view. When I go to the Properties of the roads layer, then Labels, then Placement Properties and select Remove Duplicate Labels and click OK, the duplicate labels do not disappear.

Also to note, I am using the Standard Labeling Engine and not the Maplex Labeling Engine and have to use the Standard Labeling Engine due to Maplex not being supported by our vendor application where the mxd is loaded.


texts: Selecting latest date with QGIS field calculator?

texts:
My shapefile contains columns like "Street", "number" and "date".

What I want to do is, to find the latest date (where we have done house services). If we have 2 dates for one house number ("number"), we want to get rid of the older one (this feature should be deleted). I also added a column "duplicates" which shows 1 when we have more then one date for a house number.

Has anyone an idea how to select the latest dates (or how to delete the oldest ones) with the field calculator?


texts: Can I determine hydraulic load on river bank with hec-ras?

texts:
I want to calculate the force exerted by wave action on river bank.

If I know the discharge, flow rate and water level, can I calculate the hydraulic load acting on river bank with hec-ras?

If not, how can I do this?


texts:
In my Openlayers 3.9.0, I use Bing Maps and OSM layers. They both are EPSG 3857.

var layer = new ol.layer.Tile({
        source: new ol.source.OSM({})
    });

    var bingMapsAerial = new ol.layer.Tile({
        source: new ol.source.BingMaps({
            key: '////',
            imagerySet:'AerialWithLabels'
        })
    });

I set the extent that the layers mast fit in

var textent = ol.proj.transformExtent([2297128.5, 4618333, 2459120.25, 4763120], 'EPSG:900913', 'EPSG:3857');

I also set a layer from PostGIS

    var ait = new ol.layer.Tile({
    extent: textent,
    source: new ol.source.TileWMS({
      url: 'http://localhost:8080/geoserver/mymap/wms?',
       params: {'LAYERS': 'mymap:polygon, mymap:line', 'TILED': true,  'VERSION': '1.3.0','FORMAT': 'image/png' ,'CRS': 'EPSG:3857'},
       serverType: 'geoserver'
     })
   })

My view is

var view = new ol.View({
    center: center,
    extent : textent,
    zoom: 6,
    maxZoom:20
});

and my map

var map = new ol.Map({
    target: 'map',
    layers:[bingMapsAerial, layer, ait],
    view: view
});

map.getView().fit(textent, map.getSize());

The raster tile layer is in EPSG 900913 in PostGIS and Geoserver says EPSG:900913 as Native SRS and Declared SRS. Both polygon and line layers that compose the tile layer have slightly different Native Bounding Box and Lat/Lon Bounding Box, but they render the same spatial region.
Both polygon and line layers have Published zoom levels and Cached zoom levels set to 0 -20 in Geoserver to match the Openlayers View settings

Aren't EPSG 900913 and EPSG 3857 supposed to be the same?

Both Bing and OSM layers look ok, but the tile layer appears smaller than the other two.

How do I proceed? Am I missing something, a setting or var? Is there a hack? Do I have to convert the PostGIS data?


texts: How to import and display a road network in XML to ArcGis?

texts:
I am trying to convert xml data to shapefile features using ArcGis 10.1.

Im using the open source SUMO traffic simulation software. Part of the model process is to convert OpenStreetMap road network data to an readable xml format.

The xml format describes the "lanes" between nodes in the network. The below example is for a straight horizontal road/edge with 2 lanes between nodes "A" and "B". The length is in meters, and the "shape" give the start/end X and Y co-ords on a basic XY grid.

<edge id="AtoB" from="A" to="B" priority="-1">
        <lane id="AtoB_0" index="0" speed="13.90" length="50.00" shape="0.00,51.65 50.00,51.65"/>
</edge>
<edge id="BtoA" from="B" to="A" priority="-1">
        <lane id="BtoA_0" index="0" speed="13.90" length="50.00" shape="50.00,48.35 0.00,48.35"/>
</edge>

Ideally I would be able to generate shapefile features for each lane in the network. Also the above example is for a test I am conducting and is not spatially referenced, for the end network I will use WGS 84.


texts: Splitting polygon into separate polygons using table of individual lines

texts:
I have been trying for a while but have not been able to split a polygon at multiple locations, using a table of lines, to get what is on the right side of the above image.
The best I can achieve is the large polygon being split at one of the locations. My query is below:
drop table if exists new_polygons;
CREATE TABLE new_polygons AS
SELECT ST_Split(a.geom, b.geom)as geom
FROM large_polygon as a, blades as b
WHERE ST_Intersects(a.geom, b.geom);


texts: How to view the class files of QuantumGIS for further development?

texts:
I have been searching class files of QGIS for a while for customisation and edit it for a project..But I am unable to get or download the files...can anybody give some references if we can really view the files or its not availabe..
Thanks


texts: How to create an rgb GeoTIFF file raster from bands using the GDAL python module?

texts:
I'm trying to create a new raster from three bands (rgb). I got the script mostly from here. I can create a single band raster with no problems, but I'm having trouble inputting a tree band numpy array. I'm greeted with this error message when I attempt it.

  ValueError: expected array of dim 2

If I try to create a three band GeoTIFF image I end up with a grey,alpha,alpha image instead of an rgb image, nor am I able to populate these with my three-dimensional numpy array. The following script works for creating a greyscale raster. with three bands as grey,alpha,alpha.

from osgeo import gdal, osr
import os
import io
import numpy as np

def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array,format,spatialRef):

    bands = array.shape[0]
    cols = array.shape[1]
    rows = array.shape[2]

    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName(format)
    outRaster = driver.Create(newRasterfn, cols, rows, bands,gdal.GDT_UInt16,options = ['PHOTOMETRIC=RGB', 'PROFILE=GeoTIFF',])
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

    for band in range(bands):
        outRaster.GetRasterBand(band+1).WriteArray(array[band,:,:])
        outRaster.GetRasterBand(band+1).FlushCache()
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(spatialRef)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())

def get_rgb_bands():
    source = gdal.Open("C:\\RGBtest.tif")

    red = source.GetRasterBand(1).ReadAsArray()
    green = source.GetRasterBand(2).ReadAsArray()
    blue = source.GetRasterBand(3).ReadAsArray()

    #rgbOutput = source.ReadAsArray() #Easier method
    rgbOutput = np.zeros((3,16,16), 'uint16')

    rgbOutput[0,...] = red
    rgbOutput[1,...] = green
    rgbOutput[2,...] = blue

    #Clear so file isn't locked
    source = None
    return rgbOutput

exampleRGB = get_rgb_bands()

exampleOutput = "C:\\ExampleOutput.tif"
rasterOrigin=[0,0]

array2raster(exampleOutput,rasterOrigin,1, 1,exampleRGB,'GTiff',3879)

createdFile = gdal.Open("C:\\ExampleOutput.tif")

print 'old file info:'+gdal.Info(source)
print 'new file info:'+gdal.Info(createdFile)

createdFile = None

exampleOutput = None

EDIT: Luke answered this question perfectly I have updated the code to reflect a working example to output a 16-bit rgb raster.


texts:
I am new to GIS domain. Installed ArcGIS Desktop 10.7.1 on Windows 10. There is no "Add basemap" option even after proper internet connection.

What version of .NET framework do I need to have in my system for this?


texts: SQL Server - Select all polygons nested within a larger polygon

texts:
This is a seemingly simple SQL Server geometry question that I thought would have an out-of-the box solution, but I'm not having any luck finding one.

My intent is to select all the records within one table that have polygons that are nested (contained) within a larger polygon from another table. I had expected functions STWithin and STContains as being the solutions I needed, but unfortunately, both identify only the inner polygons within those nested in the larger polygon, not those nested polygons that are touching the boundary of the larger polygon. See image for example.

An alternative option that worked to my needs was STIntersection. The problem with this function, however, is that it only returns the geometry column! I would like to obtain the record ID instead. Does anyone have a suggestion on how this can be done?

STWithin:

select a.bg10 from
gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STWithin(b.shape) = 1
where b.mktname = 'Loop'

STContains:

select a.bg10 from
gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on b.shape.STContains(a.shape) = 1
where b.mktname = 'Loop'

STIntersection:

select a.shape.STIntersection(b.shape)
from gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STIntersects(b.shape) = 1
where b.mktname = 'Loop'

Edit:

One suggestion was to omit STIntersection and use solely STIntersects as follows:

STIntersects:

select a.bg10
from gis.usa_10_block_group a
join gis.usa_10_mkt_definition b
on a.shape.STIntersects(b.shape) = 1
where b.mktname = 'Loop'

The problem with this approach, is that STIntersects appears to select all polygons either within or outside and touching the larger polygon, not just those strictly within. See image for example.


texts: Update geometries with ST_UnaryUnion at TRIGGER

texts:
We are currently working on database with multipolygons and we're using QGIS 3. Sometimes it's very useful to work with parts (each page of document is separate part of multipolygon), however at the end I need as few parts as possible.
When I tried to update geometries with St_UnaryUnion I had some problems and finally ended with 'temporary workaround' like this:

update workspace.tests
set geom = geom2 from (
select ST_Multi(ST_UnaryUnion(ST_Multi(geom))) as geom2
from workspace.tests
where id_number = 23) as ge
where id_number = 23;

It works fine, but I have to control each row separately (or write some kind of  loop). So I read many stackexchange topics and tried to write my first trigger.

CREATE OR REPLACE FUNCTION workspace.mergpoly() RETURNS trigger AS $mergpoly$
    BEGIN
        NEW.geom := ST_Multi(ST_UnaryUnion(ST_Multi(OLD.geom)));
        NEW.howmanyparts := ST_NumGeometries(NEW.geom);
   RETURN NEW;
    END;
$mergpoly$ LANGUAGE plpgsql;

CREATE TRIGGER mergpoly AFTER INSERT OR UPDATE ON workspace.tests
    FOR EACH ROW EXECUTE PROCEDURE workspace.mergpoly();

The problem is: it does nothing. I tried also other ways but I failed every time. I'm quite new to PLSQL and PostgreSQL. Am I close to good answer or would you write it completely different?


texts: How to specify pixel resolution when converting feature collection to image in Earth Engine

texts:
I have a feature collection and I would like to convert it into a raster image. I would like the resulting image to have 30m spatial resolution.
How can I use the reduceToImage method to assign resolution when saving the resulting image as an asset?

Here's my script:

var table = ee.FeatureCollection("users/mattelearning85/exe/2009_shapefile");
print(table)
var table_raster = table.reduceToImage(['Area_ha'], 'mean');
print(table_raster)

Code Editor link


texts: Letting user draw AOI or draw points for tool input In ArcMap script tool?

texts:
I have a python script which uses a point file to do some spatial analysis. After opening the tool in arcmap, I would like the user to be able to just click points on the map to generate the data. I found something about "Interactive Feature Input", which I believe is what I want. However, the information I found was pretty much useless. Currently, the tool just uses an existing point file.

I use ArcGIS Desktop 10.4 or 10.6.


texts: When I am trying to resample a raster layer or to clip it by mask layer the pixel value get automatically changed in the output file irrespective of every setting. Why it is happening and how can i get rid of that? I am using QGIS.


texts: OGR_LINK_FIELD to filename in a Geospatial PDF

texts:
I'm using gdal_translate to create a Geospatial PDF. I'm trying to add file path to my feature when clicking on them by using:

-co OGR_LINK_FIELD=<column with path>

I am creating my PDF in the main folder and the files that I link to are in main/images/unique.jpg

In my VRT layer is a column with "Filepaths".

My code is:

gdal_translate -of PDF \main\input_background.tif \main\output.pdf -co OGR_DATASOURCE=\main\data\vrt_input.vrt -co OGR_LINK_FIELD=Filepaths -co TILED=YES

In the shapefile, used to create the VRT is column Filepaths filled with:

/images/unique.jpg

In the resulting PDF I see that file:// is added to the link, but when I click on it it open the browser, but doesn't show the image.

If I place the output PDF in the same folder as the images and set the filepath to

unique.jpg

it opens the image.

What can be wrong?
Or is it not possible to use relative paths in this?


texts: Can rasters in same CONSTRAINED table have different alignments in PostGIS?

texts:
Context

I have a table filled with ASTER GDEM data for the state of Colorado tiled at 100x100 pixels. I want to query the table and return a single raster that covers all of the supplied bounding box (it surrounds Fort Collins, in this case).

Then, I export that raster in the selected GDAL format using a custom function written in Python/PL called write_raster.

I am attempting to do so using this query:

SELECT write_raster(ST_AsGDALRaster(ST_Union(rast), 'Gtiff'),  -- Raster bytea to write
  '/home/nronnei/gis/data/postgis-raster/colorado_elevation.tif'::text  -- Path to write to
)
FROM public.co_dem_tiled
WHERE ST_Intersects(
  ST_MakeEnvelope(-105.4330444336, 40.3789822705,
    -104.7024536133, 40.7170785158, 4326),
  ST_Envelope(co_dem_tiled.rast)
);

I'm using PostGIS v2.1, and the 'same_alignment' constraint is not set on the table. When I try to add the constraint it fails because the data aren't aligned.

The Problem

When I run this query, I get the following output:

ERROR:  rt_raster_from_two_rasters: The two rasters provided do not have the same alignment

Everything runs just fine if I remove ST_Union from the query above. What really makes this weird is that I'm only pulling raster from a single, constrained table.

Naturally, I tried using ST_NotSameAlignmentReason like so:

WITH base AS (
  SELECT rast
  FROM public.co_dem_tiled
  WHERE ST_Intersects(
    ST_MakeEnvelope(-105.4330444336, 40.3789822705,
      -104.7024536133, 40.7170785158, 4326),
    ST_Envelope(co_dem_tiled.rast)
  )
),
top AS (
  SELECT rast
  FROM public.co_dem_tiled
  WHERE ST_Intersects(
    ST_MakeEnvelope(-105.4330444336, 40.3789822705,
      -104.7024536133, 40.7170785158, 4326),
    ST_Envelope(co_dem_tiled.rast)
  )
)
SELECT ST_NotSameAlignmentReason(base.rast, top.rast)
FROM base, top;

However, this doesn't give me anything useful. It simply reports that all rasters have the same alignment, as expected.

What's going on here?!


texts: Using R to calculate the area of multiple polygons on a map that intersect with another overlaid polygon

texts:
I have a shapefile downloaded from the Ordnance Survey that gives electoral ward (division) boundaries for a county of the United Kingdom. I have successfully used R to load the shapefile and plotted various maps using ggplot2 as described in this question. It's all working rather well.

Now I would like to create a new polygon of arbitrary shape, add it to the map, then calculate the population living in the area lying under the shape, which might cover or partially cover multiple divisions. I have the population for each electoral division and I can make the simplifying assumption that the population in each ward is uniformly distributed. That suggests the following steps.

1) Overlay a new shape on the map that partially covers multiple electoral divisions. Let's say there are 3 divisions, for the sake of argument. It would look something like this. [Edit: except that in the image below the shape straddles 5 divisions rather than 3]

2) Calculate the percentage of the area of each of these 3 divisions that intersects with the overlaid polygon.

3) Estimate population by getting the percentage of the area of each division covered by the overlaid shape and multiplying this by the population of each division.

I think I can probably work out how to create the polygon and overlay it on the map i.e. add it to the existing data frame using the useful answer to this and other questions. The bit that worries me is the task of working out the percentage of each division that is covered by the overlaid shape. The lat and long columns in the data frame are those strange Ordnance Survey OpenData figures (Eastings and Northings or something).

So my first question is: How would I go about finding the area (or a subset of the area) of the polygons that define the borders of an electoral division using this data? Because even a meaningful subset of this data frame is large I have used dput to create a 500k file (which can be copied and pasted or downloaded from here) rather than posting it in this question. The map that forms the base for the image above was created with the following:

require(ggplot2)
ggplot(smalldf, aes(x = long, y = lat, group = group)) +
    geom_polygon(colour = "grey50", size = 1, aes(fill = smalldf$bin))

My second question is: am I using the right tools? Currently I am using readShapePoly from the maptools package to read the shapefile. I then use fortify to create a data frame of about 130k lines, suitable for use in ggplot. Maybe I should be using a different package if there is one with useful tools for such processes?


texts: Calculating a proper region with insets for my MapView

texts:
In my app I'd like to create a snapshot of a Map and as I have a custom tile-source I need to create that on my own.
I am currently stuck creating a proper region that fit's a track I want to display.

So my given values are:

width/height of the resulting image
the region (coordinates) of the track
insets (top, bottom, left, right) I want in px

So now I need to know how to calculate a region for my needs. So basically northwest / southeast coordinates for that region.

Update:

So the end-result would be a picture of a map (using my own tiles) that would show a track a user went. And I'd like to know how to exactly calculate the bounds I need (currently I do some rule of thumb calculation).

So as a sample: I want to have an image in the end, that is 500x500px. I want there to be at least 50px of margin on each side and the track should be centered inside.

What I know are the coordinates of my track (lat/lng). So I can calculate a bounding box.

Now I need to figure out a proper zoom level (already know how to do that) but I don't know how to say: In the final image the coordinate on the top left will be this, and the on on the bottom right will be that (which would be great, cause then I would EXACTLY know which tiles to download)

Currently I hat north/east and south/west coordinates. I take the center, and use max(lat-span, lng-span) as the size. Then I take this span 1.5 times and have my coordiantes. This seems to give me reasonable images, but I was just wondering how to really calculate it.

On iOS, a MKMapView has the option to say: setVisibleMapRect:edgePadding: So basically say: this is my region and I want to have at least the given padding. So I basically need to know how this MapView calculates stuff. :) But I am not a gis guy and don't know much about calculating stuff like that :)


texts: What does "Do you want to reset the field mapping" mean in QGIS Refactor Fields dialog?

texts:
I want to convert the type of some fields in a vector layer from String to Real. I am trying to use the Processing -> QGIS geoalgorithms -> Vector table tools -> Refactor fields tool to do this. When I open this up and select the layer I'm trying to edit, a dialog pops up that asks "Do you want to reset the field mapping?" (Yes/No) ... I can't find documentation anywhere on what this means.

Can someone explain?

(See lower center in image for dialog I'm referring to)


texts: Converting set of points to polygon using C#?

texts:
I want to convert a set of gps ponints to polygon.I used a simple loop to convert it to a WKT and than used makevalid() in sql server and save it in a gemetry column in databasae.Is this way right?or I missed something.
for example I have this points

<x1,y1>,<x2,y2>,<x3,y3>,<x4,y4>

I used a loop to creat a string like

Geometry::STFromText('Polygon(x1 y1,x2 y2,x3 y3,x4 y4,x1 y1)',4326).Makevalid()

Is this method right?Do I need any convertin on points .


texts: Road layers appearing as points in GeoServer previews - both OpenLayers and GeoExp

texts:
I am wondering why my roads layer keeps appearing as points when previewed inside GeoServer - layer preview using both the OpenLayers and GeoExplorer options. The data are imported from PostGIS.

I can't figure out where the problem might be originating from as I am also new to both PostGIS and GeoServer. In PostGIS, I have added a Geometry table using:

CREATE TABLE geometries (name varchar, geom geometry);
INSERT INTO geometries VALUES
  ('Point', 'POINT(0 0)'),
  ('Linestring', 'LINESTRING(0 0, 1 1, 2 1, 2 2)'),
  ('Polygon', 'POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))'),
  ('PolygonWithHole', 'POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(1 1, 1 2, 2 2, 2 1, 1 1))'),
  ('Collection', 'GEOMETRYCOLLECTION(POINT(2 0),POLYGON((0 0, 1 0, 1 1, 0 1, 0 0)))');

SELECT Populate_Geometry_Columns();

SELECT name, ST_AsText(geom) FROM geometries;

However the geometry column in my tables is named "wkb_geometry". Could this be an error? Please help!


texts: File creation error, something wrong with bin, from QGIS?

texts:
I've been trying to use the Union tool on a layer of building polygons, with a layer of output area polygons. The building layer has the largest number of features I've used so far, 106716, if that's relevant.

Yesterday, QGIS was just freezing up trying to do the Union, and I ended up cancelling it with a view to giving it loads of time today.

However today I've been getting an error that reads:

  No output created. File creation error:  creation of layer failed (OGR
  error:Failed to open Shapefile `G:/New folder\Test.shp'.)

Can anyone help me resolve this new error?


texts: How can one retrieve all the raster cell values for Kernel Density output?

texts:
This question is related to the replies in this earlier thread (How to interpret GRASS v.kernel results?).

I am new to GIS mapping (arcGIS) and have been trying to work with the Kernel Density output. I have read on many threads that a good way to check Kernel Density output is to multiply the density values by the area of the cell and then sum those values over the grid. That sum should be equal to the number of points (in my case) in the original data.

To multiply the density values by the area of the cell: The only density values I see are when I 'Identify' a cell producing a single density value for a single cell and the density classification legend used to symbolize the raster in the Table of Contents. How do I extract the density values from each raster cell? (So I can multiply those values by the area).

Or should I be generating density values for each cell from the density classification legend presented in the Table of contents?


texts: Generating wind direction raster from wind's components?

texts:
I have a point shapefile, with each point representing the location of a weather station and storing values for wind SPEED and DIRECTION (in degrees).

My goal is to create a continous surface (raster) for wind direction by interpolating the two wind's components (X and Y) derived using the sine and cosine of the wind direction as described here: http://mmc2.geofisica.unam.mx/cursos/geoest/Articulos/Geostatistics/Non-Linear%20Surface%20Interpolations.htm:

X=H*sin(q)

Y=H*cos(q)

where H is the wind speed, q is the wind direction (in degrees).

I did some web searches, and I found a very good comment from @whuber (Generating wind direction raster using its U and V components in ArcGIS Desktop?), who was suggesting to use the ATan2 function.

Would the following formula (where v and u are the two wind's components, corresponding to the Y and X components in the first cited web source) suffice to achieve the mentioned goal?:

(ATan2(v, u) / (2*3.14159265) + 0.5) * 360


texts: Cannot run semi-automatic classification (SCP) on Mac

texts:
When I click the run button, a window appear with an Open Option. When I select it, the QGIS Software reopens over the one that's already open.
There's this message:
2020-12-17T13:47:15     WARNING    Traceback (most recent call last):
      File "/Users/rodrigocarvalho/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/SemiAutomaticClassificationPlugin/maininterface/landsatTab.py", line 640, in performLandsatCorrection
      self.landsat(cfg.ui.label_26.text(), o)
      File "/Users/rodrigocarvalho/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/SemiAutomaticClassificationPlugin/maininterface/landsatTab.py", line 331, in landsat
      LDNmList = cfg.utls.findDNmin(tPMDN, NoData)
      File "/Users/rodrigocarvalho/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/SemiAutomaticClassificationPlugin/core/utils.py", line 2006, in findDNmin
      o = cfg.utls.multiProcessRaster(rasterPath = inputRaster, functionBand = 'No', functionRaster = cfg.utls.rasterUniqueValuesWithSum, nodataValue = noDataVal, progressMessage = cfg.QtWidgetsSCP.QApplication.translate('semiautomaticclassificationplugin', 'DOS1 calculation'), deleteArray = 'No', parallel = cfg.parallelArray)
      File "/Users/rodrigocarvalho/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/SemiAutomaticClassificationPlugin/core/utils.py", line 4675, in multiProcessRaster
      manager = cfg.MultiManagerSCP()
      File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/multiprocessing/context.py", line 56, in Manager
      m.start()
      File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/multiprocessing/managers.py", line 547, in start
      self._address = reader.recv()
      File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/multiprocessing/connection.py", line 250, in recv
      buf = self._recv_bytes()
      File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/multiprocessing/connection.py", line 407, in _recv_bytes
      buf = self._recv(4)
      File "/Applications/QGIS3.10.app/Contents/Frameworks/Python.framework/Versions/Current/lib/python3.7/multiprocessing/connection.py", line 383, in _recv
      raise EOFError
     EOFError


texts: Removal of drainage channels (defined by steep slopes) from LiDAR DEMs

texts:
I have LiDAR DEM tiles (geoTIFF) that feature modern (19th century?) drainage channels but also remains of the pre-modern landscape. To perform further analysis of the latter, I want to remove the channels in the same way buildings and other off-terrain objects are removed from DSMs to generate DEMs.

The channels are well defined by steep slopes which otherwise don't occur on the relatively flat landscape, so that an algorithm should be able to do this. My Question is: Is there a tool/plugin for ArcGIS or Whitebox which is able to do it or will I have to program it myself?

LiDAR has a Resolution of 1m.


texts:
I am trying to ellipses around polygons I have in im map, in ArcGIS 10.4. This what I need in the further steps for displacement of polygons. These polygons are from a geologic map and these specific polygons are small ones which I need to displace by drawing ellipse around them. One may ask, why not circle? Well, there are very different shapes in the map such as very round, very narrow and long, short, thus ellipse would fit them well than circle.

Update:

I could use two ways to find centroid, first, from feature to point (with inside) and second, centroid from "Calculate Geometry".


texts: Force Python GDAL/OGR to flush changes to disk

texts:
I am trying to write a simple list of features and export it as a KML file, adding one extra attribute: in this case trafficability. However, ogr is failing to flush the data to disk when I run the script. Instead, the XML is clearly cut off mid-element. I have even tried explicitly calling SyncToDisk, but it has still failed to write the element. Are there any other calls I need to make to ensure that the full file is indeed written to the disk?

driver = ogr.GetDriverByName('KML')
output_file = 'traffic.kml'
driver.DeleteDataSource(output_file)

data_out = driver.CreateDataSource(output_file)
fldDef = ogr.FieldDefn("trafficability", ogr.OFTReal)
data_in = ogr.Open(osm_file)
in_layer = data_in.GetLayer(1)
srs = osr.SpatialReference()
dest_layer = data_out.CreateLayer('layer1',
                            srs = in_layer.GetSpatialRef(),
                            geom_type=in_layer.GetLayerDefn().GetGeomType())
feature = in_layer.GetFeature(1)
[dest_layer.CreateField(feature.GetFieldDefnRef(i)) for i in range(feature.GetFieldCount())]
dest_layer.CreateField(fldDef)
for feature in master_feature_list:
    feature_new = ogr.Feature(in_layer.GetLayerDefn())
    feature_new.SetFrom(feature)
    try:
        feature_new.SetField('trafficability', trafficability_dict[feature.GetField('osm_id')])
    except:
        feature_new.SetField('trafficability', -1)
    print feature.DumpReadable()
    assert(0 == dest_layer.CreateFeature(feature_new))
    dest_layer.SyncToDisk()
    data_out.SyncToDisk()
print("Syncing to to disk", data_out.SyncToDisk())
print('Done conversion!')

This is probably a really silly bug, but the documentation on OGR can be lacking at times.


texts: Can one parameter choice on Python toolbox tool set Display Names of other parameters?

texts:
I would like to be able to create a Python toolbox tool (or a Python script tool) where a choice made from the first parameter (Category in the images below) changes the Display Name of other parameters (which initialize as Sub Category 1 and Sub Category 2 in the images below).

These images are mockups.

After opening the tool I choose Category B for Category.

I would like that choice to change the Display Names of:

Sub Category 1 to Sub Category B1 (choosing Category A would set Sub Category A1)
Sub Category 2 to Sub Category B2 (choosing Category A would set Sub Category A2)

Something like below (where I have mocked up the Display Name changes) but with the first parameter green rather than red.

I tried tool validation of a Python script tool first, but seemed to be blocked by the read-only Display Name there. I then moved to a Python Toolbox tool instead because that adds some extra functionality and I was hoping might remove the read-only Display Name limitation.

My attempt to code this as a Python Toolbox tool is below, but it seems not to be working because displayName is read-only there too.

Is there any ArcPy way to meet the requirement that I describe above?

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Tool"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        # First parameter
        param0 = arcpy.Parameter(
            displayName="Category",
            name="category",
            datatype="String",
            parameterType="Required",
            direction="Input")
        param0.filter.list = ['Category A', 'Category B']

        # Second parameter
        param1 = arcpy.Parameter(
            displayName="Sub Category 1",
            name="subcategory_1",
            datatype="String",
            parameterType="Required",
            direction="Input")

        # Third parameter
        param2 = arcpy.Parameter(
            displayName="Sub Category 2",
            name="subcategory_2",
            datatype="String",
            parameterType="Required",
            direction="Input")

        params = [param0,param1,param2]
        return params

    def isLicensed(self):
        """Set whether tool is licensed to execute."""
        return True

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        if parameters[0].value == "Category A":
            parameters[1].displayName = "Sub Category A1"
            parameters[2].displayName = "Sub Category A2"
        elif parameters[0].value == "Category B":
            parameters[1].displayName = "Sub Category B1"
            parameters[2].displayName = "Sub Category B2"
        return

    def updateMessages(self, parameters):
        """Modify the messages created by internal validation for each tool
        parameter.  This method is called after internal validation."""
        return

    def execute(self, parameters, messages):
        """The source code of the tool."""
        return

I am using ArcGIS 10.2.2 for Desktop.


texts: Performing multiple regression of time series data using R?

texts:
I am working with MODIS time series data, and I want to stack a stack object to prepare my data for multiple regression.

I have AMSR (soil moisture) as a response variable and I have some explanatory variables such as NDVI and LST and Albedo for a period of 6 months (183 images for each variable).

I want to perform multiple regression analysis between these variables. For a one particular day I can easily use “stack” command to stack the variables for one day and then create data.frame from the stack object, but to apply to the all-time series data I could not find a solution. I guess I should stack the list.files for each variable first, and then stack the stack objects again.

However, this has not worked.


texts:
I am analyzing Germany's high voltage electricity grid. I extracted the data from OSM, imported it in QGIS and simplified its structure by hand, specifically creating single nodes where power lines are connected in contrast to having a lot of high voltage substations with complicated geometries. Ultimately, I want to export that network and use it in GAMS to run an optimization model on it. But long story short.
I have two .shp files, standard WGS 84 CRS.

a file containing the network as lines (about 250)
a file containing all the nodes (about 350) in the network (where the power lines intersect)

Now I want to extract a shortest distance matrix containing the distances between all nodes in file 2, having used file 1 as my "road" network.

I have read through a couple of answers suggesting python, GRASS or PostGIS with pgrouting but none of them worked for me. I'm new to Python, GRASS and PostGIS, which sure does not really help.

(almost forgot. Mac OS 10.7.5, GRASS 7.0, QGIS 2.8.1, PostgreSQL 9.0.13, PostGIS 2.1.7, pgrouting 2.0.0)


texts: How to get ESRI Metadata into GeoNetwork?

texts: Has anyone had success in transferring metadata created in ArcGIS 10.2 into GeoNetwork?  Currently I am creating metadata with ArcGIS using the Item Description menu item in ArcCatalog.  I would like to store and publish that metadata in GeoNetwork but I have not been able to get this to work.  I have tried to import that xml into GeoNetwork using a variety of formats and style sheets but I either get a validation error or incomplete information.


texts: What is the maximum size supported for ArcGIS Server Workgroup?

texts:
I have a database with a size of 200 gb and expect to grow near of a 1 Tb in two years.

It´s possible to install an ArcGIS Server Workgroup 10.2 with this size of information? I´m afraid not understand the information about the size supported for ArcGIS Server Workgroup I saw that was 4gb (very poor) but in  functionality matrix Esri said 4gb but show this:

  At any given time, the firm may be using 3 GB of vector data stored in
  SQL Server Express and configured for concurrent edits from three
  desktop seats, over 50 GB of imagery stored as TIFF files, and 15 GB
  of data stored in file geodatabases.

It´s possible that the capacity of the database is different for the maximum size of ArcSDE?


texts: How to get JSON from Geoserver using AJAX request

texts:
I am building a web mapping application using Leaflet and a sort of customized WFS-T API to place and retrieve features from PostGIS. I have the script for putting features in PostGIS tables down solid, and I've managed to publish my tables as WFS services using Geoserver. My problem--which seems simple enough--is how to import those vector features back onto my map. I figured it would be easiest to go with the JSON option Geoserver provides, but when I make the below request using jQuery:

    $.ajax({
    url: "http://localhost:8090/geoserver/cite/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=cite:markers&maxFeatures=50&outputFormat=json&callback=?",
    dataType: "jsonp",
    success: function() {
        //
    }
}).done(function ( data ) {
    console.log(data);
});

or:

    $.getJSON("http://localhost:8090/geoserver/cite/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=cite:markers&maxFeatures=50&outputFormat=json&jsoncallback=?", function(data){
     console.log(data);
});

I get the following error:

SyntaxError: invalid label [Break On This Error]
{"type":"FeatureCollection","features":[{"type":"Feature","id":"markers.1","geom...

Now, I know that these two functions are equivalent; what would seem to be a mismatch between specifying json in the url string and jsonp as the datatype is in fact correct because Javascript interprets the response as Javascript rather than true JSON. I've also figured out from reading other forum posts that the problem is the json script is not being properly contained in an object variable; the syntax error is due to missing parentheses on either side -- ({jsonobject}) -- like so. My question is, how the heck do I fix the syntax to get JavaScript to work with it? Obviously everything I want is coming through--it's just breaking because of a stupid syntax error that I can't figure out how to overcome. Any help on this would be much appreciated.


texts: st_union of set of polygons result in polygon with very small holes

texts:
I have a problem with a polygon dataset. Due to the way it is produced (union of small Voronois) there are not always nodes at the same position in the common border between two adjacent polygons. When I do a st_union very small holes (can be less than 1 millimeter) appear. (They appear to be lines, but actually are polygons). This dataset is exported to other GIS environments (MapInfo) and they have the same problems with a geographical union.

I know that the polygons themselves are topological correct. But the constraint 100% adjacent polygons is not met. So my aim is to clean up this dataset to create a 100% adjacent polygon dataset. But alas, so far no succes.

Some remarks:

I use PostgreSQL 12 and PostGIS 3.0 and the improved precision of PostGIS 3.0 helps (on screen they look perfect), but does not solve this issue.
Within PostGIS I can do a st_union with a small buffer. That resolves my problems, but I make this for users with other GIS tools. So that's no solution.
I tried the topology extension, but importing failed. due to the lack of adjacency I think based on the error messages.
I tried st_simplify, st_snap, st_removerepeatedpoints and st_snaptogrid, and they all solve some issues, and produce other issues. So they are not robust solutions.
I also tried to replace the geometry with the st_difference of the geometry (with a small buffer) with the surrounding geometries, and that worked for a big part, but not for a 100%. Due to the buffer mainly, but without the buffer it doesn't work.

So I am a little lost on how to proceed. My last resort is to use ArcGIS or FME, but it would be nice to solve this in PostGIS.

I would love to have a function like st_snapCommonBorder( geom1, geom2, precision ) which creates nodes in the common border of geom1 where there are nodes in the common border of geom2, within a certain precision. (st_snap doesn't create nodes, and st_difference doesn't resolve holes)

Or is there another (robust) solution?


texts: Error importing a map with osm2pgrouting

texts:
I am trying to import an Italian osm file into PostGIS with osm2pgrouting. I launched the script and it worked until it encountered an error. This is the output of osm2pgr:

.....
Nodes table created
2create ways failed:
Types table created
Way_tag table created
Relations table created
Relation_ways table created
Classes table created
Adding tag types and classes to database...
Adding relations to database...
Adding nodes to database...
Adding ways to database...
ERROR:  extra data after last expected column
CONTEXT:  COPY ways, line 681382: "934378       117     0.15808824934990467    11.1369311       45.8187991      11.1388487      45.818314100000002      37264295srid=432..."
Creating topology...
NOTICE:  CREATE TABLE will create implicit sequence "vertices_tmp_id_seq" for serial column "vertices_tmp.id"
CONTEXT:  SQL statement "CREATE TABLE vertices_tmp (id serial)"
PL/pgSQL function "assign_vertex_id" line 14 at EXECUTE statement
#########################
size of streets: 1072190
size of splitted ways : 2281124
finished

I used the italy.osm.bz2 file decompressed finded here:http://download.geofabrik.de/osm/europe/
but have also tried with the one here:
http://download.gfoss.it/osm/osm/
I have tried different times, with different files -- in fact, I have been trying from almost a month, but nothing has changed: I always get the same error. Any suggestions?

Thanks
Luca


texts: CSS and border transparency around web map markers

texts:
The issue is that the image that I have in a web map visualizes map markers and place icons with white space around them.

I checked this in the FireFox browser developer console and it seems to be an issues with the background image/color being turned on. I click off the background image/color and the white space does indeed disappear but only in the browser developer console. I have attempted to modify my code in Brackets in just about every section to make the background color=transparent, but no luck. When I review the map in Google, there is no white space around images, just seems to be an issue with FireFox. I have opened the map in my smartphone as well as other devices and the white space remains. I have checked out help on www3 schools, but have been unable to resolve. Any thoughts?


texts:
I am trying to select all categories of shops for a portion of the UK.
To select just the supermarkets I would use "shop" = "supermarket", and so forth, while to select all categories I have tried "shop" = "*" but no results are returned.
Do you know how the below statements can be fixed to select all categories of shops? thanks

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“shop”
*/
[out:json][timeout:25];
// gather results
(
// query part for: “shop”
node"shop" = "*";
);
// print results
out body;

;
out skel qt;


texts:
On shapefile data, I had used center point extractor using point mode center and add StringConcatenator with result string 2201@Value(_inside_x)@Value(inside_y).

After that I get the result 2201158707510.139225325713360.34431
where _x value is 158707510.139225 , _y value is 325713360.34431 and 2201 is constant value.

I want to remove after decimal values in the result and want result as 2201158707510325713360.

Any way to do this?


texts: How can I apply correct colours to an elevation raster?

texts:
I have a problem to display elevation data with pseudocolours. The channel 1 is from -3.40282e+38 to 149.717 . A colleague used 1.7. and it looks like it should (from blue to Red)
when I use QGis 1.8 the complete raster is red with an unreasonable classified legend with four colours (three with minus values, one with maximum value red)
any ideas? Bug, hidden button, secret?

Thank you for help!!

Sorry, I am yet not allowed to attach pics for better description


texts:
I have a set of GPS coords with X,Y and Course. I need to interpolate additional points based on the course of each point. It would be a simple line intersection based on lines drawn using the 2 locations and course.

EDIT:
Let me clarify .. i have a set of gps logged points, each with lat/long, heading and speed. I need to interpolate an additional point between each known location based on the 2 known points and there headding. Kind of like dead-reconing but after the fact. It'll be a linear interpolation.. here is an image of said problem 1...


texts:
I have two datasets that contain tree canopy mortality classifications for aerial photographs. These two datasets were completed about 10 years apart, and therefore have different classification schema. Since consumers of these datasets will want to do comparisons of the datasets, I'd like to re-bin these classification data into a 'crosswalk' of sorts, where the severity of mortality is uniform for both surveys. I'm wondering how I might use ArcPy to accomplish this.

(Above) 2009 system of rating mortality

(Below) 2019 system of rating mortality

I would like to create a crosswalk for these two datasets so I can bin them using a schema common to both of them for comparison's sake.


texts: Atmospheric correction for RapideEye without atmospheric parameters

texts:
I'd like to perform an atmospheric correction of my RapidEye images but metadata do no include atmospheric parameters required to use the 6s algo implemented in the Orfeo ToolBox. Is there a solution to still apply the corrections?

Browsing on the web, I found that in grass it was possible not to specify them (http://grass.osgeo.org/grass64/manuals/i.atcorr.html) but rather define if the site is midlatitude summer, midlatitude winter, subarctic summer, subarctic winter, etc. Do those areas have predefined parameters and if yes where can I find them?


texts: Saga GIS at command line in OSGeo4W - 'shapes to grid' tool - sqlite input

texts: I have been using the Saga 'Rasterize' plugin in QGIS to rasterize spatial data from sqlite format (road networks, and intersections) to tiff. I find this tool more flexible than the gdal_rasterize equivalent. Whilst i can batch convert files in qgis, i would like to implement this tool (actually saga 'shapes to grid' tool - which i believe is one and the same as the qgis saga rasterize tool) at command line in windows batch scripts. However, I am struggling to find the correct syntax to do this - in terms of selecting the sqlite file and then the table within it. How do i enter the file and table into the saga tool at command line? There only seems to be the facility to enter an input file and a table attribute at command line, no input table name? The tool must surely support sqlite, as the equivalent in qgis works fine?


texts:
I got the GeoExt QueryPanel working so far, it recognizes the layers available on my Geoserver but it is not able to get the attribute columns for the dataset. I configured an osgeo proxy on my Apache Server, I use the override-ext-ajax.js and have the OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url="; class configured.

A DescribeFeatureType request is sent when selecting the layer I want to filter. The same happens when trying to choose the attribute column in the query panel, but there are none displayed.

The sent request URLs look like this:

http://localhost/cgi-bin/proxy.cgi?url=
http%3A%2F%2Flocalhost%2Fcgibin%2Fproxy.cgi%2Fwfs%3F
VERSION%3D1.1.0%26REQUEST%3DDescribeFeatureType%26TYPENAME%3Dcases%26query%3D

The strange part here is, that the request is sent via the proxy.cgi twice. The request URL should look like this, going passed geoserver´s WFS:

http://localhost/cgi-bin/proxy.cgi?url=
http%3A%2F%2Flocalhost%3A8080%2Fgeoserver%2Fwfs%3F
VERSION%3D1.1.0%26REQUEST%3DDescribeFeatureType%26TYPENAME%3Dcases%26query%3D

Does anybody have a clue what´s going wrong here and why the request URL is wrong?

The strange thing is that even if I remove the Open.Layers.ProxyHost class and the override-ext-ajax.js from my application, the request is sent that way.

Can it be related to my server?

Other requests sent in the application (GET & POST) work fine. The (relevant) code looks like this, maybe you can spot the error.

OpenLayers.ProxyHost = "/cgi-bin/proxy.cgi?url=";

    var owsUrl = "/geoserver/wfs";

    var symbolizer = OpenLayers.Feature.Vector.style["default"];

    symbolizer.graphicName = "square";

    var store;

    Ext.onReady(function() {

                store = new GeoExt.data.WFSCapabilitiesStore({
                url: Ext.urlAppend(owsUrl,
                   "VERSION=1.0.0&REQUEST=GetCapabilities"),
                protocolOptions: {version: "1.1.0"},
                autoLoad: true,
                listeners: {
                    load: initApp
                }
                });
            });

    function initApp() {
            var vecLayer = new OpenLayers.Layer.Vector("Query", {
                eventListeners: {
                    featuresadded: function() {
                        var extent = vecLayer.getDataExtent()
                        mapPanel.map.zoomToExtent(extent)
                    }
                }
            });

var queryPanel = new gxp.QueryPanel({
                title: "Query",
                region: "center",
                width: 350,
                bodyStyle: "padding: 10px",
                autoScroll: true,
                layerStore: store,
                bbar: ["->", {
                    text: "Query",
                    handler: function() {
                        queryPanel.query();
                    }
                }],
                listeners: {
                    query: function(panel, store) {
                        var rule = new OpenLayers.Rule({
                            filter: panel.getFilter()
                        });
                        rule.symbolizer[queryPanel.symbolizerType] = symbolizer;
                        var style = new OpenLayers.Style("", {rules: [rule]});
                        wms.mergeNewParams({
                            sld_body: new OpenLayers.Format.SLD().write({
                                namedLayers: [{
                                    name: panel.selectedLayer.get("name"),
                                    userStyles: [style]
                                }]
                            })
                        });

                        store.bind(vecLayer, {
                            initDir: GeoExt.data.FeatureStore.STORE_TO_LAYER
                        });
                    }
                }
            });


texts: Intersects query with where clause and max time

texts:
I am having an issue with this query. I am trying to select the buildings (gta_bldgs) that intersect with gps points (wherez) but when the speed is equal to zero and the most recent (Max devicetime). I have been struggling to make it work.

SELECT pkuid, name, type, a.geometry
FROM gta_bldgs AS a, (SELECT deviceid, devicetime, speed, geometry FROM wherez WHERE speed = 0) AS b
WHERE ST_Intersects(a.geometry,b.geometry)
GROUP BY b.deviceid
HAVING b.Max(devicetime);


texts: Currently i am serving map tiles with gunicorn and tilestache with nginx reverse proxying gunicorn. This works, although for some reason i am unable to do anything with the appearance of the map as it keeps defaulting to this modest maps index. Currently i am using mbtiles for this setup however i would like to scrap gunicorn and tilestache and serve the .png files directly from nginx, but i currently have no idea how to go about doing that. I seek help in the configuration for this kind of setup. I have used mbutil to convert my mbtiles into a directory for the .png files.


texts: Export Web Map is not exporting layout template

texts:
I'm trying to share the Export Web Map Task as a service in arcgis server so that I can integrate the printing functionality in my CMV(Configurable Map Viewer).
To do so, in Arcmap I'm passing a layout template I set up. The Export from Arcgis desktop does not render any error and also while publishing to arcgis server it tells me that the layout is going to be exported but when I consume the service URL, the Layout is not added in the choice list in the CMV. I repeated the export task many times yet the issue is presisting.
If you have any idea about the plausible causes of this behaviour it would be great to let me know.


texts: Query an MS Access geodatabase for specific field shape data in ArcMap

texts:
I'm trying to help my wife with something in ArcMap.  She works with ArcMap/ArcGIS every day but has no programming or database knowledge beyond creating basic tables, queries and forms in Access.  I, on the other hand, have good knowledge of programming and databases, but none with ArcMap (so I apologize in advance if my terminology is misleading or flat out wrong).

We have an Access database which is basically a manure spreading log and we're trying to create a report that will pull a particular days spreading activity (ie, which fields were spread, acres, how much was spread, etc.).  Along with that, she'd like to generate/pull the shapefile data for only those fields so that when working in ArcMap she doesn't have to pull in the shape data for every single field (there are a LOT).

We've successfully exported the shapefile data into an Access database, but we're stuck at being able to query it and get the query to work within ArcMap.  What I mean by that is, the data we exported shows up in the catalog in ArcMap and the "table" can be dragged over to the map layers to draw the shape data for all the fields.  We're trying to create a query we can drag over to draw only the fields she's currently working with.  Is there any way to do this?


texts: Assign nearest raster cell value to point

texts:
I came across to a similar issue described in Finding nearest raster cell value based on vector point?

I understand the proposed solution with Raster to Point and then Spatial Join with (for example) Closest Geodesic option. However, is not clear how to use the coordinate formula described by @whuber:

(x',y') = ( min(max(x,x0),x1), min(max(y,y0),y1 )

For example on my global raster x0=-180, x1=180, y0=-90 and y1=90
If I apply the above formula for each given (x,y) then I get no difference (no shift). That is obvious, because all the (x,y) points for which I want to extract value from raster are within the extent of the raster, just that they fall in NoData cells. I presume the formula makes sense when knowing the x0,x1,y0,y1 of the nearest cell that contains any valid value different from NoData (which I need to find in the first place).


texts: Can't load pgrouting extension, could not load library librouting_dd.ddl error 998 windows

texts: I have a windows 7 machine 64-bit, I installed Postgresql 9.3 via the opengeo suite 4.0.2 in the program files (x86) folder, so I guess I have to install the 32-bit version of the extension. I am trying to load the 2.0.0 version of pg routing and I always get the mentioned error message. Anyone has suggestions in order to get work?


texts:
I am a newbie working on QGIS 3.6 / Windows 10

I need to style a legend for map with headings and sub headings
The main divisions of the legend item need to be subdivided into sub-classes

i.e. The topmost items each represent a group of biotope types represented as the categorised colours of the polygon features in the map.
- see herefor details: https://www.stadtentwicklung.berlin.de/umwelt/umweltatlas/ed508_04.htm

the "sub-classes" of each of these group of biotopes have individual biotope-codes that I need to label as shown in the other question I posted here:
rule based coloured background for labeling in QGIS

I have attached a pic of the legend as it should look below:


texts: Mapping matrix of origins and destinations in QGIS?

texts:
I have a central point on a map, and I have other points in a separate layer (the location of a venue, and the origin cities of the venue visitors). Is there a plugin which would generate output similar to airline destination maps like this;

They don't have to be "great circle" type of lines, just straight lines will do. This is for visual representation of point of origin for visitors of certain events.

Any tips or plugin recommendations?


texts:
I'm taking some photos of my city and at the same time I'm recording the track so that I can geotag the images when done.

The problem was that couldn't do it with all the software offered (HoudahGeo, Geotag...). I started wondering why this couldn't work and I finally discovered that when I opened my gpx, the track points were one hour less than the creation date of the images. I don't know why? The timing should be the same because Im using the same device.

If you have any idea about how to deal with this i would be pleased.


texts:
I'm trying to retrieve the informations about points around the mouse click. However if the click is not on the inside the symbol but just outside, the WMS returns an empty response.

To simplify the question here are two screenshots:

Here the response of the server contains the information

Here the response of the server contains nothing because the point is in the middle.

I would like to know if there is a way to ask the GeoServer WMS to return the features with more tolerance (given that the parameter feature_count is already set at 50)


texts: Error returned when harvesting from GeoNetwork instance – bad-xml-response

texts:
I’m attempting to harvest metadata records from a target GeoNetwork (v3.10.0.0) instance and am receiving a “Harvest error bad-xml-response” with the following returned response:

Please see full response file here.
Target Env:   Deployed from .war file running on Tomcat 9 and CentOS 7.
The CSW test calls (csw-GetCapabilities, csw-DescribeRecord, csw-GetRecords) from within the target instance UI seem to work ok, but not sure if they actually use the external URL Get interface or if they just exercise the GeoNetwork library calls internally.
We have other harvesters set up and are harvesting records from similar(but not same) instances just fine.  There have been some minor customizations done to the landing page of the target instance (noted in yellow highlight in the full response file), but I don’t understand how/why that would affect the harvest request.


texts:
1)  Is there any web API which can show 3D trajectory below earth surface and can consume the ArcGIS feature service which can work closely with ArcGIS web platform.

2)  Is there any way to publish 3D feature service in ArcGIS Server 10.1 and can be consumed in web application and can be displayed on map.


texts: Interpolate and export the following points

texts:
I have two crossections of the same position.
I am looking for a way to interpolate the blue points from the attached image.

I need the xyz values of the interpolated points, how can I do this?

I use QGIS. But the points are exported from Terrascan. I tried to come up with a calculation in Excel.


texts: Google Maps mixed results for ATM and Banks

texts:
When I search for nearby ATM at a particular latitude and longitude I get top 20 results but, the results are mix of ATM and the Banks. Why do I get a mixed result and, as I am not interested in bank results, how can I modify my search?

Below is the nearby search result.


texts:
How can I get the symbol from a FeatureLayer?
I want to use the symbol of a FeatureLayer to symbolize a graphic layer.


texts:
I've got some linestrings covering a road network, and need to split them at the crossroads.
Does anybody have any idea on how to do it?

I'm using QGIS 2.4, with a Postgis 2.1.3.


texts: How to undo misclicks while drawing a polygon?

texts:
I'm trying to draw some polygons in QGIS 2.2 and I lose it, because everytime I misclick, I don't find a way to undo that point and start drawing again.

I'm sure there's a way. Does anyone know it?


texts: GDAL change spatial resolution via function

texts: i have an geoTIFF with 30 arc Seconds spatial resolution. I would like to change its spatial resolution to 60 arc Seconds. Can i use an gdal function to change the 30 arc Seconds resolution to 60 arc seconds?