texts: Compute 3D seismic fold with QGIS and Python

texts:
Computing 3D seismic fold is quite straight forward, for those unfamiliar with the concept, below is a brief explanation.
Background information (Seismic 101) :
Seismic Fold:
Basically, I have (several hundred thousands of) Receiver Points (RP) and Source Points (SP). Each Source, when activated, will send a signal that will be picked up by (several thousands of) Receivers.
For every pair (SP, RP) I will compute the Midpoint, i.e. the coordinates of the middle of the segment (the ray) joining the 2 points and find out in which bin this Midpoint falls.
The bins grid is a regular grid (usually orthogonal/rectangular) with cell dimensions being RP_i/2 and SP_i/2 where _i denotes the interval between points (these intervals are fixed as they are taken on the theoretical RP/SP grid).
The Fold Map will then simply be a heatmap with the values being the number of Midpoints falling into each bin of the bins grid.
Ideally, I'll need to be able to filter the Midpoints based on the Offset length. The offset being the length of the ray (i.e. the distance between SP and RP for the pair considered).
A typical fold map may look like this.
SPS (Shell Processing Support) files:
The information source is usually in the form of 3 SPS files (see format definition here). I am able to load the 2 Point files (R and S) into QGIS so I get all Receivers and Sources positions.
The Relation (X) file gives the list of all RP listening to each SP.
The problem:
The problem lies mainly in the numbers. I have to compute several million rays so every nano-second will count here!
The X file lists the points by name (or rather, number) so I'll have to lookup the coordinates of these points every single time as I'm reading through the X file.
Is there an efficient way to load my R and S files in memory (ideally caching only a portion of these files at a time) and how should I hold the information of the computed Midpoints ? In memory ? In a numpy array ? Is there a QGIS object that would be particularly suitable for this kind of extensive computation ?


texts: I'm looking to use the geojson module in Python, as well as Pilllow to convert a .geojson to a .png image. I know how to extract the data from the file but I don't know what to do from there.


texts: ogr2ogr MSSQLSpatial truncating attribute fields

texts:
I have a Shapefile with a number of attribute fields defined in its corresponding .dbf file. I am trying to transfer this data to an Azure SQL Server DB using the following command:

ogr2ogr -f MSSQLSpatial "MSSQL:server=<server>;database=<db>;uid=<username>;pwd=<password>;Driver={ODBC Driver 17 for SQL Server}" <shp file>

This command adds all the features to the database. However, after closer inspection I have noticed that all non-numerical fields have been truncated to a single character.

Now, the strange thing is that the types of the SQL columns were instantiated such that they should be able to fit the data. Could this be an encoding problem, as some of the fields contain French characters?

I am using GDAL 2.2.3, released 2017/11/20


texts: Map extent based on atlas feature attribute

texts:
I am using Atlas on QGIS 2.18.17.
The city I am working on is divided by wards and these wards are divided by blocks.

Each page of the atlas is focusing on one block (layer coverage: block.shp).

There is a small map on each page displaying the ward where is located the block.
To show the ward where is located the current atlas feature (i.e. the current block) the extents are set as follow:

A:  x_min( geometry(get_feature('ward layer','WardName',attribute(@atlas_feature,'WardName'))))

B:  y_min( geometry(get_feature('ward layer','WardName',attribute(@atlas_feature,'WardName'))))

C:  x_max( geometry(get_feature('ward layer','WardName',attribute(@atlas_feature,'WardName'))))

D:  y_max( geometry(get_feature('ward layer','WardName',attribute(@atlas_feature,'WardName'))))

To improve a bit the visibility, I will add a 15% margin around feature as follow:

x_min = A - (C - A) *15%

y_min = B - (D - B) * 15%

etc.

Would it be interesting to improve the option "controlled by atlas" by adding an option such as "set extent to another layer than atlas coverage layer"? (it's a kind of request for feature improvement in fact but I don't know where to post this)


texts: I am new to QGIS and am trying to convert a .dxf file to a .dgn file.  This is my first time working with .dxf files so I'm not even sure what the difference is.  Working with Version 3.8.1.


texts: Geoserver VM with NAS, Should I store entire VHD on the NAS or just map the data directory? Or just go with internal SSDs w/ weekly backups?

texts:
Our humble geoserver instance needs to be able to host GEOTIFF image pyramids ranging in size from a few gigs to a few TBs. A NAS appliance with raid 10 seems smart from a maintenance/expandability perspective but wouldn't a few internal SSDs beat the snot out of it performance wise?

related: Disk access characteristics of accessing vs building an image pyramid.

host: 12 core xeon, 16GB ram, windows server 2008 running hyper-v and an ftp server. It has two gigabit ethernet ports and a few USB 3.0 ports.  One ethernet port is reserved for the uplink. The other is not being used.

vm: Ubuntu 14.04 running Geoserver 2.7.2 with 8GB ram and 4 dedicated processors.

uplink: 320/320 Mbps or 40/40 MBps

edit:

We need a NAS regardless for imagery backups.   I am leaning toward putting a SSD in the host machine for the VM's VHD, Then I will mount the geoserver data directory on the NAS and perform periodic backups of the VHD to the NAS


texts: shp2pgsql creates one dimensional shapes

texts:
I use shp2pgsql to import an .shp ESRI shapefile into PostgreSQL database. I use brew installed shp2pgsql. I used this command first:

shp2pgsql -I -s 4326 /path/to/shape_file.shp public.os_gb | psql -U postgres -d postgres

It shows this message on terminal:

It imports data into database but PostgreSQL Geometry Viewer says that it is not able to render "3D" shapes.

Then I tried this command:

shp2pgsql -I -t 2D -s 4326 /path/to/shape_file.shp public.os_gb | psql -U postgres -d postgres

Geometry Viewer shows 1 geom like this:

So forcing in 2D did not work. Could you please help me with this and tell me what's wrong here.

Edit:

Here is the result of "ogrinfo -so -al my_shape_file.shp":

INFO: Open of `GB_GreenspaceSite.shp'
using driver `ESRI Shapefile' successful.

Layer name: GB_GreenspaceSite
Metadata:
  DBF_DATE_LAST_UPDATE=2018-09-24
Geometry: 3D Polygon
Feature Count: 143822
Extent: (9819.840000, 8274.570000) - (655230.650000, 1214133.190000)
Layer SRS WKT:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","27700"]]
id: String (38.0)
function: String (40.0)
distName1: String (254.0)
distName2: String (254.0)
distName3: String (254.0)
distName4: String (254.0)


texts: ArcMap "Polyline to Raster" output is empty greyscale raster

texts:
I have a polyline shapefile of streams that I am attempting to convert to a raster using the "Polyline to Raster" conversion tool. However, rather than outputting a raster with distinct stream pixels, my output is a seemingly empty raster with a greyscale value scale of 0 to 255. My methodology is as follows:

1) I have created a field within "streams" called "rastervalue"
2) I assigned values of "1" to each stream feature in hope that all stream
pixels within my output raster are represented by this value.
3) I ran "Polyline to Raster" using the following script:

    arcpy.PolylineToRaster_conversion(Streams,"rastervalue")

I have attempted this with varying cell sizes, and my results seem to be the same. How can I generate the desired output?


texts: Which format to display the gsmap downloaded data in QGIS?

texts: Where can I download gsmap data especially for precipitation / rainfall and which format can be supported by QGIS software?


texts: Convert X and Y values from an Israeli source to match the WGS 84 / Pseudo-Mercator CRS in QGIS

texts:
I have an Israeli database with  X and Y values that I need to upload into QGIS.

I have converted the Excel file into a CSV UTF-8 (Comma delimited) file and upload it as Delimited Text Layer.

I have tried all the Israeli CRS listed in QGIS, and that data is still not matching.

It sure do look like the Israel outlines, but it is about 700kM off somewhere in the Egyptian desert.

What is needed in order to fix this problem?


texts: No data values collars error - delete values from image

texts:
I have a satellite image that has collars outside of the image. When I set no data value as 0 and set this colors as no display it also removes the pixels from the image with 0 - RGB.
The image is jpg2 format and is an satellite image (Pleiades).
I have also used copy raster tool and the result is the same (the dark outbound area is disappear but also the pixels from the image with the 0 value).
How can I deal with this  issue?


texts: Problem with Sentinel 1 images in Google Earth Engine

texts:
I am trying to apply Otsu algorthm to S1 images starting from its launch to till date on my study area for waterbody mapping. I perform otsu without applying any date filter for Satellite pass so I can get more imageries. But I get the error

ImageCollection (Error) Error in
map(ID=S1A_IW_GRDH_1SSV_20160525T235440_20160525T235507_011422_011607_60C4): Dictionary.get: Dictionary does not contain key: bucketMeans.

Then I tried to omit the above error by applying
.filter(ee.Filter.neq('system:id', 'COPERNICUS/S1_GRD/S1A_IW_GRDH_1SSV_20160525T235440_20160525T235507_011422_011607_60C4'))

and again I get error

ImageCollection (Error) Error in
map(ID=S1A_IW_GRDH_1SSV_20160805T235445_20160805T235511_012472_0137E8_08F2): Dictionary.get: Dictionary does not contain key: bucketMeans.

and I apply another filter
.filter(ee.Filter.neq('system:id','COPERNICUS/S1_GRD/S1A_IW_GRDH_1SSV_20160805T235445_20160805T235511_012472_0137E8_08F2'))

But nothing works it just goes on. I came to understand the issue happens only on Descending passes. So when I put a filter to obtain Ascending pass the output progress without any errors but it won't apply the same in Descending pass.
Here is my code..


texts:
The application that I work with consumes WMS services created in Geoserver via MapProxy. I used ysld to style my layers. As long as I want just display my layers in the application, the scales that I configured via ysld are preserved. When I print the same WMS service scales change.
I found in the MapProxy documentation (Scale vs. resolution: https://mapproxy.org/docs/1.9.0/configuration.html?highlight=scale) that the scales are changing depending on dpi. The one that I use is 300. Is there any way to keep MapProxy from changing the scale?
Here is an example of my ysld style that looks different if it's displayed or printed.

feature-styles:

- name: Ortsname
  rules:
  - name: Ortsname (1:2000 - 1:9500)
    scale: [2000,9499]
    symbolizers:
    - text:
        label: ${name}
        fill-color: '#000000'
        font-family: Arial
        font-size: 12
        font-style: normal
        font-weight: normal
        placement: point


texts: Converting GOES x,y coordinates from a 2D fixed grid into latitude and longitude

texts:
I am working with some GOES data in Python and I want to grab latitude/longitude coordinates from the data. GOES' stores their projection data in radians under x,y variables that are used to project onto their fixed grid. Their PUG has a formula to convert between the two, however when I try to implement this I run in to an issue where a lot of data is lost due to a math domain error.

Here is the formula from the PUG and my code:

proj = nc.variables['goes_imager_projection']

major,minor = float(proj.semi_major_axis), float(proj.semi_minor_axis)

H = float(proj.perspective_point_height)+major

lambda_0 = float(proj.longitude_of_projection_origin)

a = (math.sin(x))**2 + ((math.cos(x))**2)*(((math.cos(y))**2) +
((major**2)/minor**2)*((math.sin(y))**2))

b = -2*H*(math.cos(x))*math.cos(y)
c = H**2 - major**2

r = (-b - math.sqrt(b**2 - 4*a*c))/2*a

s_x = r*math.cos(x)*math.cos(y)
s_y = -r*math.sin(x)
s_z = r*math.cos(x)*math.sin(y)

lat = math.atan(((major**2)/minor**2)*(s_z/math.sqrt((H-s_x)**2 + s_y**2)))
lon = lambda_0 - math.atan(s_y/(H - s_x))

However, I keep getting a math domain error when it tries to compute r.

I know why it's tripping that math domain error, what I want to know is how do other programs bypass this? I know there is software that can convert from 2D grids to latitude and longitude (i.e Basemap), so what is it that my code is missing to do this?

If you wish to retry this, the values I am getting from the goes_imager_projection variable are as follows:

perspective_point_height: 35786023.0
semi_major_axis: 6378137.0
semi_minor_axis: 6356752.31414
inverse_flattening: 298.2572221
latitude_of_projection_origin: 0.0
longitude_of_projection_origin: -75.0
x = -0.106483996
y =  0.135044


texts: Summarizing (Statistics_analysis) with field list in ArcPy?

texts:
I want to summarize a field in 5th order in the shapefile. For that reason, I write a script below

output_lines = r"C:\Outputs\lines.shp"
sum_table = r"C:\Outputs\summary.dbf"

Listfield = arcpy.ListFields(output_lines)

for field in Listfield:
        arcpy.Statistics_analysis(output_lines, sum_table, [["Total_Leng", "SUM"]], field[5])

I get the error which is "TypeError: 'Field' object does not support indexing". I'm using ArcMap 10.3


texts:
I have used the software Ecognition to create a land cover & land change map. I used the class Barren in two images. However, I do not want to use the class Barren anymore. I have reclassified the image objects classified as Barren to a separate class. I cannot remove the class Barren from the Class Hierarchy window.

When I chose delete or delete class I get the following error message:

cannot delete class because it is in use

My only other option is to delete all the classes which I cannot do because I did some manual classification. Please help me figure out how to delete one class without deleting the entire classification in the image.


texts: arcpy.da.SearchCursor keeps running, without any response or error

texts:
I am using 'arcpy.da.SearchCursor' in Python script for ArcGIS 10.5, while the script is always keeping running, without any result and error. I use print to see where the problem is; it seems that 'arcpy.da.SearchCursor' cannot be finished, as it does not print 2. However, this works well if I use these syntaxes in the Python window of ArcMap.
Here is the code:
import arcpy

fc = r'C:\\Users\hsdfyan1\Desktop\add_point_attributes\Feb_8_1_e_pth_split_join.shp'
fileds = ['time_first', 'time_last', 'SHAPE@']
print(1)
cursor = arcpy.da.SearchCursor(fc, fileds)
print(2)
for row in cursor:
    print(row[0],row[1])

The only result is 1


texts: Deleting polygon defined area from multi band raster layer in ArcGIS Desktop?

texts:
This is a similar question to Deleting polygon defined area from raster layer in ArcGIS Desktop?, but I need the answer to be able to apply to multi bands since I am working with a RGB file.  I have a Landsat 8 file and have used SAGA to classify the clouds (wasn't to fond with ArcMaps).  I  I want to use arcmap to remove the polygon clouds from my original landsat image (SAGA is too slow and crashes constantly).  From the post linked above, I have converted the shape file to a raster.

I was thinking may changing all the cloud raster values to zero (not sure how to do this) and then multiplying with the original landsat 8 raster.

Any suggestions?


texts: How to render Multi-Response GeoJSON tiles

texts:
I am using TileStache to combine and server multiple vector layers in a single Mutli-Response GeoJSON tile, which is basically all my individual tiles/layers (GeoJSONs containing either roads, pois, buildings) merged into one.

When TileStache combines the layers, the original layer name (like "roads" or "pois") is added before the features of that layer are included in the GeoJSON, so when I parse the tile to render it in leaflet (using the TileLayer.GeoJSON plugin) it doesn't recognize the objects in the GeoJSON file. It just skips over them and doesn't render anything.

This is what the sample response looks like:

{
    "roads": {
        "type": "FeatureCollection",
        "features": [{
            "geometry": {
                "type": "LineString",
                "coordinates": [
                    [73.12500, 33.61371],
                    [73.12432, 33.61462]
                ]
            },
            "type": "Feature",
            "id": "ff7274b342",
            "clipped": true,
            "properties": {
                "priority": 5.00000,
                "name": "Islamabad Highway"
            }
        }]
    },
    "boundary": {
        "type": "FeatureCollection",
        "features": []
    },
    "pois": {
        "type": "FeatureCollection",
        "features": [{
            "geometry": {
                "type": "Point",
                "coordinates": [73.10332, 33.60876]
            },
            "type": "Feature",
            "id": "0d39282b10",
            "clipped": true,
            "properties": {
                "priority": 1,
                "name": "BENAZIR BHUTTO INTERNATIONAL AIRPORT ISLAMABAD"
            }
        }]
    },
    "cartographic": {
        "type": "FeatureCollection",
        "features": [{
            "geometry": {
                "type": "MultiPolygon",
                "coordinates":
                    [

                        [
                            [73.10095, 33.61462],
                            [73.10065, 33.61438],
                            [73.09731, 33.61462],
                            [73.10095, 33.61462]
                        ]
                    ]
            },
            "type": "Feature",
            "id": "963cec5810",
            "clipped": true,
            "properties": {
                "category": 6,
                "priority": 5,
                "kind": "port places",
                "name": "BENAZIR BHUTTO INTERNATIONAL AIRPORT"
            }
        }
}

And this is the code I am using in Leaflet:

var geojsonURL1 = 'http://localhost:8080/composite/{z}/{x}/{y}.json';
var geojsonTileLayer1 = new L.TileLayer.GeoJSON(geojsonURL1, { maxZoom:22,
        //clipTiles: true,
        unique: function (feature) {
            return feature.id;
        }
    }, {
        style: style,
        onEachFeature: function (feature, layer) {
            if (feature.properties) {
                var popupString = '
';
                for (var k in feature.properties) {
                    var v = feature.properties[k];
                    popupString += k + ': ' + v + '
';
                }
                popupString += '
';
                layer.bindPopup(popupString);
            }
            if (!(layer instanceof L.Point)) {
                layer.on('mouseover', function () {
                    layer.setStyle(hoverStyle);
                });
                layer.on('mouseout', function () {
                    layer.setStyle(style);
                });
            }
        }
    }
);

map.addLayer(geojsonTileLayer);

How can I render them in Leaflet? Are there any viable alternatives to render the GeoJSON tiles?


texts: Opening Sentinel 2 images with GDAL in Python

texts:
I'm trying to open Sentinel-2 images with python, and as far as I can tell, I have GDAL 2.1 installed, which (so I understand) ships with a Sentinel 2 driver. At least, the lines

from osgeo import gdal
for i in range(gdal.GetDriverCount()):
    drv = gdal.GetDriver(i)
    if drv.GetMetadataItem(gdal.DCAP_RASTER):
        print(drv.GetMetadataItem(gdal.DMD_LONGNAME), drv.GetMetadataItem(gdal.DMD_EXTENSIONS))

include "Sentinel 2" on the list of outputs.

However,

gdal.Open(filename, GA_ReadOnly)

returns "None" for filenames that point at the zip-files from scihub, the xml-file or the jp2-files. Similarly, the gdalinfo tool returns that it is not able to recognise the file format.

Any pointers on where I am doing something wrong?


texts:
My measure dialogue box seems to have taken leave of my QGIS. I've reviewed the previous answers about resetting the interface, checking for the dialogue box outside of view, all kinds of things.

The measure line does appear, but no measurement appears. Any suggestions on how to restore it?


texts: How to write log messages from arcpy to ArcGis server?

texts:
Published geoprocessing services write logs that are available at ArcGis Server Manager (arcgis-server/arcgis/manager/log.html)

Arcpy has AddMessage, AddError, etc methods for logging, but they are doesn't works when script is published as geoprocessing service. As I understand, they works only locally in ArcMap.

Is there some way to write server logs from arcpy?


texts: Openlayers 5 VectorSource getextent only of certain features inside source

texts:
i have a VectorSource built from GeoJSON, and i know how to get the extent of all the features inside the vectorSource, but i don't know how, or what the best method is, to get the extent only of certain features.  Basically instead of

const featExtent = featEntrySource.getExtent();

I'm looking for something like (in pseudo code)

const featExtent = featEntrySource.getExtent(<where each feature.get('referenceid') === certainValue>);

whats the best way to do this?


texts: QGIS2WEB outputs data layers as JScript "Script File" ... can I convert them back to .shp?

texts: I have been given a set of folders with a .html file,  that was exported from QGIS plug-in "qgis2web".  Normally you can just click on the .html file and it will open a  local version of a QGIS authored website.   I don't know the origin of these folders and I don't have access to the original data.  The spatial layer is stored in a 'data' folder is a JScript Script File (.js).   Is it possible to pluck these files out and view them in another program or ideally convert them to shape files (.shp)?


texts: Conditionally Map.addLayer() in Google Earth Engine

texts:
My app has a dateSlider that searches multiple collections, and selects the image nearest to a date within a window. I want to conditionally display those images, with error handling for empty images.
var image = l8sr
    .filterDate('2016-01-01', '2016-07-01')
    .first()

var nullImage = s2sr
  .filterDate('2016-01-01', '2016-07-01') // returns empty collection
  .first()
print(nullImage) // null

Map.addLayer(nullImage) // Error! Expression does not evaluate to an image.

// Check before trying to add the image to map
if(nullImage !== null){Map.addLayer(nullImage)} // Same error

// Dosen't work server-side either
ee.Algorithms.If(nullImage !== null, Map.addLayer(nullImage)) // Same error

I want my script to check each image before it tries to add a layer, to avoid the error message. In other words, I want my script to be robust to handle empty image collections, so it can flexibly add layers to the map. How can I accomplish this?
https://code.earthengine.google.com/0502525f29a25d7f1d4ab7b63acde203


texts: MODIS HDF file without Coordinate System?

texts:
I am using the MOD04_L2 MODIS Aerosol Deep Blue product from the LAADS API. These files are in hdf4 format, but I am transforming them to hdf5. I did transform the files because I'm unable to compile gdal with hdf4 support, and I found this great command-line tool from the HDF Group that allows to easily translate files between versions.

All is good until there. But when I ran a gdalinfo I get this:

Files: Downloads/MOD04_L2.A2000055.0005.061.2017205182825.h5
Size is 512, 512
Coordinate System is `'
Metadata:
  ArchiveMetadata.0=
GROUP                  = ARCHIVEDMETADATA
  GROUPTYPE            = MASTERGROUP
(...)

There is no CRS or coordinate information. The only relevant information is the Corner Coordinates:

Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

But there is not much I can do with that.

It is my first time working with the hdf format.

Is this missing data due to the modification of the file to hdf5?

The command-line utility is already compiled without documentation.


texts: Find distance between two points by following a polyline in R

texts:
I have a polyline shapefile that follows the center of a stream, and I have coordinates for animal positions within the stream (not necessarily immediately on the line), and I am trying to calculate the distance between points by going from point1 to the line, following the line to point2, then going from the line to point2 (see illustration below). This should return a single value that (given the width of my stream) is a fairly close approximation to the distance and animal would have to swim to get between two points. I need to be able to do this in R. (note that points are often near the middle of the line, so I need to find the path to go to the nearest point on the line).
I have spent several hours trying to find a solution, and while I have found some similar questions, I have yet to find an answer that actually does what I need.
Thanks
In the image below, the blue dotted line is the total path for which I want a distance


texts: Writing GDAL/OGR point shapefiles over 4GB size?

texts:
I am writing point shape files for a Monte-Carlo predictive model. My shapefiles are exceeding 4 GB in size. I see corrupt attribute table entries.

Is this a result of the large file size?

I'm using ArcGIS Pro to examine the attribute tables. I've been searching my C++ code for other memory overwrites (clobbering) and other possible problems.

Could it just be a matter of trying to write too big a file?


texts: Select the features within a specified number of attribute groups closest to a point

texts:
I am having a difficult time wording this question properly so I apologize for the vague title. I have a table of subway stops in NYC and what I am trying to do is select the points from the five closest subway stops to a given point. However, each subway 'stop' has multiple points associated with it due to there being many entrances for each stop. The table looks like this:

ogc_fid
line
name
objectid
url
wkb_geometry

1
2-5
Birchall Ave & Sagamore St at NW corner
1734.0
http://web.mta.info/nyct/service/
POINT (-73.86835600032798

2
2-5
Birchall Ave & Sagamore St at NE corner
1735.0
http://web.mta.info/nyct/service/
POINT (-73.86821300022677

3
2-5
Morris Park Ave & 180th St at NW corner
1736.0
http://web.mta.info/nyct/service/
POINT (-73.87349900050798

4
2-5
Morris Park Ave & 180th St at NW corner
1737.0
http://web.mta.info/nyct/service/
POINT (-73.8728919997833

5
2-5
Boston Rd & 178th St at SW corner
1738.0
http://web.mta.info/nyct/service/
POINT (-73.87962300013866

6
2-5
Boston Rd & E Tremont Ave at NW corner
1739.0
http://web.mta.info/nyct/service/
POINT (-73.88000500027815

If I split the 'name' field by "at" and select the first slice I can group all rows at the same 'stop', however, how do I then return all rows associated with the 10 closest groups (stops)?
For instance, when I use this query:
 select wkb_geometry
 from osm.subwaypoints s
 where s.line like '%C%'
 order by ST_Distance(ST_SetSRID( ST_Point(  -73.96833020663865,40.68396650310555), 4326),s.wkb_geometry ) LIMIT 10

I get the nearest 10 entrances, but it only accounts for 2 actual stops


texts:
I'm working with data points (x,y) and timestamp t as a separate attribute in my table. (x, y) are in WGS 84 latitude,longitude format while timestamp is in Unix Time. For instance (71.8, -13.5) and 1420246483

I want to find if 2 points are close to each other both in terms of physical distance (meters) but also within a certain time. This is very simple using a query like

SELECT id, geom
FROM points p1, points p2
WHERE ST_DWITHIN(p1.geom::geography, p2.geom::geography, 1000)
AND (p2.timestamp > (p1.timestamp - 600) OR p2.timestamp < (p1.timestamp + 600))

in which we transform our geometry to geography so that we can say 1000 is the distance in meters between the 2 points, and with attribute timestamp we could restrict to points that are max 10 minutes (600 seconds) away.
The performance is not too bad but having millions of points it's not enough and I think it is possible to take advantage of PostGIS spatial indexes to achieve faster searches.

What I want to do now is to integrate time within the geometry points as (X,Y,time) to do a similar query which can improve performance.
I read in the Boundless webpage about 3D/4D Index that with 3D index this is possible and particularly with time.

I quote: "Some of the fancy new 3D functions like ST_3DDistance can make use of the 3D index for fast searches in high dimensional spaces. I imagine folks working with XYT data may also find the new indexes useful".

I find hard to see how such implementation can be easy to use and useful. In fact functions such as ST_3DDistance or ST_3DDWithin only work with distances of units. This mean that meters and seconds can't go together so well, one would need to convert dimensions in some way to be able to write the query

SELECT id, geom
FROM points p1, points p2
WHERE ST_3DDWITHIN(p1.geom, p2.geom, #some_unit_number)

I would like to know how experts in this wonderful community work with XYT data and find the new indexes useful as quoted above.


texts: Count points in polygon by attribute column for making choropleth

texts:
I have a point shapefile of complaints for NYC. I also have a polygon shapefile of community district boundaries with an attribute column for the community district number.

Is it possible to count the number of points that fall into each district number boundary?

I tried a spatial join, but my output file's count column just displayed a bunch of 1 values.

My final desired outcome is to have the counts for each district as a property in the polygon which I will then use to make a choropleth normalized to each districts polygon area.

Using Arcmap 10.4


texts: Trying to convert geometry to geography keep failing in MSSQL 2012

texts:
In MSSQL I have LineString in Geometry type.
When I try to convert it to geography using :

GEOGRAPHY::STGeomFromText(geom.MakeValid().STUnion(geom.MakeValid().STStartPoint().MakeValid()).STAsText(),4326)

I still get error that

Msg 6522, Level 16, State 1, Line 19 A .NET Framework error occurred during execution of user-defined routine or aggregate "geography":System.ArgumentException: 24200: The specified input does not represent a valid geography instance. Use MakeValid to convert the instance to a valid instance. Note that MakeValid may cause the points of a spatial instance to shift slightly.System.ArgumentException: at Microsoft.SqlServer.Types.SqlGeography..ctor(GeoData g, Int32 srid) at Microsoft.SqlServer.Types.SqlGeography.GeographyFromText(OpenGisType type, SqlChars taggedText, Int32 srid)

Any idea how I can convert it to geography ?

Any idea how I can share the LineString ? Its very long (exceeding this post body length)


texts: Extracting hyperspectral pixel values to shapefile?

texts:
I have developed a workflow for a question that might have an easier answer.

Say you have a hyperspectral image and a set of polygons overlapping the image. I want to extract all pixel values that fall in the polygons, such that the final dataset has n columns for n bands (with each observation being the pixel value for that band).

In ArcMap it would be the tool 'Extract Multi Value From Point', but is there a similar tool in GDAL?

My current workflow is as follows:

Load the hyperspectral image
Save each band as a separate tif
Convert each separate tif to points
Spatially join the points to the polygons
Rinse and repeat

While this works, it feels horribly inefficient and that I must be missing some multiband raster functionality.


texts: Empty vector menu in QGIS on PC with French as default language

texts:
I installed QGIS 3.0.0 on a PC with French as default language. After changing the language of QGIS to English, the 'Vector' menu is empty and there is a French 'Vecteur' menu.
Any ideas how to fix this?
I cannot change the system language.


texts: If statement not moving to else: when false?

texts:
I am looping through 2 sets of rasters to do a extract by mask of one raster using the other, with ArcPy.

However, if the mask raster does not exist I want to skip to the next set of rasters.

Both sets of rasters are named the same.

for inras in inras_list[5:]:
    rasx=arcpy.sa.Raster(inras)
    base=(os.path.basename(inras).split(".")[0])
    print base
    trial=os.path.join(base+"_wilderness_reclass.tif") #name of second set of rasters
    print trial
    rasy=arcpy.sa.Raster(os.path.join(workspace2, trial))
    print (os.path.join(workspace2, trial))
    if os.path.isfile(os.path.join(workspace2, trial)):
            inras_extent = rasx.extent
            print inras_extent
            maskras_extent = rasy.extent
            print maskras_extent
            arcpy.env.extent = rasy
            out_extract = arcpy.sa.ExtractByMask(rasx, rasy)
            outname = os.path.join(outdir, os.path.splitext(inras)[0] + "_opporwithoutroad_withoutwilderness.tif")
            arcpy.BuildRasterAttributeTable_management(out_extract)
            out_extract.save(outname)
    else:
            print "Wilderness raster does not exist for the state. Saving original ras x as is"
            out_same= rasx
            outname = os.path.join(outdir, os.path.splitext(inras)[0] + "_opporwithoutroad_withoutwilderness.tif")
            arcpy.BuildRasterAttributeTable_management(out_same)
            out_same.save(outname)

But the loop stops at the 9th file (because yes one set of rasters does not exists) and does not move to the next set of rasters i.e. 10th file.

The error I get-

RuntimeError: ERROR 000732: Input Raster: Dataset E:\Staging\GLCarbonScience\Reforestation_Attribution\Outputs\Anciliary_variables\PAD_US_Outputs\Wildnerness\Wilderness_raster_reclass\09_wilderness_reclass.tif does not exist or is not supported

Any solutions?


texts: How to display dynamic raster (lidar) data on a Web map?

texts:
I have some high quality LiDAR data. I need to process it dynamically, and update a field (not elevation field) depending on my input and elevation at that point. Once this is done, I need to display this raster data on Web Map (Bing Map, Google Map etc.) as real time as possible.
Here is what i have been able to do till now:
I have come up with some scripts using ArcGIS Python Lib, that process my lidar clips, add a field , update that field according to the input coming in. So now i basically have all the points on raster classified into 3 categories every 5 mins. If i view this raster in arcGIS i can view nice 3 colored raster overlapped on a base map.

But now i need to render these raster clips on a Web map. My first though was to polygonize these classified clips and display the polygons on a JavaScript Map. But polygonization takes a lot of time, and thus its eliminated.

Now Other option I have is to to use GeoServer's Raster rendering capabilities to display the images. Then use this WMS i created, to overlap it on Bing Maps.
This is where i am stuck. I have no clue how to automatically get Rasters that Geoserver supports, like Geotiff, jpg200o etc from my ADF files, render them on open layers first, then try to overlap them on bing maps.

If someone knows of any good tutorial that does this kind of rendering of dynamically changing rasters on map , will be great if you could point me to it.

Also if there are some other better approaches to get this kind of system in place, Please suggest. I am open to almost any solution.

Thanks

PS: I am pretty good at web programming, web services,databases  and scripting etc. But new to working with Lidar and Geo-Rasters.


texts: I can't export in QGIS a Shapefile to PostGIS

texts:
When I try export a pipe network shape to PostGIS I get this error:

Problema inserindo feições do arquivo:
C:/PP/gvSIG_Projects/Shapes/Abastecimento/QGIS_VATROCOADUTOR.shp

Ocorreu um erro no banco de dados ao executar este SQL:
SELECT AddGeometryColumn('public','QGIS_VATROCOADUTOR','the_geom',-1,'NULL',2)
O erro foi:
ERROR:  Invalid type name "NULL(2)" - valid ones are:
    POINT, MULTIPOINT,
    LINESTRING, MULTILINESTRING,
    POLYGON, MULTIPOLYGON,
    CIRCULARSTRING, COMPOUNDCURVE, MULTICURVE,
    CURVEPOLYGON, MULTISURFACE,
    GEOMETRY, GEOMETRYCOLLECTION,
    POINTM, MULTIPOINTM,
    LINESTRINGM, MULTILINESTRINGM,
    POLYGONM, MULTIPOLYGONM,
    CIRCULARSTRINGM, COMPOUNDCURVEM, MULTICURVEM
    CURVEPOLYGONM, MULTISURFACEM, TRIANGLE, TRIANGLEM,
    POLYHEDRALSURFACE, POLYHEDRALSURFACEM, TIN, TINM
    or GEOMETRYCOLLECTIONM
CONTEXT:  SQL statement "SELECT AddGeometryColumn('',$1,$2,$3,$4,$5,$6,$7)"
PL/pgSQL function addgeometrycolumn(character varying,character varying,character varying,integer,character varying,integer,boolean) line 5 at SQL statement

I can't interpret this error can you help me?


texts:
I am a complete beginner with CartoDB.

Somebody else in our organisation has created a CartoDB map of New Zealand.  They created the map by interacting with the CartoDB web site.  They didn't upload a config file or use POST /api/v1/map/named as described in the Carto named maps documentation (sorry, insufficient reputation to link to that).

And yet it is a map and it has a name - so in a plain English sense, it's a named map.  Also the config (its CSS and SQL queries) seems to be stored on the server - which ties in with what I've read about CartoDB named maps.

I need to get information from this map.  The method I use is partly determined by whether it's a named map.

The map I'm starting with looks like this when I browse our organisation's maps:

When I click on that thumbnail I'm show stuff like this:

Question

Is it a "named map"?

Bonus question (for additional gratitude or up-votes, whichever I can provide)

How can I add auth tokens for this map?


texts: Elevation API: Same coordinates, different elevation result

texts:
I'm using google maps elevation API, but while asking a single request(single coordinates) and a multi-request(many coordinates), I'm getting different elevation result of the same latitude and longitude coordinates.

Any explanation please?

The result of the single request:

The result of the multi-request with the same lat and lng coordinates:


texts: How to delete the data of one state from the TIGER geocoder

texts:
I'm getting wrong values for the coordinates using the tiger geocoder.
I have some issues with this states because some missing values and I want to upload all the data again, but now it says 'xxxxxx' relations already exist. How can I delete those data and relations for a specific state?


texts:
I used the toolbox to refactor attribute field names which worked well. However I then seem to have lost the changes as the layer was only a temporary scratch file.

How do I save to make these changes permanent for my project?


texts: Placing axis labels inside graticule on data frame of ArcMap?

texts: I'd like to add a graticule in DMS to my map in ArcMap 9.3, having tick marks and label axes inside the main figure/frame (similar to the below map). Although setting it for tick marks is very easy, I couldn't find a solution for the labels. Any thoughts?


texts: I have a problem with setting up zoom levels and max/min scale in QGIS web client. It looks like that values used in GlobalOptions.js file does not apply. Right now ZOOM_LEVELS is 70, maxScale 50 and minScale 40000.


texts: Spatialite ST_area(geom) won't compute for GeomObject

texts:
I extracted a dataset from Wikimapia via API using .kml download --> ogr2ogr .kml to GML -explodecollections --> GML to .shp only LINESTRING part of geometry collection.

Running QSpatiatlie in QGIS, I'm trying to run

SELECT name, ST_Area(geom)

but I'm getting only zeroes under the ST_Area(geom) column. No SQL errors.

Is this a geometry issue?


texts:
I am writing an ArcPy script that creates a Point Density Raster file, and then changes the following properties in the resulting layer:

Symbology

Classified

Classes = 10

Classification Method = Geometrical Interval

Excluded Value = 0
Display

Resample during display using:  bilinear interpolation

I know I can set lyr.symbology.numClasses = 10. and lyr.symbology.excludedValues = 0, and even create a dummy layer with the above properties (symbologyLayer) and arcpy.ApplySymbologyFromLayer_management(lyr, symbologyLayer).  Is there a way to do this simply by setting attributes of lyr.symbology, do I have to do it in one fell swoop with ApplySymbologyFromLayer_management, or is there a better way?


texts: Construct curve of equal distance from two straight lines

texts:
Given two straight lines defined by two points each, I want to get the line which separates points nearer to line a from points nearer to line b. The attached PostGIS code is an example case, it finds points of equal distance by creating a lot of buffers around the lines and intersecting them:

Is there an algorithm to calculate this line more direct, without sampling points or rastering?

I thought the line for this example would have a linear central part, but the dots are placed slightly under that line. Are they placed correctly?

CREATE TABLE tv_lines1 (
    id text PRIMARY KEY,
    geom geometry(linestring, 4326)); -- set meaningless CRS for easy handling in QGIS
INSERT INTO tv_lines1 (id, geom)
    SELECT 'a', ST_SetSRID(ST_MakeLine(ST_Point(0.0, 0.0), ST_Point(30.0, 0.0)), 4326);
INSERT INTO tv_lines1 (id, geom)
    SELECT 'b', ST_SetSRID(ST_MakeLine(ST_Point(16.0, 2.0), ST_Point(40.0, 12.0)), 4326);
CREATE INDEX ON tv_lines1 USING gist (geom);

-- create buffers in Steps of 0.1
CREATE TABLE tv_buffers (
    line_id text,
    size float,
    geom geometry(linestring, 4326),
    PRIMARY KEY (line_id, size));
INSERT INTO tv_buffers (line_id, size)
    SELECT
        id,
        generate_series(1, 1000) / 10.0
        FROM tv_lines1;
UPDATE tv_buffers bf
    SET geom = ST_ExteriorRing(ST_Buffer(l1.geom, size))
    FROM tv_lines1 l1
    WHERE l1.id = bf.line_id;
CREATE INDEX ON tv_buffers USING gist(geom);

-- find intersections
CREATE TABLE tv_intersections (
    id SERIAL PRIMARY KEY,
    distance float,
    geom geometry(point, 4326));
INSERT INTO tv_intersections (distance, geom)
    SELECT
        bf1.size,
        (ST_DUMP(ST_Intersection(bf1.geom, bf2.geom))).geom
        FROM tv_buffers bf1, tv_buffers bf2
        WHERE bf1.line_id = 'a' AND
            bf2.line_id = 'b' AND
            bf1.size = bf2.size;
CREATE INDEX ON tv_intersections USING gist(geom);


texts: Why will my raster only display at certain extents?

texts:
I have a large raster image which contains the florida keys and has 1 meter by 1 meter resolution.  When I change the extent and zoom in from 1:6,000 and closer, the image is displayed.  When I change the extent and zoom out from 1:8,000 and beyond, the image is black.

This is LiDAR data and my display troubles only occurred after I created a Mean raster (I created a new raster from the mean values of the pixels from the original LiDAR raster)

I set the pyramid structure to -1 for automatic default.


texts:
I'm trying to find out what tool would be best suited to address a problem  that I have.

I am trying to find the most efficient route (least amount of pipe (footage) / most amount of houses provided with fire hydrant coverage (500 ft diameter buffer) for a water pipe through a sub division.

I know this is very vague, but what tool would be best suited to this?

I know find route would find the shortest route, but I also want to provide the most amount of coverage provided by the fire hydrants, so shortest route isn't necessarily the best suited

At the moment I am just putting my own buffers on a path that I myself think would serve the most houses with the least amount of buffer overlap from the hydrants, but I'm frustrating myself because there most probably is a better way of doing it (that I cant think of at the moment!)


texts: ESRI Javascript API - Emulate 'Red Glow' Effect From Flex API

texts:
I am trying to emulate the red glow effect so common in the Flex API Examples.  Really, any type of blinking/flashing style marker will do.  I am assuming I will have to use a picture marker symbol with an animated gif?

I have probably infected my computer with who knows what downloading trial software to convert from swf to gif and all types of other non nonsensical conversions.  I am hoping someone can point me to an example or give me some insight on how to accomplish this.


texts: How can I average point features spatially in QGIS?

texts:
I have a GPX file with track points and I want to average the points spatially to improve accuracy of the point feature.  I can't just make a polygon and find its centroid as the number of overlapping points in a certain location should "weight" the result more than another location with few points.  ie. if the field data collector was standing in one spot for 2 minutes I should have the bulk of points 'drifting' around the likely location of the point, and have scattered outliers around that, by averaging spatially based on all the points the outliers would be discounted in favour of the clustered points in the averaging.

Do I have to make a raster grid with small resolution and overlay it on my points and then count the points in each cell, or is there an easier geoprocessing kind of function to simply tell me the center coordinates (or make a new point) for the average center?


texts: R gstat krige() - Covariance matrix singular at location [5.88,47.4,0]: skipping

texts:
When I want to perform the kriging it only works sometimes, depending on what values I use in my datatable. As a result of the krige function I get for var1.pred: NA NA NA ... and var1.var: NA NA NA ... (but only when I use the "wrong" values in my datatable.)

For example:

it works always (so far) when I only use 10 values
it works when I use 50 values, but only with certain ones
it works not when I use 50 values and the "wrong" values
it works when I use 25 values and the before mentioned "wrong" values

I don't get it why it's sometimes working and sometimes not. The strange thing is that when I add Zwiesel;49.02999878;13.22999954;2.2 to the datatable it works when I use less then ~20 values, but it doesn't work when I use more then 50 values ...

Where is my mistake?

myWeatherTable.csv:

Place;Latitude;Longitude;Temperature
Aachen;50.77999878;6.09999990;3
Abbikenhausen;53.52999878;8.00000000;7.9
Adelbach;49.04000092;9.76000023;3.1
Adendorf;51.61999893;11.69999981;1.9
Alberzell;48.45999908;11.34000015;4.6
...
...

My code to perform kriging interpolation

WeatherData <- read.csv(file="myWeatherTable", header = TRUE, sep ";")

coordinates(WeatherData) = ~Longitude + Latitude

vario <-  variogram(log(Temperature) ~1, WeatherData)
vario.fit <- fit.variogram(vario, vgm("Sph"))

min_lon <- min(WeatherData$Longitude)
max_lon <- max(WeatherData$Longitude)
min_lat <- min(WeatherData$Latitude)
max_lat <- max(WeatherData$Latitude)
Longitude.range <- as.numeric(c(min_lon,max_lon))
Latitude.range <- as.numeric(c(min_lat,max_lat))
grd <- expand.grid(Longitude = seq(from = Longitude.range[1], to = Longitude.range[2], by = 0.1),
  Latitude = seq(from = Latitude.range[1],to = Latitude.range[2], by = 0.1))
coordinates(grd) <- ~Longitude + Latitude
gridded(grd) <- TRUE

plot1 <- WU_data_spatial %>% as.data.frame %>%
  ggplot(aes(Longitude, Latitude)) + geom_point(size=1) + coord_equal() +
  ggtitle("Points with measurements")

plot2 <- grd %>% as.data.frame %>%
  ggplot(aes(Longitude, Latitude)) + geom_point(size=1) + coord_equal() +
  ggtitle("Points at which to estimate")

grid.arrange(plot1, plot2, ncol = 2)

kriged <- krige(Temperature~ 1, WeatherData, grd, model=variogram_fit)

Warnings:

1: In predict.gstat(g, newdata = newdata, block = block,  ... :
Covariance matrix singular at location [5.88,47.4,0]: skipping...
2: In predict.gstat(g, newdata = newdata, block = block,  ... :
Covariance matrix singular at location [5.98,47.4,0]: skipping...
3: In predict.gstat(g, newdata = newdata, block = block,  ... :
Covariance matrix singular at location [6.08,47.4,0]: skipping...
4: In predict.gstat(g, newdata = newdata, block = block,  ... :
Covariance matrix singular at location [6.18,47.4,0]: skipping...
...
...


texts: connection from ArcMAP to Postgresql, Is ArcSDE required?

texts:
I wan to make a direct connection form Arcmap to Postgresql.

I have ArcGIS for Desktop 10.1 on windows (64bit).

So I need to put "pg_client_windows64" libraries in C:\Program Files (x86)\ArcGIS\Desktop10.1\bin. Is this right?

Where to get "pg_client_windows64"libraries from?

Do I need to install ArcSDE?


texts: Is PhD in GIS useful in non-educational career path?

texts:
A bit of an open debate - I am curious for your opinions. I am considering commencing a PhD in GIS which may concern the application of surveys within the natural environments. Eventually, I wish to continue my career in professional consultancy and not to pursue academia beyond these studies.

I wonder whether anybody has any experience or can just share their thoughts about how useful is to have a GIS PhD in a professional industry. My concern is that people have told me that for some professions (e.g. Law), having a PhD qualification can in fact reduce your employability.

How is it seen by industry employers to have a PhD in GIS? Also, are there any specific sites or tips anyone can suggest?


texts:
I am new to the Quantum Plugin development and wonder if there is more information about the following.

I have added a shapefile and am (off course) able to select a point from this shapefile. In the shapefile I have stored the location (x,y) and an id.

Now I would like to be able to select a point in this layer and then my plugin should be able to recognize the id of the point, look up a database and display a graphical representation of the data in the database.

So what I am looking for is this,

how do I get the information of a selected point to my plugin
is it possible that my plugin overwrites the default information
dialog (so the plugin pops up instead of the normal dialog)
is there perhaps a plugin that already shows such behaviour and I
can learn from (haven't found one so far)


texts: Reducing space in geopandas and matplotlib.pyplot's subplots?

texts:
I use geopandas and matplotlib.pyplot's subplots to plot two subplots in a single figure, with a single colorbar, as:

How do I reduce the whitespace around the maps in each subplot (not in between subplots - I know how to do that)?

The relevant code is:

import geopandas as gpd
import matplotlib.pyplot as plt

side_length = 10
test_data = gpd.read_file(shpfile)
fig, axes = plt.subplots(1, 2, figsize=(side_length * 2.1, side_length))
fig.subplots_adjust(hspace=0.5, wspace=0.5)

# 1st subplot
test_data['data'] = data1 # assume data is in a list data1
ax = test_data.plot(column='data', axes=axes.flat[0], cmap=cmap, edgecolor='gray', vmin=vmin, vmax=vmax)
ax.set_axis_off()
ax.set_aspect('equal')

# 2nd subplot
test_data['data'] = data1 # assume data is in a list data2
ax = test_data.plot(column='data', axes=axes.flat[1], cmap=cmap, edgecolor='gray', vmin=vmin, vmax=vmax)
ax.set_axis_off()
ax.set_aspect('equal')

# Add single colorbar
cax = fig.add_axes([0.9, 0.17, 0.02, 0.65])
sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm._A = []
cbar = fig.colorbar(sm, cax=cax)
cbar.ax.set_title('Number of Cases')

fig.savefig('myFig.pdf', bbox_inches='tight', pad_inches=0.1)


texts:
I am following these instructions in order to define my projection for my map. And then use the ol.proj.transform function to transform the coordinates. So my code snippet is the following:

  var myProjection = new ol.proj.Projection({
      code: 'EPSG:2100',
      extent: [104022.946289, 3850785.500488, 1007956.563293, 4624047.765686],
      units: 'm'
       });

    ol.proj.addProjection(myProjections);

    var center = new ol.proj.transform([21.767997, 39.556105], 'EPSG:4326', 'EPSG:2100');

However if I try to console.log(center) the coordinates are not transformed. I'm still getting the same coordinates which are

21.767997, 39.556105

Does this mean that the addProjection didn't work??


texts:
I want to create a sortable lists of layers using jquery ui. I'm trying to use .setZindex() and  .index() method. I think that index of
 does not refresh after order changing. I'm not famililar with jquery any help?

here is the basic concept:

<!doctype HTML>

    <meta charset="utf-8">

    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />

    <link rel="stylesheet" href="jquery-ui.min.css">

Layers

rivers

landuse


texts: Proper assignment of QML legend to raster failed in QGIS 3.0.1

texts: I downloaded raster data sets together with the legend qml files from the SoilGrids page (ftp://ftp.soilgrids.org/). Assignment of the qml legend (TAXOUSDA_1km/TAXN_WRB_1km) didn't work in QGIS 3, as the legend was discretized in floating-point numbers and not integers (to represent the xx unique identifiers of the raster).  However, in QGIS 2.18. loading and proper assignment of the legend worked. Thus, I was thinking that this might be a technical issue of the new QGIS version, which could be addressed/solved in the future. Or is there another way to solve the problem in the legend/qml properties?


texts: Update Hosted ArcGIS Feature Service Schema Issue

texts:
I am trying to update a field length of a hosted feature service in ArcGIS Online.  I load the update definition page by going to the Admin REST page of the service, and clicking Update Definition:

I then take the JSON and copy/paste it into an editor, like Notepad++.

I find the field I want to update, and change the length.  In this case I want the Category field to have a length of 70, rather than 50:

I then copy/paste the JSON back to the Update Definition page and click Update Service Definition, and receive a success message:

But when I go back to the REST page, CATEGORY is still 50 in length:

Not sure what the issue is here.  Do I need to set ASync to True?


texts: Subset a lascatalog point cloud within a function while calculating grid_metrics

texts:
I am trying to create subsets of a point cloud based on the ScanAngleRank and compute some metrics for each subset. Then return a stack of values(rasters) for the metric. I modified the rumple index example as shown below:

met_calc = function(cluster, res)
{
  las = readLAS(cluster)
  if (is.empty(las)) return(NULL)

  outlist <- list()
  for(i in y){
    ls <- lasfilter(las, abs(ScanAngleRank) >= min(i) & abs(ScanAngleRank) < max(i))
    if (is.empty(ls)) return(NULL)
    out <- grid_metrics(ls,~mean(Z),res)
    outlist <- c(outlist,out)
  }
  #bbox   <- raster::extent(cluster)
  #outlist <- raster::crop(outlist, bbox)
  return(outlist)
}

opt_chunk_buffer(ldr) <- 1
opt_chunk_size(ldr) <- 300
opt_select(ldr) <- "*"
opt <- list(raster_alignment = 20, automerge = TRUE)
output <- catalog_apply(ldr, met_calc, res = 20, .options = opt)

y <- list(0:10, 10:20, 20:30, 30:40, 40:50)

In the first iteration, for example, I tried to subset las for minimum and maximum scan angle of 0 and 10 degerees respectively, and so on. In this case, that means 5 rasters.

bbox and outlist are comments because they were generating errors.

I get the following warning:

The list returned by 'catalog_apply' contains unsupported objects. Merging is impossible. A list has been returned.

Output is a complicated list of RasterLayers for small areas.

I tried using grid_metrics() to return a list of metrics as shown in the example in the documentation. Have not been able to figure that.

I understand how to work with a single las file. My confusion is with the usage of LAScatalog and related options for this purpose, efficiently.


texts: I would like to draw a line between few GPS coordination. I tried to use OpenLayers.Geometry.LineString for that, but it uses a OpenLayers.Geometry.Point. Is there any  function in OpenLayers that converts GPS coordination (latlon) to Geometry.Point? Or maybe I'm doing this wrong and shouldn't use lineString for drawing a line between latlon points?


texts: Changing the projection transformation for a feature class after adding it to an mxd

texts: Is there a way to select a different projection transformation after you've added a feature class to an mxd?  I have a pretty complex mxd set up to use as a map service where the original data is in CA State Plane and the mxd is set up as web mercator.  When I use this service on top of another service where the data is originally in web mercator, the datasets don't line up due to the transformation selected originally.  I'm trying to avoid having to re-add all of the datasets and set them all up again.


texts:
I'm trying to change the values of nodata from my source file to a set value for my output file, but all I am getting is an error.
This is my code:

    from osgeo import gdal
    from subprocess import call

    input_raster = r"c:\blah\someinputtiff.tif" #This file I want to clip

    output_raster = r"c:\blah\someoutputtiff.tif" #The resulting clip file

    input_shape = r"c:\blah\someshapefile.shp" #The file I want to use to clip the.tif

    call(['C:\Program Files (x86)\GDAL\gdalwarp.exe','-
    dstNodata"value[-9999]"', '-cutline' , input_shape , '-crop_to_cutline' ,
    input_raster , output_raster, '-overwrite'])

This produces the following error:

    FAILURE: Unknown option name '-dstNodata"value[-9999]"'

I have tried multiple permutations of '-dstNodata' followed by a value but no luck.


texts: Points in Polygon missing from QGIS 2.16.x?

texts:
Looking to do a basic Hex analyis with some random points within the Gulf of Mexico.

I created a hex grid from MMQGIS

Under my Vector tool - Analysis tools, I dont have Points in Polygon.

Is it missing on 2.16.3 (2.16.x)?


texts: How to calculate distances in QGIS field calculator?

texts:
So I'm using qgis 2.6.0 and I have a csv (converted into a .shp) with two geometry columns, I want to create a new column with the distance between the two geometry:

shot.lat    shot.lon    haul.lat    haul.lon
57.563      -12.986     57.534      -12.985
57.7        -13.304     57.67       -13.299
57.611      -13.529     57.582      -13.538
57.594      -13.294     57.565      -13.315
57.432      -13.314     57.402      -13.304
57.339      -13.053     57.311      -13.072

I am trying to use the function geometry function "distance" with the following expression:

distance( geomFromWKT( 'POINT("shot.lon" "shot.lat")' ) , geomFromWKT('POINT( "haul.lon" "haul.lat")' ))

This should then create a new field with the distance between the points. When it runs though, it just returns Null. Am I missing something obvious?


texts: Which free/cheap tools can open ESRI Geodatabase Files?

texts:
Ok I don't need anything advanced, I just don't know anything about this GIS business.  I'm trying to get local parcel information and the government site states it is in "ESRI geodatabase format."  I've tried Google and get a wiki article that seems to be an advertisement and an esri.com that won't load a single page... so I'm stuck for information.  I don't know what software I would need, and I'm really not hoping to plunk down some ridiculous amount of cash for a proprietary format our lovely gov't reps decided to use...

Really, I just want to be able to correlate parcel listings with a location on a human (me) readable map.  An added [super helpful] bonus would be the capability to feed in a large listing of parcels and see their locations at the same time, as an overview.

Any ideas/help?


texts: GeoServer layer doesn't show up if admin is logged out

texts:
I imported a shapefile successfully to GeoServer.
When I query the WMS service for its capabilities, i.e.:
https://<domain>/geoserver/ows?service=wms&version=1.3.0&request=GetCapabilities
If the user is logged in as admin, the layer itself and all its capabilities show up fine.
If the user is logged out (using another browser, a private window, another computer, etc...), the Layer simply doesn't seem to show up.
A more specific example - When I specifically query for the layer while having an admin session, i.e. https://<domain>/geoserver/wms?service=WMS&version=1.1.0&request=GetMap&layers=<my layer> [other parameters omitted for brevity]
The layer is shown just fine if I have admin logged in.
If the user is logged out, I get the error Could not find layer.
A few more details about my setup:
My GeoServer is based on the Docker image: https://hub.docker.com/r/geonode/geoserver - Version 2.18.2
I checked security/layer.properties and they seem to be:
#Thu May 20 10:14:42 GMT 2021
*.*.r=*
*.*.w=GROUP_ADMIN,ADMIN
mode=MIXED

Which I think means "Anyone" can read any layers, right?
Does anyone have any idea what could be causing this? Is there any privacy parameter that I might have overlooked? Any specific Docker parameter that I missed? I ask this because this doesn't seem to happen with my other non-Docker deployment.


texts: Why is there a shift between OpenLayers and vector layer in print?

texts:
Sorry for my poor english.
I'm a very recent user of QGIS (QGIS 1.8; ubuntu 12.04).
Does anybody knows why the .tif or .jpeg or .pdf files I generate with the QGIS composer have different layers completely shifted?

For example, 2 layers: 1 open street map landscape and 1 layer of my GPS points.
In the composer window, the points are exactly at the right place, the map is nice and ready to be saved or printed...

But when I do export to pdf or .jpeg or .tif the picture that I obtain shows the map very small in a corner and the GPS point in the middle of a white picture. When I try to print, a white sheet of paper is produced.

Any ideas?

thank you very much


texts:
I have a list of coordinates in Lambert Conformal Conic projection. Those coordinates draw the border of specific area.

  (179506,161069),(180016,161633),(180375,162142),...

I need to create a topojson file to use this as a map in D3.js

I found how to insert them in QGIS but after that I do not understand how to generate my topojson.

Can I do it directly in QGIS or is there an other way ?


texts: How do I create an attribute table in a web page from GeoServer using JavaScript?

texts:
I'm trying to create a custom web application using Leaflet and Geoserver. I've got my layers visible in a map viewer, but I'd also like to show an attribute table of the visible map layers. I can't seem to figure this out.

I tried going this route:

var url = "http://serveraddress/geoserver/store/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=store:layer&outputFormat=text/javascript&format_options=callback:processJSON”";

    $.ajax({
        jsonp: false,
        jsonpCallback: 'processJSON”',
        type: 'GET',
        url: url,
        async: false,
        dataType: 'jsonp',
        success: function(data) {
           // useMyJson(data);
        }
    });

But it just crashes my browser, probably from data overload...

Any ideas?


texts: I am looking for a single (paid) provider of postal (zip) code shapes (boundaries) for countries in Europe. I hope to find reliable partner with reasonable pricing and update policy. What I need are polygons (shp and shx, or geojson format) to draw service coverages and host it in commercial web application.


texts: Drawing Circles with set radius using QGIS

texts:
I have a site and want to draw an indicative line in a 5km buffer around the centre of the site. So I want to draw a circle with a radius of 5km. I have a tool for drawing circles and ellipses and rectangles squares and spirals.
None of these let me state the centre point and the radius.


texts: Problem loading raster to PostGIS using QGIS

texts:
I am having trouble loading a raster (.asc) file to PostGIS 2.0.1 using the plugin Database - PostGIS Raster - Load raster to postgis in QGIS 1.8.0 on Win Vista 64 bit. I get the following error message:

  Ett fel har inträffat vid Pythonexekvering:

  Traceback (most recent call last):   File
  "C:/Users/Magnus/.qgis//python/plugins\wktraster\DlgRasterLoader.py",
  line 126, in browseRaster
      self.getMetadata(str(fileName))   File "C:/Users/Magnus/.qgis//python/plugins\wktraster\DlgRasterLoader.py",
  line 163, in getMetadata
      ds=gdal.Open(filename) AttributeError: 'module' object has no attribute 'Open'

  Python version:
  2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]

  QGIS version:
  1.8.0-Lisboa Lisboa, 6416f38

  Python sökväg: ['C:/Users/Magnus/.qgis//python/plugins\sextante',
  'C:/PROGRA~2/GIS/QGIS/apps/qgis/./python',
  'C:/Users/Magnus/.qgis//python',
  'C:/Users/Magnus/.qgis//python/plugins',
  'C:/PROGRA~2/GIS/QGIS/apps/qgis/./python/plugins',
  'C:\PROGRA~2\GIS\QGIS\bin\python27.zip',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\DLLs',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\plat-win',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\lib-tk',
  'C:\PROGRA~2\GIS\QGIS\apps\qgis\bin',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages\PIL',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages\win32',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages\win32\lib',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages\Pythonwin',
  'C:\PROGRA~2\GIS\QGIS\apps\Python27\lib\site-packages\wx-2.8-msw-unicode',
  'C:\Users\Magnus\.qgis\python\plugins\mmqgis/forms',
  'C:\PROGRA~2\GIS\QGIS\apps\qgis\python\plugins\fTools\tools']

QGIS opens the file just fine. Any suggestions on what to do?

Magnus


texts: Openlayers - Vector Layer - Combining hide/show features using style property

texts:
I am trying to filter vector layer features according to their attributes. I have 4 check boxes: type_1, type_2, shape_1 and shape_2 Note that I am using Extjs for interface.

at beginning I could filter type attribute setting style to 'none' or '' ... like so :

switch (f) {
case 'type_1':
if(checked)
filter('show_type_1');
else filter ('hide_type_1);

case 'type_2':
if(checked)
filter('show_type_2');
else filter ('hide_type_2);

}

function filter(value)
{
 switch (value){

case 'hide_type_1':
for (i=0;i<=features.length;i++)
if(features[i].attributes.type == 'first')
features[i].style = 'none';
layer.redraw();

case 'show_type_1':
for (i=0;i<=features.length;i++)
if(features[i].attributes.type == 'first')
features[i].style = '';
layer.redraw();

case 'hide_type_2':
for (i=0;i<=features.length;i++)
if(features[i].attributes.type == 'second')
features[i].style = 'none';
layer.redraw();

case 'show_type_2':
for (i=0;i<=features.length;i++)
if(features[i].attributes.type == 'second')
features[i].style = '';
layer.redraw();

}
}

all above works great, if I uncheck type_1, all type_1 features will diappear. then if I uncheck type_2, all type_2 features will diappear. then if I check type_1 again, all type_1 features will appear, and type_2 features will stay hidden.

The problem is when I tried to do the same with shape attribute, by adding:

case 'shape_1':
if(checked)
filter('show_shape_1');
else filter ('hide_shape_1);

to the first function.

and:

case 'hide_shape_1':
for (i=0;i<=features.length;i++)
if(features[i].attributes.shape == 'first')
features[i].style = 'none';
layer.redraw();

case 'show_shape_1':
for (i=0;i<=features.length;i++)
if(features[i].attributes.shape == 'first')
features[i].style = '';
layer.redraw();

to the second one.

it works when I uncheck the shape_1 checkbox. all shape_1 features hid. but if check it to display them, all features will be displayed even type_1 and type_2 which I was keeping them unchecked and hidden!

I don't understand why that happens. because it works fine if I handle one attribute (type in my case) but failed when tried to combine type filtering with another feature attribute (shape).

I hope my explanation is clear.


texts:
The code snippet below creates a map variable locally inside an init() function.  I am getting errors such as this.map.getLayers and map.getLayers is not a function or cannot read property getLayers of undefined, depending on how I call the function updateFliter()

If I remove the var in front of map, js creates a global variable, and updateFilter() works.

How can I call  updateFliter() without creating map as a global variable?

FYI, updateFilter() is implementing a CQL filter on the geoserver layers.

function init() {
    var option1 = document.createElement("option");
    var option2 = document.createElement("option");
    option1.text = "layer1";
    option2.text = "layer2";
    document.getElementById("dropdown").add(option1);
    document.getElementById("dropdown").add(option2);

    var map = new ol.Map({
        target: 'map',
        layers: [
           new ol.layer.Tile({
                source: new ol.source.OSM(),
                name: 'OpenStreetMap'
            }),
            new ol.layer.Tile({
                source: new ol.source.TileWMS({
                    url: '../../geoserver/foo/wms',
                    params: {'VERSION': '1.1.1',
                        tiled: true,
                        LAYERS: 'foo:layer1',
                        tilesOrigin: -181.800003051758     + "," + -66.7399978637695
                   }
               }), name: 'layer1'
            }),
            new ol.layer.Tile({
                source: new ol.source.TileWMS({
                    url: '../../geoserver/foo/wms',
                    params: {'VERSION': '1.1.1',
                        tiled: true,
                        LAYERS: 'foo:layer2',
                        tilesOrigin: -181.800003051758     + "," + -66.5354995727539,
                   }
               }), name: 'layer2'
            })
          ]
        })

    var updateFunc =     document.getElementById('updateFilterButton');
    updateFunc.addEventListener('click',     function(){updateFilter(map)});
}

//Use a CQL filter to select data to display
function updateFilter(map){
    var filterParams = {'CQL_FILTER': null};
    var filter =     document.getElementById('CQLfilter').value;
         //by default, reset all filters
    if (filter.replace(/^\s\s*/, '').replace(/\s\s*$/,     '') != "") {
        filterParams["CQL_FILTER"] = filter;
    }
        this.map.getLayers().forEach(function(lyr) {
          if (lyr.get('name') == dropdown.value){
          lyr.getSource().updateParams(filterParams);
          };
        });
}


texts:
I try the new filter expression with QGIS 2.8. I want to filter only the selected objects on a layer in my map. My expression is:

import processing
from qgis.core import *
vlayer = processing.getObject('mylayer')
selected_objects = vlayer.selectedFeatures()`

I save this expression but I can't open QGIS after.

Error: "AttributeError: 'NoneType' object has no attribute 'selectedFeatures"


texts: Open Source .Net GIS software end-to-end

texts:
I need to provide a client with an end-to-end web-based GIS solution running on IIS. I will admit to not knowing much about the subject, other than having developed a Google Earth enabled version of the project.

However, I need to develop a version that works offline, in the sense that I have an internal network without internet connection.

My problem is twofold:

I don't know what I need, or rather I don't know how many separate layers make up a GIS server. For instance, I thought QGIS Server by itself would be sufficient, but I still need to get my hands on a spatial data store (or do I?) and it seems that QGIS will only run on a Windows environment if IIS is disabled.
I don't know where to get my hands on the various components I don't know I need.

Has anyone implemented a .NET web-based GIS service utilizing Open Source technology?

UPDATE: The features I need are:

Javascript API
Ability to easily import custom tiles, e.g. for places like Angola that don't have low-level tiles
Ability to add and remove place markers with balloons(like google maps pins with info window popup) using javascrip API


texts: I'm trying to merge a couple of GeoTiff maps using Raster > Miscellaneous > Merge but the merge is only coming out in grayscale. The geotifs are using palettes color. How do I get the merge to keep the color information?


texts: Fill mbtiles with hatch instead of single colour

texts:
I'm using tileserver-gl to serve tiles. At the moment I am using fill-color to fill with one colour.

"fill-color": "rgba(230, 36, 0, 1)",

But now I need to fill with stripes, e.g rgba(230, 36, 0, 1), rgba(11, 111, 202) and I can't figure out how to do this.

I tried creating an svg stripes.svg and replacing fill-color with fill-pattern, but when I do that nothing displays. I'd like to find out how to get my stripes.svg to display, or if there's a different approach I could be using to get the hatched fill.

original my_style.json (solid fill)

{
    "id": "my_style",
    "version": 1,
    "name": "my_style",
    "sources": {
        "mbtiles1": {
            "type": "vector",
            "url": "mbtiles://my_mbtiles.mbtiles"
        }
    },
    "layers": [
        {
            "id": "layer1",
            "type": "fill",
            "source": "mbtiles1",
            "source-layer": "cat1",
            "paint": {
                "outline-color": "rgba(230, 36, 0, 1)",
                "fill-color": "rgba(230, 36, 0, 1)",
                "fill-translate-anchor": "map"
            }
        }
    ]
}

updated my_style.json (attempt to hatch fill, but nothing displays)

{
    "id": "my_style",
    "version": 1,
    "name": "my_style",
    "sources": {
        "mbtiles1": {
            "type": "vector",
            "url": "mbtiles://my_mbtiles.mbtiles"
        }
    },
    "layers": [
        {
            "id": "layer1",
            "type": "fill",
            "source": "mbtiles1",
            "source-layer": "cat1",
            "paint": {
                "fill-pattern": "stripes.svg"
                "fill-translate-anchor": "map"
            }
        }
    ]
}

config.json

{
    "options": {
        "paths": {
            "root": "/data",
            "fonts": "fonts",
            "styles": "styles",
            "mbtiles": "/data"
        }
    },
    "styles": {
        "my_style": {
            "style": "my_style.json",
            "tilejson": {
                "bounds": [
                    -180,
                    -85.06,
                    180,
                    85.06
                ]
            }
        },
    }
}

Directory layout

├── my_mbtiles.mbtiles
├── config.json
├── fonts
│   └── __empty__dir__
├── stripes.svg
└── styles
    └── my_style.json


texts:
I am trying to add the PyQGIS in PyScripter, I have followed the tutorial provided in QGIS tutorial: IDE settings for writing and debugging plugins. I have created the bat file as follows:

SET OSGEO4W_ROOT=C:\OSGeo4W
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
SET QGIS_PREFIX_PATH=%QGIS%
SET PyScripter="C:\Program Files\PyScripter\PyScripter.exe"

CALL %OSGEO4W_ROOT%\bin\o4w_env.bat

SET PATH=%PATH%;%QGIS%\bin
SET PYTHONPATH=%QGIS%\python;%PYTHONPATH%

REM PYTHONHOME: %OSGEO4W_ROOT%\bin
REM PYTHONPATH: %OSGEO4W_ROOT%\apps\Python37\lib;%OSGEO4W_ROOT%\apps\qgis\python;%OSGEO4W_ROOT%\apps\Python37\Lib\site-packages

start "script" /B %PyScripter% %*"

The error message appears again:ModuleNotFoundError: No module named 'qgis'.

Note: I am using QGIS 3.10  and python 37.


texts: Creating points projected on line using QGIS?

texts:
How can we create points perpendicularly on a line with QGis ?

Otherwise, how can we draw a line between one point and a line ?

Here is a picture to show what we want :


texts: Seeking function to calculate mean for missing or NA values in dataframe using R?

texts:
I am trying to develop a function that will calculate a mean lat/long for my missing or zeros location (lat/long) data in my dataset.

I have locations for many SETS within the same TRIP, which tends to be close in space. However, in many cases, I am missing the location information for some of my SETS. I would like to replace these NAs or Zeros values with the mean lat/long calculated from the remaining SETS within that TRIP.

> head(newdata[,c(3,12,26:27)])#This is my data
        TRIP_ID SET_NO SH_LATITUDE SH_LONGITUDE
13365 100020848      3    41.88983    -66.29183
13373 100020848      6    41.93317    -66.29767
13430 100020848     11    42.01200    -66.45217
13442 100020848      7    41.92383    -66.29733
23207 100020848      4    41.94817    -66.29750
24669 100020848      9    41.88300    -66.30567

So for example, SH_LONGITUDE has 17239 SET_NO with missing values from different TRIP_ID.

> look<-newdata[newdata$SH_LONGITUDE==0,]
> dim(look)
[1] 17239   130
> head(table(look$TRIP_ID))

100020997 100021109 100021204 100021306 100021337 100021340
        1         2         1         1         2         1

I am a newbie with user-defined functions, but this is what I would like to be able to call for.

replacebymean<-function(dat, x...?) {
    if dat[dat$x==0, ] or dat[is.na(dat$x), ]
    then "find these NAs/zeros associated TRIP_ID" and
          calculate the mean lat/long for each TRIP_ID as
             mean(dat$x[dat$TRIP_ID=="e.g.,100022478"], na.rm=T) and
          replace these NAs/zeros by this mean value.
}


texts: AutoCAD Map import MapInfo TAB file georeferencing info

texts: Is there a way to import the TAB file created for georeferencing of an image by Mapinfo into AutoCAD so I don't have to do all the georeferencing again?


texts: Why does the order of queries matter in overpass?

texts:
I am new to overpass / OSM. I want to locate Dodo Pizza in Moscow, a specific chain of fastfood. Since these restaurants are sometimes tagged by brand and sometimes by name, I ve included a or statement in my query as explained here. I then realized that the returned results are different depending on the order of my conditions.

If I use first the tag brand, overpass doesn't return any matches

[out:json][timeout:2500];
area["name"="Москва"];
(

node(area)["brand"="Додо Пицца"];
node(area)["name"="Додо Пицца"];

);
out body;

But if I use first the tag name, overpass return 20 matches

[out:json][timeout:2500];
area["name"="Москва"];
(

node(area)["name"="Додо Пицца"];
node(area)["brand"="Додо Пицца"];

);
out body;

I am so confused, what am I doing wrong?


texts:
I am currently working on the different coordinate system.
Specifically, I am concerned with the correct presentation of EPSG:25833 (ETRS89/UTMZone89) in Leaflet.

If I am right Leaflet used as default this system:

https://spatialreference.org/ref/epsg/wgs-84/

As an example I use this WMS Server https://www.orka-mv.de/geodienste/orkamv/wms?version=1.3.0&service=WMS&request=GetCapabilities that uses this system:

https://spatialreference.org/ref/epsg/25833/

I had a look to Proj4Leaflet (https://kartena.github.io/Proj4Leaflet/) that makes it possible to use projections and CRS not built into Leaflet. All major projections (https://en.wikipedia.org/wiki/List_of_map_projections) are supported by Proj4Leaflet.

Then I created a WMS layer without projection and in Leaflet everything looks OK. This is my code:

var testlayer = L.tileLayer.wms("https://www.orka-mv.de/geodienste/orkamv/wms?version%3D1.3.0%26", {
    format: 'image/png',
    transparent: true,
    opacity: 0.6,
    layers: 'orkamv',
    attribution: "todo"
}).addTo(map);

I then looked at the exact differences and found that the "difference/error" is less than 1 meter.

I think about whether the effort of a projection for display is worthwhile.

But, it is important for me to output the coordinates correctly. I think I can do that by converting the standard Leaflet coordinates (for example, 53.22 12.18) into my format (for example, 347560 5922583) using a formula.
If I see that correctly, I must also exchange lat and lon

My question: Whats the best way to convert EPSG:25833 to WGS:84 in Leaflet? Is there something I can rely on or need to reinvent the wheel myself


texts: Using a colour wheel as a colour ramp for displaying directional data in Arcmap

texts:
I would like to define a colour ramp for displaying direction data in a 'smooth' fashion. So for data in degrees there wouldn't be a jump in colour between 359 and 0 degrees. It seems to me a colour wheel would be great for this.  I haven't seen this done, the preset ramps in arc don't do it, but I would think someone must have.

My question is how to implement this. My thought is to combine a ramp from red through yellow to green with one from green through blue to red (or whatever the proper quadrant colours are) and map that to 0 through 90 to 180 and 180 through 270 to 360. But I don't know how to set this up. Any ideas on how to define the ramps?


texts: PHP MapScript API Development for MapServer

texts:
PHP MapScript Tutorial
I'd like to use the PHP MapScript API and am curious how others have successfully used it to produce a site that serves up remote sensing data.  I've got land surface temperatures, emissivities, and reflectances I'd like to share via a mapserver.  I've already installed mapserver on an apache webserver and have also installed the MapScript shared library.  What is a good development environment to use for PHP MapScript development?  Any example sites that use PHP MapScript people know of?  How useful is the MapScript API and how easy is it to use?  Thanks in advance!


texts:
I have worked with 2D geometry for CAD - My questions are general and relate to understanding operations on geometric entities (lines).

Should geometric operations such as line-line intersection be done in the projected spaces (such as EPSG  3857) ?

I've had a brief look at the geos library, but it is unclear whether the inputs should be in a linear projected space.

If operations are done in projected space - does this introduce significant error in the result when the output is un-projected?


texts: Fixing dangles (lines) with v.clean and "fix geometries"

texts:
I need help with fixing my two dangle-errors:

there are 2 dangles on the top right and one small hole on the top-left corner.

I can deduct them with the Geometry Checker, but I cannot fix it (automatically).
This is only an example-file. I need the skill for a big file with hundreds of lines that need fixing before polygonizing...
Has someone an idea why v.clean (I ticked line and break, snap, rmdangle with tresholds 0,5,5) or geometry fixer cannot do that?
I already studied useful information like that:

https://www.geodose.com/2017/10/how-to-clean-topology-error-in-qgis.html
https://www.geodose.com/2017/09/how-to-check-topology-error-in-qgis.html
Finding Invalid Geometry using QGIS?
https://www.qgistutorials.com/en/docs/3/handling_invalid_geometries.html
What am I doing wrong with QGIS v.clean?
https://www.qgistutorials.com/en/docs/3/handling_invalid_geometries.html

... but still it is not working.


texts: PostGIS UPDATE Statement with ST_Contains(table1.geom, table2.geom)

texts:
I'm strugeling with a query I should write to UPDATE my dataset.

First of all, I have 2 tables containing data and I need to mix it.

Tabel 1 -- lets call it names:

name_str
the_geom (points)

Table 2 -- lets call it objects:

the_geom (polygons)

Now i want to do a cross selection and update my objects table. If a point is in the polygon it should update the objects.name_str with the value in the names.name_str.

Please keep in mind that I have over 2 mil datasets in the table!

I tried a selection query like this without success:

SELECT * FROM name, object WHERE st_contains(the_geom, way)='t' LIMIT 10

Please help! I have no clue where the error is hidden!


texts:
I'd like to install QGIS. I'm on  a system that defaults to Python 2.5x and I downloaded python 3.7x with homebrew. I created a virtual environment with this code.

brew install pyenv
brew install pyenv-virtualenv
exec $SHELL
eval "$(pyenv init -)"
eval "$(pyenv virtualenv-init -)"
pyenv activate my-virtualenv

I install python 3.6 and the modules plotly, numpy, pygments, psycopg2, owslib, and PyYaml.

Each of them are confined to the virtual environment. How do I install the package so that it loads and uses the Python from only the partcular environment.


texts: Why does OGR driver.Open() for File Geodatabase give AttributeError: NoneType object has no attribute Open?

texts:
There are a number of existing Questions that hint at a possible Answer to this Question, but none seem to be providing me with the solution to my problem:

How to access feature classes in file geodatabases with Python and GDAL?
How to add GDAL to existing Python installation?
Gdal Python Bindings and File Geodatabase API

To get GDAL/OGR 1.10 working on a Windows 7 SP1 (64bit) laptop with Python 2.7.3 from ArcGIS 10.2 for Desktop I did the following steps:

Installed ArcGIS 10.2 for Desktop which installed Python 2.7.3 - there are no other versions of Python present
Downloaded and installed gdal-110-1600-x64-core.msi (compiled 2/6/2014 6:11 AM) from http://www.gisinternals.com/sdk/PackageList.aspx?file=release-1600-x64-gdal-mapserver
Downloaded and installed gdal-110-1600-x64-filegdb.msi (compiled 2/6/2014 6:17 AM) from http://www.gisinternals.com/sdk/PackageList.aspx?file=release-1600-x64-gdal-mapserver.  As expected this placed an ogr_FileGDB.dll in C:\Program Files\GDAL\gdalplugins.
Added a new system environment variable called GDAL_DRIVER_PATH and set it to C:\Program Files\GDAL\gdalplugins, and I appended (using ';') the same folder onto the end of my Path variable.  This was from advice offered in one of the Questions cited above.
Downloaded and installed GDAL-1.10.1.win32-py2.7.exe from http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal

This enabled me to run some Python scripts from IDLE like the one below to count the number of features in a shapefile.

import os
from osgeo import ogr

daShapefile = r"C:\temp\test.shp"

driver = ogr.GetDriverByName('ESRI Shapefile')

dataSource = driver.Open(daShapefile, 0) # 0 means read-only. 1 means writeable.

# Check to see if shapefile is found.
if dataSource is None:
    print 'Could not open %s' % (daShapefile)
else:
    print 'Opened %s' % (daShapefile)
    layer = dataSource.GetLayer()
    featureCount = layer.GetFeatureCount()
    print "Number of features in %s: %d" % (os.path.basename(daShapefile),featureCount)

However, when I tried to run the code below:

import os
from osgeo import ogr

daFileGeodatabaseFC = r"C:\temp\test.gdb\testFC"
daFileGeodatabase = r"C:\temp\test.gdb"
featClass = "testFC"

driver = ogr.GetDriverByName('FileGDB')

dataSource = driver.Open(daFileGeodatabase, 0) # 0 means read-only. 1 means writeable.

# Check to see if feature class is found.
if dataSource is None:
    print 'Could not open %s' % (daFileGeodatabase)
else:
    print 'Opened %s' % (daFileGeodatabase)
    layer = dataSource.GetLayer(featClass)
    featureCount = layer.GetFeatureCount()
    print "Number of features in %s: %d" % (os.path.basename(daFileGeodatabaseFC),featureCount)

I got an error:

  Traceback (most recent call last):   File
  "C:\temp\GetFileGDBFeatureCount.py", line 17, in
      dataSource = driver.Open(daFileGeodatabase, 0) # 0 means read-only. 1 means writeable. AttributeError: 'NoneType' object has no
  attribute 'Open'

After some more reading I decided to use ogrinfo --formats in a DOS window to see if FileGDB was recognized and it is!

Am I missing a step or two that would complete enabling OGR FileGDB support for my GDAL/Python/ArcGIS configuration?

I do not class myself as a developer so compiling a new version of GDAL, or starting from scratch and installing OSGeo4W instead are options I would prefer to avoid, if at all possible.


texts: How to rotate labels according to an angle stored in a field in the attribute table?

texts: I have a street name layer. I want to have the name of each street be fitted exactly into the street shape. There is a field in the attribute table which has a degree of rotation, which is the angle at which each label must be placed to be inside the street shape. How is it possible to do this with Python programming?


texts: Iterate feature class by Date created to select recent shapefile

texts:
I would like to select only newly created shapefile (today's). Every time new shapefiles are created from Trimble pathfinder office which I need to upload to ArcGIS server feature class.

Model below iterates all shapefiles and appends but I need to append only newly created so as not to duplicate the data.


texts:
I am running into an issue where the arcpy.Exists() function returns a false for a known good layer in the python interactive window.  We are using a map document template for the layers set up in the Table of Contents.

I usually would drag the layer from the Table Of Contents into the interactive window so it shows like r'Work Layers\Replacement Tree' since it is part of a group in the Table of Contents.  If I manually add in the data instead of use a template and resource the layers, the arcpy.Exists() function works and returns true, but on occasion, it will return a false.

When it returns false, I can check the layer's source, navigate to the location and specified feature class, and see that it is there.  I can add it in manually and run it, too.

For some reason, it seems using an existing map document as a template and resourcing the Table of Contents breaks this on occasion.  The code I am using is a bit more complex than the arcpy.Exists(), but I have narrowed it down to that being the issue.  My code throws an exception saying a particular layer doesn't exist so that is the easiest way to troubleshoot.

I have even gone as far as to drag multiple layers from the Table Of Contents into a list in the interactive window and iterate over it with arcpy.Exists() to the same result of false.

I can set the workspace to the file geodatabase we're using and the using the feature class name (ie. layer = 'Replacement Tree') links to the correct feature class.

Is there a bug or a disconnect somewhere when using a map document as a template and resourcing the data from the Table of Contents to the python interactive window so that layers (which show up on the data frame, too!) aren't connected as if they were added in manually?

Also posted to GeoNet..

Update From GeoNet

Following a suggestion to check the dataSource of each of the layers in my list, I tried a little different code.

I ended up with this copied from the interactive window.. Formatting is a little weird from the copy/paste.

>>> for lyr in arcpy.mapping.ListLayers(mxd, "*Replace*"):
...     print lyr.dataSource
...
C:\***\*****\F_Geodatabase\*****.gdb\ReplacementTree
>>> for lyr in arcpy.mapping.ListLayers(mxd, "*Replace*"):
...     print lyr.dataSource
...     print arcpy.Exists(lyr)
...
C:\***\*****\F_Geodatabase\*****.gdb\ReplacementTree
False
>>>

Now from there, I hit add data and added in the same feature class, from the same geodatabase and re-ran the above script to get:

>>> for lyr in arcpy.mapping.ListLayers(mxd, "*Replace*"):
...     print lyr.dataSource
...     print arcpy.Exists(lyr)
...
C:\***\*****\F_Geodatabase\*****.gdb\ReplacementTree
False
C:\***\*****\F_Geodatabase\*****.gdb\ReplacementTree <--Put in here by hitting add data button (same file)
True
>>>

So, I don't exactly get why two TOC items, both sourced to the same feature class, shows as one existing and one not.  The only difference is that one (that doesn't exist) was re-sourced via a python script, while the other (that exists) was added manually.


texts: MapServer rendering of points takes too long

texts:
I have MapServer configured with a Mapfile and things are mostly fine except for the fact that my shapefile has many many points in it.

This means MapServer takes a while to render the image at far zoom levels.
I tried MAXFEATURES to lessen the processing but it does not distribute the points evenly across the areas where the points should be.

Is there some other Mapfile directive that anyone can enlighten me with?


texts: PyQgis does not load a vector layer to form

texts:
I will want create "Hello World" project at PyQgis. I will want to show map on    form.ui from .shp-file, but it does not  download and show.
But it is not loaded and does not appear. It`s value variable "layer" after load .shp

    file = QFileDialog.getOpenFileName(self,
                   "Open Shapefile", ".", "Shapefiles (*.shp)")
    fileInfo = QFileInfo(file)

    # Add the layer
    layer = QgsVectorLayer(file, fileInfo.fileName(), "ogr")

    if not layer.isValid():
      return

    # Add layer to the registry
    QgsMapLayerRegistry.instance().addMapLayer(layer);

    # Set extent to the extent of our layer
    self.canvas.setExtent(layer.extent())

    # Set up the map canvas layer set
    cl = QgsMapCanvasLayer(layer)
    layers = [cl]
    self.canvas.setLayerSet(layers)

layer.isValid() = false
map files have extensions:
.cpg
.dbf
.prj
.qix
.shp
.shx

I tried to change. Nothing changed


texts:
I am working with ArcGIS Android..Offline Geodatabase. I have Point as a graphics and Polygon layer as a feature layer.I want to find out Point in Polygon?

I have used the following code :

query = new QueryParameters();
               query.setOutFields(new String[] { "*" });
               query.setSpatialRelationship(SpatialRelationship.INTERSECTS);
               query.setGeometry(mapPoint);
               query.setInSpatialReference(mMapView.getSpatialReference());
               QueryTask qTask = new QueryTask(**featureLayerEight**);

I am getting following error:

  Error:(276, 34) error: no suitable constructor found for
  QueryTask(FeatureLayer) constructor
  QueryTask.QueryTask(String,UserCredentials) is not applicable (actual
  and formal argument lists differ in length) constructor
  QueryTask.QueryTask(String) is not applicable

(actual argument FeatureLayer cannot be converted to String by method invocation conversion)


texts: Unable to change polyline segment to circular arc

texts:
I have a poly line connected by a variety of vertices. I think the polyline is a straight line rather than a proper arc between each vertex. I was trying to follow the instructions here but the cursor never came up to select the individual segment. I even split the polyline at the vertex prior to the line segment I wanted to alter, but I find no option to "Change Segment".

Is there another way to do this? I have ArcGIS10.3 (specifically the ArcMap application).


texts: Trouble returning wms tile in Cesium, from geoserver, vnd.ogc.se_xml

texts:
I am running Cesium in sandbox, and trying to access wms layers from geoserver I am running locally on a virtual machine but having trouble returning the wms data.

The request hits the tile image, but I think nothing is being returned, or there is some issue with processing the data on the way back when trying to be processed in Cesium.

ex: the request url: https://192.168.99.100/geoserver/geonode/wms?version=1.1.0&service=WMS&request=GetMap&styles=&format=image%2Fjpeg&layers=ne_50m_admin_0_countries&srs=EPSG%3A4326&bbox=-180%2C0%2C-90%2C90&width=256&height=256

reaches an image/tile:

the response headers:
Access-Control-Allow-Origin:*
Connection:keep-alive
Content-Type:application/vnd.ogc.se_xml;charset=UTF-8
Date:Mon, 10 Aug 2015 19:46:08 GMT
Server:nginx
Strict-Transport-Security:max-age=31536000; includeSubdomains;
Transfer-Encoding:chunked
X-Frame-Options:SAMEORIGIN
X-XSS-Protection:1; mode=block

how do i interpret the type: vnd.ogc.se_xml ?

the errors are: An error occurred in "WebMapServiceImageryProvider": Failed to obtain image tile X: N Y: N Level: N.
The relevant code is as follows:

var viewer = new Cesium.Viewer('cesiumContainer');

var provider = new Cesium.WebMapServiceImageryProvider({
url: 'url_to_wms_server',
layers : 'layer_name'
});
viewer.imageryLayers.addImageryProvider(provider);

any idea what's going wrong or what to change?

also, when i change the version to 1.3.0, the response type is text/xml but the url does not hit an image when i check.