texts: Problem with Leaflet Control Layers using Vector Layers

texts:
I'd like to add a Vector Layer and configure it in a Leaflet Map with Leaflet Control Layers.

The code is ...

          <!DOCTYPE html>

              <meta charset="utf-8" />
              <title>Example</title>
              <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css" />
              <link rel="stylesheet" href="http://netdna.bootstrapcdn.com/font-awesome/4.0.3/css/font-awesome.css" />

It doesn't work: my vector layer doesn't appear on the map, but if I remove these rows...

                var baseMaps = {
                    "Basemap": baseLayer
                };

                var overlayMaps = {
                    "Mapillary Sequences": mapillarySequences
                };

                L.control.layers(baseMaps, overlayMaps).addTo(map);

... all works fine.

Any suggestion, example, jsfiddle, alternatives?

Cesare


texts: Merge polygons separated by road spaces while keeping the boundaries of separations

texts:
I made a Merge of cadastral parcels to have localities.
The problem is that I have the polygons of 6 localities but with voids/spaces in the middle that represents roads.

I would like to merge/union these 6 polygons of localities in order to have polygons by localities but with good defined boundaries glued to each other, without separations by roads.


texts: How to export csv/excel in QGIS3.0 model processor?

texts: In the model processor for QGIS 3.0, how do I have the attribute table export as an excel? I have five different outputs/excels I need and would like to automate the process instead of needing to gobto each layer individually.


texts: How severe is a "10001 layer's data source has a different projection value than the data frame projection" in ArcMap?

texts:
I am currently practicing using the generate near tool with some older point shape files and would like to see if I can use it in distance analysis between the points and other shape files like lines and polygons that represent features of interest. When I uploaded the point shape files and the features of interest, all the features are projected close to each other. Since the points and features of interest are all near by each other in the same region in the real world, I assumed that everything was well. However, when I analyze the map, I get a 10001 error of High severity.I am not sure how grave this issue is in my case. From what I read, it may have to do with drawing performance of the layers, but I want to make sure I am not making a mistake.

The data frame properties says the coordinate system is the following:

Projection: Transverse_Mercator
False_Easting: 500000.0
False_Northing: 0.0
Central_Meridian: -99.0
Scale_Factor: 0.9996
Latitude_Of_Origin: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_North_American_1983
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_North_American_1983
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

The features all have the same projection, one of the features of interest layer properties says the following:

Data Type:  File Geodatabase Feature Class
Database:   C:\Ento Files\BollWeevilEradication\BollWeevilEradication.gdb
Feature Class:  T2011_LRGV
Feature Type:   Simple
Geometry Type:  Polygon
Coordinates have Z values:  No
Coordinates have measures:  No

Geographic Coordinate System:   GCS_North_American_1983
Datum:  D_North_American_1983
Prime Meridian:     Greenwich
Angular Unit:   Degree

I can see that they both have the same Geographic coordinate system but the 10001 warning has me worried that my distance analysis might end up being wrong.


texts: What sources of free multispectral satellite imagery are there in Europe other than Sentinel 2?

texts: I am looking for free multispectral satellite imagery for north Germany prior to 2015. Does anyone have a list of sources?


texts: Statistical downscaling of raster with sum function uisng R

texts: Usually raster downscaling methods exploit interpolation techniques. However, I am not sure how to proceed in R if a raster needs to be downscaled (=i.e. its resolution increased, e.g. from 1 km to 100 m) such that the sum of the underlying pixels yields the value of the original lower resolution pixel. Of course I would like to to this in a way that not all pixels underlying the coarser resolution pixel have the same value, but based on the distribution of pixels on another raster layer (i.e. with a statistical regression approach).


texts:
I am using GDAL to convert a raster dataset into a numpy.array and then using numpy to buffer the data in the raster.  For example, I need to know all the regions within 100ft of a school point feature and I need it represented as a raster.

So, I convert my points to a binary raster (school points are 1 and all other pixels are) using GDALRasterize in Python.  Then, I open the raster dataset and use GDAL.ReadAsArray to create a numpy.array for the raster.  Using scipy.ndimage.morphology.binary_dilation, I buffer all the non-zero pixels.

The issue is: some of my input datasets can be rather large and so I am running out of memory.  The first thing that came to mind was to tile my input data and process individual tiles.  I'm not sure how to handle any edge effects if I were to tile.  What if a buffer spans the edge of a tile?  The next tile wouldn't know to include those buffered cells.

Does anyone have any suggestions for how to handle these kinds of edge effects in numpy/GDAL?


texts:
I've made a map with openlayers and have the vector features clustered with a clustering strategy.

Now I want my users to be able to toggle this behaviour with a control similar to the layerswitcher.

I know the cluster strategy has activate() and deactivate() methods. But how can I combine those with a custom control?
here's a jsFiddle with a working clustering example http://jsfiddle.net/BzcAD/5/

many thanks

EDIT:
I tried the sollution proposed below, but somehow adding more strategies to the vector layer forces OL to look for a protocol. I've got no idea how to fix that or to work around. Still looking for a better suggestion. Take a look at: http://jsfiddle.net/BzcAD/6/

If you look at http://www.openlayers.org/dev/examples/strategy-cluster-extended.html there they just throw away everthing and build the complete layer again. I don't like that.

The only work around I can think of is setting the clustering threshold one above the number of features on the layer. I'de like to have a layerswitcher style folding menu calles 'options' with a selectbox toggling this clustering behaviour.


texts: qgis python: select on screen and change (pre-set) properties accordingly

texts:
I have a map with different zones where I have drawn isochrones -where I will arrive within 15 min, within 30 min and so on. The picture attached shows the 15-120 min isochrones.

This is set as a matrix table, where each column represents a zone (D2000, D2001, and so on) and every row is the time between the zone selected and the rest of the zones.

The plot currently shows travelling distances from D2000, and if I want to select a different origin, I go to properties > style > column and I chose a new one (i.e.: "D2001"). Then, I will get the isochrones for this new origin "D2001" as the properties (eg: Colours) will remain the same (so I don't need to change anything just ask qgis to refresh the layer to get the new values).

The problem is that I need to present something more interactive. The idea is to avoid the "properties > style > column" but just to select a value from the screen (using "select features by freehand" -only one feature each time-) and automatically run the "properties > style > column" process to show the selected zone as the origin for the new isochrones.

I believe it could be done in python, but unfortunately I don't possess that knowledge yet (it is my goal for 2016!) but I need this map done the sooner the better.

Basically, it would be something like this:

layer selection
select features by freehand
selection of (just) one zone ("D2000")
go to properties > style > column
change default column ("ZONE") to select "D2000"
apply
allow any type of reset to set the default column "ZONE" in order to start the process again.


texts: Is there a name for a situation when the place is clearly or unclearly named?

texts:
This is not specific to any system. Suppose there's a large red rock on a terrain which can be easily found and clearly seen and there's an artificial object such as a mall near that large red rock.

That object could have a name that matches its location such as "Large Red Rock Mall" and so it would be intuitive to look for that mall near that rock or it could have some mismatching name that would make the search much harder.

I'm pretty sure this arises in GIS systems and causes people to place some kinds of metainformation next to such badly named objects.

Is there a widely used term for these two situations that would mean something like "how well the object name matches its location"?


texts:
Is there a way to add locations and info windows to my own map in .tif format? I'm using MySQL and PHP. Can I use Google Maps API pins but not their map? I thought I can put my georeferenced map (WGS84) and add Google Maps API pins on it. Is that possible and if it is not with Google Maps API, maybe you can recommend me some other way?

Thank you.


texts: Locating and using Full Motion Video (FMV) toolbox for ArcGIS Desktop?

texts:
I've just installed version 1.2 of the Full Motion Video Add-In, along with the new geoprocessing toolbox. After running the installation for the tools, it creates a folder on the C drive but does not contain a toolbox. There is no link from the FMV Video Manager pane to the toolbox, nor does it appear in the system toolboxes. I tried to "Add Toolbox" and looked in all the normal places (Documents, Users, C drive etc) and came up empty.

Where is the toolbox located?

I would like to test out the new tools as mentioned in the blog post. I'm running ArcGIS 10.2.1.


texts: How can one transform XY and Z coordinates to another coordinate system?

texts: How can we convert X ,Y and Z values from one coordinate system to another coordinate system using node js or any other similar packages or equations without directly using the convertor tools available online?


texts: Seeking GIS data for fast food industry?

texts: I am working a small project related to fast food industry. I am trying find GIS data for fast food chains like subway, MC Donald's, Arby's, etc in state of Ohio. Basically, I am looking for point location of each facility (for each fast food chain). I know that I could digitize them but that would be painful. Is there a faster way of getting this information? Excel files with the address details will also help.


texts: How to install PostGIS shp2pgsql-gui on openSUSE 13.2?

texts:
I installed PostgreSQL 9.4.4 together with pgAdmin III (1.20.00) and the PostGIS extension (2.1.7, via the Stack Builder) from http://www.enterprisedb.com/products-services-training/pgdownload on my openSUSE 13.2 (32 bit) machine.

Unfortunately, the GUI of the Shapefile-Loader (shp2pgsql-gui) has not been installed. After endless Google queries and reading similar questions here, I gave up and hope that I can get a helping hand here.


texts: How to categorize features based on hstore values in QGIS?

texts:
I have a table with the columns of id,tags(hstore data type) and geometry.
I want to categorize the table based on selecting the specific key form tags:
for example:

select * FROM table_name WHERE tags @> 'public_transport=>stop_position'


texts: Postgis load raster to DB running on AWS RDS with raster2pgsql

texts:
I want to perform very common task:

I have Postgresql + Postgis running on Amazon Web Service cloud, using RDS service.
I have there some vector layers and can execute SQL (Postgis) commands via connection using Pgadmin3 SQL console. So I know for sure that Postgis is installed OK.
Also I can execute psql console with SELECT PostGIS_Version(); running OK.

The next is to execute something like raster2pgsql -G; in order to test that the raster2pgsql tool installed. The problem is it doesn't.

Basically When one using Potsgresql on AWS RDS the Postgis is pre-installed. And if so, tools like "raster2pgsql" should be automatically also be installed.

How now I can install this tool on such remote DB? I have access to psql and that's it.
My task is to load raster data (topo) in the way I can leveraged Postgis spatial queries. What are my other alternatives to load into Postgis DB DTED files?


texts:
I unregistered moveend event.

map.unByKey(keyMoveend);
console.log('keyMoveend Unregistered');

What if I need reRegister that event later. What do I have to do?
or
How can I suspend moveend event temporarily later to use?


texts:
I wanted to join a CSV file to a shapefile.

When I identified a field to facilitate the join, e.g.,

map <- spChFIDs(map, as.character(map$ID))

it returned

  Error in spChFIDs(SP, x) : lengths differ

Can anyone advise?


texts: json issue: GetFeatureInfo with raster images

texts:
I have succesfully typed a GetFeatureInfo request with the following structure directly into my browser:

http://192.168.1.20/geoserver/workspace/ows?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetFeatureInfo&CRS=CRS:84&BBOX=39.32295698174142,-8.4236378760969881,39.363077494358599,-8.3933494101086144&WIDTH=1114&HEIGHT=841&LAYERS=layer&QUERY_LAYERS=layer&FORMAT=image/png&INFO_FORMAT=application/json&FEATURE_COUNT=50&I=396&J=579

It returns me onto the screen:

{"type":"FeatureCollection","features":[{"type":"Feature","id":"","geometry":null,"properties":{"GRAY_INDEX":162}}]}

This image was a 8-bit jpeg.

Now, when I try to do the same over a 16-bit raster image, it does not return any value:

{"type":"FeatureCollection","features":[]}

I was wondering whether the raster type is a key-fact0r on this issue.

Do I have to run the query in a different way?


texts: Mapping relation nodes or way nodes to lat/lon coordinates?

texts:
The goal is to construct a data structure that takes the id of a node/way/relation and returns the pair of lat/lon or the collection of lat/lon coordinates assigned to it.
So far, I have downloaded the osm file of my city, gotten rid of superfluous nodes like those pointing to amenities, and then used osmconvert with the all-to-nodes feature to construct a csv file that contains the columns @Id@name@lat@lon. This approach does work to construct a look up system for id(of node)->lat/lon. But not for id(of relation/way)--> List of lat/lons.
Does anybody know if there is a osmconvert feature or another way to retrieve this kind of list or csv file.


texts: GeoJSON to shapefile conversion with GeoTools

texts: Can someone please provide me an example of how to convert GeoJSON to shapefile using GeoTools?


texts: Why are xy coordinates displayed in WGS84 project instead of lat/lon?

texts:
I have a project in QGIS comprised of shapefiles and raster maps. The origin CRSs are mostly on ETRS89 / Portugal TM06. I am trying to convert the project to WGS84 to have lat/lon instead of metric values on the map coordinate labels.

However, I cannot do it. I have enabled on-the-fly transformation after choosing WGS84 EPSG:4326 but still cannot get lat/lon coordinates in the coordinates window in the bottom of QGIS screen or on map coordinates in Print composer. All my layers are on WGS84 (EPSG:4326).


texts: QGIS DEM (Terrain models) Mode Slope creates wrong data

texts:
I have a DEM raster file downloaded from the USGS site. The file can be downloaded here.

I am trying to create a slope raster in percent with the DEM tool in QGIS. The GDAL command is the following:

gdaldem slope /Users/Downloads/n41w107/imgn41w107_13.img /Users/Downloads/slope.tif -p -s 1.0 -b 1 -of GTiff

The tool runs fine, but the created file slope.tif has strange data. The statistics for the band are the following:

STATISTICS_MAXIMUM=29834248
STATISTICS_MEAN=1845376.751065
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=1658350.0283881

So they are obviously way to big for percent. What am I doing wrong?

Update: I calculated the Slope in degrees and got the following statistics:

STATISTICS_MAXIMUM=89.999809265137
STATISTICS_MEAN=89.583139343302
STATISTICS_MINIMUM=0
STATISTICS_STDDEV=5.9592219763749

A mean of 89 degrees seems really high to me. Is that reasonable?


texts:
I want to perform snapping in all the lines of a shapefile that don't connect to the nearest other lines. It would be okay to have a tolerance of 0.5 meters. So between 0.5 it should connect to the near lines. I have done this:

import geopandas as gpd
from shapely.ops import snap
s=gpd.read_file(r'C:\Users\user\Desktop\ag_eus')
result=snap(s,s,0.5)

It gives:

AttributeError: 'GeoDataFrame' object has no attribute '_geom'

Some more info about the shapefile:

s.columns

Index(['Layer', 'SubClasses', 'ExtendedEn', 'Linetype', 'EntityHand', 'Text',
       'geometry'],
      dtype='object')

UPDATE from the answer that provided:

for index, row in s.iterrows():
    tmp_gdf = s.copy()
    tmp_gdf['distance'] = tmp_gdf.distance(row['geometry'])
    closest_geom = list(tmp_gdf.sort_values('distance')['geometry'])[1]
    # I took 1 because index 0 would be the row itself
    snapped_geom = snap(row['geometry'], closest_geom, 0.5)
    s.set_value(index, 'geometry', snapped_geom)

    #new addition
    o = r'C:\Users\user\Desktop\ag_eus'
    s.to_file(o+ '\\new')

and the file is exactly the same with the first?


texts:
I'm new with editing so I need help. I have several lines (one feature class) which share the same location (cables) but I want to represent them as lines that are next to each other because they have different properties. How exactly do I do that? Using representation rules? What if we're talking about two different feature class?

I use ArcEditor licence, ArcGis 10.

Thank you


texts: Problem with projection of layers in QGIS

texts:
I'm having a problem with projections in QGIS that I hope someone here can help me with. I have been searching and found other people having the same problem, but have not been able to solve mine based on those topics.

I have a shapefile, converted from a .SOSI file, and a .csv file made by me. Both layers are displayed, but they are two different places like shown in this picture:

I have also two pictures of the some general info of the shapefile and .csv file. Hope they can give some info.

Hope you can help me out here, can't figure this out.

On the fly reprojection is enabled, I have tried to choose both same CRS and different CRS, nothing is working.

If you need any more info to solve this, just ask me and I will provide it.


texts: How to add a field to a shapefile using R

texts:
This question is trying to make sense of a post that already exists that looks good, but is failing in practice.

The second ranked answer uses the foreign package to read a shapefile attribute table as a .dbf, then allows @mdsummer to add a field to the table as an additional attribute.

I am using Raster and Shapefiles packages (shapefiles borrows the foreign read and write dbf).

library(raster); library(shapefiles)
shp<-shapefile(ZoneShape);

shp

>class       : SpatialPolygonsDataFrame
nfeatures   : 213
extent      : 577829.3, 582592.8, 3837297, 3839300  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=15 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0
nvariables  : 1
names       : OBJECTID
min values  :    24571
max values  :    28496

Now I'll read my attribute table:

shp.AT<-read.dbf(gsub(".shp", ".dbf", ZoneShape))

class(shp.AT)
>[1] "list"

shp.AT is a list that contains an index produced by R and the Shapefile's OBJECTID. Their are 213 records, but the length of the list is 2.

I'm not a great programmer, but I know I can't add a field to a list.

shp.AT$GID<-1:nrow(ZoneShape)
>[3] ERROR: nrow(ZoneShape)

nrow(ZoneShape)
>nrow(shp.AT)

# "nothing is returned!"

@mdsummer includes an "as.is = TRUE" in the read.dbf statement. Not sure what that might have done in the past, but R returns an error when this is included.

If I modify my new field slightly:

shp.AT$GID<- 1:length(ZoneShape)
shp.AT

This time, no error.
My list of ID's is again reported along with my new GID field, but:

>$GID
[1] 1

If this worked as expected, I would overwrite the old DBF with the write.dbf command.

So, a couple questions really:
Has the source changed, or is @mdsummer's post incorrect?
or, am I totally missing something?

and (most importantly) how can I add fields to my shapefile attributes.

I also have a cross-posting on Overflow


texts: Could any one give me some good flow in diagram of how lat, long information from a csv file is changed into a point in cartoDB. I am assuming it is being done with the help of SQL queries, but still is there any architectural diagram of the data flow in cartoDB?


texts: Dissolve features in ArcGIS using multiple attributes

texts:
This may sound simple, apologies if I missed some earlier discussion on same subject.

I have a big set of road data with road numbers and also zone numbers . The road names can repeat across the zones, ie Road XYZ will be there in Zone A and Zone B.
I would like to dissolve the road line features using "Dissolve" command, but see that roads in tow zones are note merged together.

Can the dissolve command be tweaked to achieve this ?


texts: View streets with only 2019 Google street view

texts:
Is there a way to view the streets that have a specific year street view, or is there a way to find the routes the Google Street Car travelled in each year?
eg. A layer that can show streets that have street view imagery from 2019


texts: Error coercing raster with one non-NA value to SpatialPixelsDataFrame

texts:
I'm getting a strange error when trying to convert a raster with only 1 non-NA pixel to SpatialPixelsDataFrame.

r0 <- raster( matrix( 1, 1,1))
as(r0, 'SpatialPixelsDataFrame')

  Error in (function (classes, fdef, mtable)  :
    unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric"’

And it doesn't seem to depend on the size of the raster

r1 <- raster( matrix( c(1, NA), 1,2))
as(r1, 'SpatialPixelsDataFrame')

  Error in (function (classes, fdef, mtable)  :
    unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric"’

r2 <- raster( matrix( c(1, rep(NA,99)), 10,10))
as(r1, 'SpatialPixelsDataFrame')

  Error in (function (classes, fdef, mtable)  :
    unable to find an inherited method for function ‘coordinates’ for signature ‘"numeric"’

Yet when I include 2 or more non-NA pixels ...

r3 <- raster( matrix( c(1,1, rep(NA,98) ), 10, 10))

as(r3, 'SpatialPixelsDataFrame')

Everything works fine.

Object of class SpatialPixelsDataFrame
Object of class SpatialPixels
Grid topology:
   cellcentre.offset cellsize cells.dim
s1              0.05      0.1        10
s2              0.05      0.1        10
SpatialPoints:
        x    y
[1,] 0.05 0.95
[2,] 0.05 0.85
Coordinate Reference System (CRS) arguments: NA
...blah blah blah

I realize this is a boundary case, but any idea what is going on here?

I'm getting this error only occasionally as part of a larger automated process, so it would be great if there was a way to convert my 1-pixel raster to a SPDF without an elaborate hack.


texts: Modifying several hundred geojson files by adding a feature / attribute

texts:
I have just downloaded several hundred identically formatted geojson files. The problem is none of them have any attributes in their file structure just coordinates. I want to add the following lines at the top of every file using either python2 or python3:

   content = {'"type": "FeatureCollection",'
      '"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },'
             '"features": ['
             '{ "type": "Feature", "properties": { "id": "%s" }, "geometry":' %i
             }

Where i is just the filename. I have tried a variety of things but nothing seems to get all the way there. There always seems to be an error with the string, the combination or the output. I would also like to add a ]} at the end but that is fairly straightforward append. Any help would be awesome.

import json
with open('TestLayer.geojson', 'w') as geojson_file:
data = json.load(geojson_file)
# data is now a dict,
content = {'"type": "FeatureCollection",'
           '"crs": { "type": "name", "properties": { "name":  "urn:ogc:def:crs:OGC:1.3:CRS84" } },'
           '"features": ['
           '{ "type": "Feature", "properties": { "id": "%s" }, "geometry":'
           }
data = content+data
json.dump(data, geojson_file)


texts: Does leaflet support wms tile layers in applications/openlayers formats?

texts: I am trying to create a tilelayer.wms overlay on top of a wfs base map and I wasn't sure whether leaflet allows the applications/openlayers type of format which I served up from geoserver or only images/pngs.


texts: Optional parameters in PointSymbolizer href

texts:
I am trying to build a dynamic url for a PointSymbolizer:
/api/icons?foo=${strUrlEncode("foo")}&bar=${strUrlEncode("bar")}

The problem I am running into is that the strUrlEncode will return null if the "foo" property is null.  When that happens I get the following exception:
DEBUG [function.string] - Failed to concatenate string in Concatenate function
java.lang.NullPointerException
    at java.net.URLEncoder.encode(URLEncoder.java:204)
    at org.geotools.filter.function.string.URLEncodeFunction.evaluate(URLEncodeFunction.java:87)
    at org.geotools.filter.expression.ExpressionAbstract.evaluate(ExpressionAbstract.java:45)
    at org.geotools.filter.function.string.ConcatenateFunction.evaluate(ConcatenateFunction.java:66)
    at org.geotools.filter.expression.ExpressionAbstract.evaluate(ExpressionAbstract.java:45)

I imagine that is happening because the result of my first strUrlEncode (which is null) is being concatenated to the second parameter.
Is there any way I can provide a default non-null value in a case where foo is null to avoid the NPE?


texts:
Is it possible to disable OsGeo Live 7.0 Welcome window? It annoys me a bit appearing each time I Run the system.

Solved by:
menu->settings->settings manager->session and startup->"Application Autostart"->uncheck the message appearance


texts: MapServer computes wrong map scales -- off by factor 0.77. Or is it just me?

texts:
I've noticed in a current project that MapServer seems to get the  map scale wrong, or then I'm wrong, but don't see how.

My mapfile contains the following definitions:

RESOLUTION 72   # this is the default, but I set it explicitly, just to make sure

PROJECTION
    "init=epsg:21781"  # = CH1903 LV 03, which uses meters as coordinate units
END

LAYER
    NAME "one"
    MAXSCALEDENOM 7500
    ...
END

LAYER
    NAME "two"
    MINSCALEDENOM 7502
    ...
END

If I only display these two layers, I should definitely be able to get an empty map at scale 1:7,501.

I have written a simple map viewer application to verify this. My viewer also assumes a screen resolution of 72 DPI, just like MapServer. However, it turns out that both layers are invisible around map scale 1:5,803. I would have to adjust my viewer to 93.062 DPI to get the empty map at 1:7,501. Why is this?

Just for reference, this is how I calculate the map scale in my viewer:

double imageWidthInPixels = ...;   // passed to mapserv.exe (via 'map_size' arg.)
double extentWidthInMeters = ...;  // = extent's xmax-xmin (via 'mapext' arg.)

// constants:
const double dpi = 72.0;
const double inchesPerPixel = 1.0 / dpi;
const double inchesPerMeter = 1.0 / 0.0254;

double inchesOnScreen = inchesPerPixel * imageWidthInPixels;
double inchesInReality = inchesPerMeter * extentWidthInMeters;

double scaleDenominator = 1 / (inchesOnScreen / inchesInReality);

How could this go wrong?

P.S.: I am quite certain that I have the right extentWidthInMeters since I can easily derive it from the extent. (I'm using the Swiss projected coordinate system CH1903 LV03, which already uses meters as coordinate units.)

UPDATE: It seems as if the inaccuracy was partly caused by differing aspect ratios of the map image control in my viewer and the extent sent to MapServer. I've now made sure that only map extents are passed to MapServer that have the same aspect ratio as my map image control. But some small error still remains: The map is now empty around scale 1:7'489 instead of 1:7'501.


texts: OpenStreet Map - Download power data for ArcGIS

texts:
I want to find power data of OSM. After search I found this source: https://download.osmdata.xyz/

It gives the power data of all the world (6GB file, hard to edit). Is there any other source to download only power data per country (Greece specifically) or at least per continent?


texts: QGIS displaying polygon that is crossing the 180 meridian

texts:
I have the following table

create table schema.table
(
    id bigint,
    geography geography
);

with the following geography

INSERT INTO schema.table (id, geography) VALUES (1, '0103000020E61000000100000006000000820DAB1F7BC86440D91A1EB605D733C024C27A3479CB61409763729D583728C0755E7C3479CB614091CC699D583728401184AA1F7BC864404AA108B605D73340808702B3EB7166C095D626E80B2E31BE820DAB1F7BC86440D91A1EB605D733C0');

This describes a simple hexagon that is crossing the 180th meridian, but QGIS displays is as follows.

Permise: I cannot change geography to geometry. And I have thousands of polygons that have the same issue, manual fixing is not an option.

I have looked at ST_ShiftLongitude and ST_WrapX as I have worked with these functions before, but they only apply to geometry data.

Question: What settings do I have to change in QGIS to get the right visualization?

Note: In PostGIS, e.g. ST_Intersects and ST_DWithin work just fine.

QGIS version
3.4.3-Madeira

Update 1:

I was able to get a visually acceptable solution with

SELECT id, st_wrapx(st_shiftlongitude(geography::geometry), 180, -360) as geom FROM schema.table;

Essentially the tip that is crossing the 180th is cut off and placed on the 'left' side of the map. However, this only works for polygons that are actually crossing the 180th.

Also, I made a mistake, it is a pentagon, not a hexagon in this example.

ST_IsValid returns true and I could not see a self-intersection.

POLYGON((
166.265029749 -19.839930899,
142.358545532 -12.108097954,
142.358545535 12.108097953,
166.265029748 19.839930894,
-179.560021882 -0.000000004,
166.265029749 -19.839930899))


texts: Convert from spherical angles to longitude and latitude

texts:
I have some images extracted from a radar from Austria and a formula for the projections (the exercice says The projection formulas map pixel locations to geographical coordinates), that transforms pixels into (phi, lambda) which the exercise states that are the spherical angles.

PRODUCT-PH-FORMS   =     "(1.570796326795-2*atan(0.404026225835*pow((1/sin(acotan(5899179.2/\
((5676531.8*pow((tan((1.570796326795-(v1*0.017453292520))/2)/\
0.373884679485),0.737361597616)*sin(0.012869387656*(v0-13.333333333))-\
v3*v4)+COL*v4)-((6380000.0*(0.924636243305-0.889738536848*pow((tan((\
1.570796326795-(v1*0.017453292520))/2)/0.373884679485),0.737361597616)*\
cos(0.012869387656*(v0-13.333333333)))+v2*v5)-ROW*v5)/((5676531.8*\
pow((tan((1.570796326795-(v1*0.017453292520))/2)/0.373884679485),\
0.737361597616)*sin(0.012869387656*(v0-13.333333333))-v3*v4)+COL*v4))))*\
((5676531.8*pow((tan((1.570796326795-(v1*0.017453292520))/2)/\
0.373884679485),0.737361597616)*sin(0.012869387656*(v0-13.333333333))-\
v3*v4)+COL*v4)/6380000.0*0.737361597616/sin(0.767944870878),1/\
0.737361597616)))"

PRODUCT-LA-FORMS   = f(pixel row, pixel col, v0, v1, v2, v3, v4, v5)

But this gives very small values which cannot be (latitude, longitude) so I am looking for the transformation between spherical angles (phi,lambda) to (latitude,longitude)

For the result I need to provide (latitude, longitude) so how this conversion could be made?
The exercise provides me as well with some data I do not understand, which may be necessary:

v0=8.194 (reference lon)
v1=50.437 (reference lat)
v2 = 0
v3 = 0
v4=1018.18 (meter / pixel)
v5=1018.18 (meter / pixel)

This variables are used for the transformation formulas but may be useful as well to obtain (latitude, longitude)

EDIT

I have tried to directly consider the values I obtain from the transformation as (latitude, longitude) after *180/pi but the coordinates I get are not in Austria


texts: Polygon area from Google Earth so different to QGIS & Shapely/PyProj?


texts: Converting DMS coordinates to decimal in R

texts:
I have a vector of DMS coordinates similar in format to the following:

coords <- c("43°46′50″N 79°24′53″W", "43°46′06″N 79°24′46″W")

and I would like to convert them into two numeric decimal vectors for latitude and longitude to eventually mark them in leaflet.

coords_lat <- c(43.780556, 43.768333)
coords_lon <- c(-79.414722, -79.412778)

My question is: Is there a function available out there in R that does this conversion automatically, or would I need to rely on regex and some outside tool to convert them myself?


texts: Automating subtracting two layers in QGIS?

texts:
Using QGIS, I would like to calculate the difference in depth between several geological horizons. I have one horizon that covers a large area (horizon1) and several smaller horizons (horizon2, horizon3, ...) within that area. All of these horizons have X,Y,Z values.

So basically I'd like to make new horizons which show the difference between the Z-values between horizon1 and the other smaller horizons. I have been able to do this in other (seismic) software manually, but I have so many horizons that I would like to be able to automate this. I can convert these horizons into shapefiles or rasters if necessary. I also have bounding polygons for the smaller horizons.

Edit: An example in this figure: the green dots (200x200 spacing) show horizon1, the dark red dots show horizon2 (50x50 spacing). After interpolating and using the raster calculator (horizon1-horizon2) I get this funny looking raster in black and white on the figure. Not sure why this isn't working. The black lines show bounding polygons of the smaller horizons (horizon2, horizon2, etc).


texts: How can I query and categorize a shapefile using MS Access?

texts: Is there a way to categorize locations by lat-long in MS Access using caparisons to shapefiles?


texts: Generating urban blocks by road centerlines using PyQGIS

texts: I work on cartographic generalization project and I need to generate urban blocks by road centerlines.


texts: Load feature class with tab delimited file that is missing Lat/Long

texts:
I am working on a Python script that receives a nightly tab delimited flat file. This flat file has Latitude and Longitude fields and the script loads the data using the gp tools Make XY Event Layer and Append management. The flat files sometime have missing Lat/Long values and when this occurs the script still successfully loads the feature class with the records even though they are missing Lat/Long.

The problem occurs if every record on the flat file is missing Lat/Long data. So it seems that as long as one record on the flat file has Lat/Long, the code works and loads the other records to the feature class as well.

Since receiving a flat file where all records have blank Lat/Long data is a realistic scenario, is there a different gp tool/approach that will allow me to load a full file of records to the existing feature class when all records are missing Lat/Long data?


texts: Which SRS is computational more effective to store continental scale data? (given some conditions)

texts:
For a Precision Agriculture application I need to make calculations in UTM projection. The data will be stored in Postgis and is made of soil samples results (point data) and field boundaries of many farms. But the farms are scattered through many utm zones (ex. a farm is in Z23 and the other is in Z22, so I can't store the data in UTM, it needs to be stored in a SRS without zones.

Since the calculations will be made in UTM the data will always need to be converted back and forth to the storage SRS.

So which of those inputs is cheaper in terms of computations of SRS conversion when the output SRS is UTM Wgs84:
-Geometry vs. Geography Postgis types
-Geometry projected (ex World Mercator) vs. Unprojected (latlon)

I think that it involves the internal Proj4 math but I couldn't understand it myself...


texts:
I used below codes to save IMap into byte[] (to store in my database), but I can't load it again, are there anyone know what's wrong ?

public byte[] SaveMap(AxPageLayoutControl pageLayOut)
{
    IGraphicsContainer graphicsContainer = pageLayOut.GraphicsContainer;
    graphicsContainer.Reset();
    IElement pElement = graphicsContainer.Next();
    while (pElement != null)
    {
           if (pElement is IMapFrame)
           {
               IMapDocument mapDocument = new MapDocumentClass();
               if (File.Exists(@"D:\test.mxd"))
                        File.Delete (@"D:\test.mxd");

               mapDocument.New(@"D:\test.mxd");
               mapDocument.ReplaceContents((IMxdContents)pageLayoutTo2_3.ActiveView.FocusMap);
               mapDocument.Save(true, true);

               System.IO.FileStream objFileStream = new System.IO.FileStream(@"D:\test.mxd", FileMode.Open);
               byte[] result = new byte[objFileStream.Length];
               objFileStream.Read(result, 0, (int)objFileStream.Length);

               return result;
            }
            pElement = graphicsContainer.Next();
      }
}

public void LoadMap(AxPageLayoutControl pageLayOut, byte[] result)
{
    IGraphicsContainer graphicsContainer = pageLayoutTo2_3.GraphicsContainer;
    graphicsContainer.Reset();
    IElement pElement = graphicsContainer.Next();
    while (pElement != null)
    {
           if (pElement is IMapFrame)
           {

               System.IO.FileStream oFileStream = new System.IO.FileStream(@"D:\test.mxd", FileMode.Append);
               oFileStream.Write(result, 0, result.Length);

               IMapControl2 map = new MapControlClass();
               map.LoadMxFile(@"D:\test.mxd", Type.Missing, Type.Missing );

               IMapFrame frame = new MapFrameClass();
               frame.Map = map.Map;
               pElement = (IElement)frame;

               pageLayOut.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);

               return;
            }
            pElement = graphicsContainer.Next();
      }
}


texts:
I have thousands of markers - I need to cache them and use show-hide tactics to only show like 100 markers at once but I have trouble with caching.

Adding all single markers/points to an array and if map is dragged or zoomed, then loop through the array to see it has already been loaded and then load all new markers that are not in that array is rather expensive.

How should I solve this caching issue?

How to load only markers/points/coordinates that haven't been loaded yet?

Could I benefit from grid system? A.k.a load all markers on that part of grid and cache grid - this would keep the array of already loaded areas/points much shorter.


texts: Why v.surf.idw returns only two value (one is nan)

texts: I am trying to perform a vector interpolation of a point shapefile with the intent of eventually creating a DEM. No matter what I try, I get only 2 values one of them is NAN, and nothing is working to fix it.


texts: Identifying area code for a particular US city

texts:
I'm trying to build a query like: "tell me what area codes (telephony area codes) overlap San Francisco" which should return 415 and 628.

Doing some google searches indicates that the US census data might be a good place to get this information from but I can't manage to make sense of the data on their website.

Any hints on how can I extract my desired information (City and overlapping area code information) from this data set?

Alternatively, are there any other, ready made, datasets with this information? Or APIs I could leverage?


texts: Create a histogram from a raster with no attribute table

texts: I converted an ASCII file to a raster. I had to use the float raster because it is continuous elevation data so there is no attribute table. I need to make a histogram in GIS of the elevations to figure out where to put transects, but am unable to do to the lack of attributes.


texts: Returned BBOX of a WMS with varying aspect ratio

texts:
WMS GetMap queries take a BBOX parameter for the map and a WIDTH and HEIGHT parameter for the image.

If the aspect ratio of the image (WIDTH/HEIGHT) doesn't match the aspect ratio of the BBOX (dx/dy) does the spec say it will distort the image into that WIDTH/HEIGHT aspect or does it say it will increase the area of the map returned in the image.

The server I'm playing with is doing the latter, but I see mentions on the web of the former behaviour.

In the case of increasing the area of the map, does the spec say that either the vertical or horizontal extent of the image will match the requested BBOX (depending on if its too tall or too wide)?

I suppose the safest thing is to always call with matching aspect ratios...

Here's some test grabs with three different aspect ratio WIDTH/HEIGHT settings and the same BBOX for your comparison:

wget -O test.png "http://www.osmgb.org.uk/OSM-GB/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=OSM-GB:__all__&SRS=
EPSG:27700&BBOX=532685,176236,542182,183918&FORMAT=image/png&WIDTH=400&HEIGHT=400"
wget -O test2.png "http://www.osmgb.org.uk/OSM-GB/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=OSM-GB:__all__&SRS
=EPSG:27700&BBOX=532685,176236,542182,183918&FORMAT=image/png&WIDTH=400&HEIGHT=800"
wget -O test3.png "http://www.osmgb.org.uk/OSM-GB/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&LAYERS=OSM-GB:__all__&SRS
=EPSG:27700&BBOX=532685,176236,542182,183918&FORMAT=image/png&WIDTH=800&HEIGHT=400"


texts: How to force login with user name and password at WMS service at a QGIS2WEB webmap?

texts:
Currently I am compiling a web map for an internal platform (15 users) via QGIS2WEB. There is also a WMS service (aerial photograph) integrated, which requires the input of user name and password, as soon as the map is clicked.

However, one problem is that the WMS service does not support "GetFeatureInfo" and therefore error messages appear when clicking with the mouse on the aerial image:

The same error message also appears when requesting information for objects that are also displayed on the map (e.g. construction projects), even if the layer with the aerial photos is deactivated:

Annoyingly, GetFeatureInfo must be activated in the QGIS2WEB dialog. If this is not the case, the query for user name and password is not started and the logon to the WMS server does not start. The aerial photo is therefore not displayed in the first place. Tested it several times already.

Is there any other way to "force" logging on to the WMS service than through a GetFeatureInfo request that causes error messages?


texts: Size based scale problems in point displacement styling

texts:
I am trying to style shapes for a continuous variable by size in point displacement. I used point displacement because I have multiple sampling events for sites. I styled my categorical data with colors. Then I ran into  problem when scaling the sizes of my points based on my continuous data (EP Chinook):

I think that my 2 data points for the site "Camping rio Llaima" are being averaged, and displayed as 1 dot size (see attribute table).

Since the values are 11 and 115 the dots at this site should be a very different size, considering the other point with a value of 58 looks similar to these points.


texts:
So far I can create a dropdown menu but in the shown example I could still choose 'coins' in "type" and '5$' in "value" which is wrong.
Are there rules I can add to make the options in "value" dependend of the choice in "type"?

And if so, is there also an option which automatically fills out "type" with the correct value when "value" is filled out?

I'd be thankful for a short example on how to set it up.


texts: Dissolve and then points along lines fails

texts:
I've a set of street data from OSM loaded into a QGIS layer.  A lot of the street names have multiple features (where OSM shows a road in a number of line segments).
I then use Geoprocessing Tools, Dissolve using name of feature and end up with each street having a unique name.  Each street name is now a single feature.
I now want to run the GDAL Points along lines to find the mid point of each street (and with the mid point lying on the street - as opposed to what you would get with Analysis Tools, Mean coordinates tool).
The GDAL Points along lines process fails (with a message saying look at the log, but I can't find any meaningful messages in the log).  I suspect that it is failing because a few of the street features are not connected lines - but are multiple lines that are not connected. The most obvious case is where the many features are un-named are joined together, but there are a small number of other features contain more than one line that is not connected. I can see several, but there are some I can't find.
Is there any way of finding the number of separated multiple lines per feature?  If there was, then I could exclude these features from the Points along lines and I suspect it would then work.
Fairly new to QGIS, and part of my problem in searching for a solution is understanding what the right terminology would be.
The log in the Points Along Lines is on Dropbox here
It ends

The following layers were not correctly
generated./Users/david/Dropbox/@Islington/Mapping/OSM
Clerkenwell Highway Points Dissolved.shpYou can check the
'Log Messages Panel' in QGIS main window to find more information
about the execution of the algorithm.

The main log only has

2020-06-29T17:22:51     INFO    ogr2ogr
"/Users/david/Dropbox/@Islington/Mapping/OSM Clerkenwell Highway
Points Dissolved.shp" "/Users/david/Dropbox/@Islington/Mapping/OSM
Clerkenwell Highways Dissolved.shp" -dialect sqlite -sql "SELECT
ST_Line_Interpolate_Point(geometry, 0.5) AS geometry,* FROM 'OSM
Clerkenwell Highways Dissolved'" -f "ESRI Shapefile"
2020-06-29T17:22:51     INFO    GDAL execution console output

Thanks for any help you can provide.


texts: Publish WMS with Geoserver and Openlayers from own server

texts:
I`m struggling with my Geoserver resp. OpenLayers.

Environment:
Ubuntu 1804, Tomcat8, Geoserver 2.13.2

I just want to publish vector and raster data from my own server. I can upload the files on my server and the Layers Preview is working correctly.

http://193.31.24.240:8080/geoserver/raster/wms?service=WMS&version=1.1.0&request=GetMap&layers=raster:500k_eng_fac&styles=&bbox=7.236923932667394,47.35615033439554,12.891923932667394,52.770837834395536&width=768&height=735&srs=EPSG:4326&format=application/openlayers#toggle

Also, the wfs-t example in the /www folder is working correctly.

http://193.31.24.240:8080/geoserver/www/wfs-t.html

But I can`t use any other map I tried the Simple Map Example from OpenLayers Example page, changed nothing but it is not working.

http://193.31.24.240:8080/geoserver/www/index.html

I tried about 15 different tutorials, the Projection system is clear EPSG:4326, checked the extent but I never get it running.

A normal HTML page is running so I assume my server is running and working.
http://193.31.24.240:8080/geoserver/www/index_work.html


texts: GeoServer - Caching with Response Headers Settings

texts:
I'm trying to find a method to cache one layer that its data is updated every 2 hours in the DB. From what I've found the best option is using Response Cache Headers timeout http://docs.geoserver.org/stable/en/user/webadmin/data/layers.html

What I've done is:

enable direct integration with GeoServer WMS.
add tiled=true parameter to WMS request
replace http://GEOSERVER_URL/wms with http://GEOSERVER_URL/gwc/service/wms
check "Response Cache Headers" with 7200 time specified as in http://docs.geoserver.org/stable/en/user/webadmin/data/layers.html

and when i look at the headers I see:

Cache-Control:max-age=7200, must-revalidate
Content-Disposition:inline; filename=geoserver-dispatch.image
Content-Type:image/png
Date:Sun, 04 May 2014 18:38:16 GMT
Expires:Sun, 04 May 2014 20:36:36 GMT
geowebcache-cache-result:HIT
geowebcache-crs:EPSG:900913
geowebcache-gridset:EPSG:900913
geowebcache-tile-bounds:12523442.712499999,2504688.5425000004,15028131.254999999,5009377.085000001
geowebcache-tile-index:[13, 9, 4]
geowebcache-tile-index:[13, 9, 4]
Last-Modified:Sun, 04 May 2014 07:06:05 GMT
Server:Apache-Coyote/1.1
Transfer-Encoding:chunked

I thought it worked but after 2 hours the layer didn't get updated! it's like that the layer is cached and frozen like forever. I think I'm missing something. any idea?


texts: Error on making a get request to a feature layer?

texts:
I am trying to query a feature layer from a script using the arcgis api for python:

from arcgis.features import FeatureLayer
lomr = FeatureLayer("https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer/1")

However, when I make a call to view properties, or actually perform a query:

lomr.properties.extent
lomr.query(where="DFIRM_ID = '08013C'",outSR=2876,outFields='*')

I get the following error (after quite a long traceback):

ConnectionResetError: [Errno 54] Connection reset by peer

I even get a similar error if I try pasting this query url through the requests library:

import requests
response = requests.get("https://hazards.fema.gov/gis/nfhl/rest/services/public/NFHL/MapServer/1/query?where=DFIRM_ID+%3D+%2708013C%27&text=&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=2876&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=false&datumTransformation=&parameterValues=&rangeValues=&f=html")

This seems like a problem on FEMA's side, as I'm guessing "peer" means the server of the client-server request. How do I avoid this error and actually get a successful get request?

[EDIT 1] - I thought this could potentially be a problem with using jupyter notebooks, but I tried this same process through python in the terminal and got the same error.

[EDIT 2] - I have also tried this through a virtual environment set up with pipenv with the same results. I've tried this on macOS and my Windows at work. Here are the versions of everything I'm working with:

pipenv - 2018.11.26
python - 3.7.4
jupyter notebook - 6.0.1
arcgis - 1.6.2


texts:
i have a database that has got tables with images and some of the tiles for the images is repeated and i want to eliminate the duplicates based on which imagery is oldest so as to remain with the newest images, i have run the following to get the duplicates from the database but do not know how to delete the row based on the age of the imagery

select "index",path from ngi_catalog_aug2012 where "index" in
(SELECT "index"
 FROM ngi_catalog_aug2012
GROUP BY "index"
HAVING ( COUNT("index") > 1 )
) order by "index"

this gives results like

index       path

2228BD_23    /data/ngi_old/images_30k_sid/2228/2228BD_23_498_469_12_0007_PAN_RECT.sid

2228BD_23   /data/ngi_old/images_30k_sid/2228BD/2228BD_23_498_469_11_0100_PAN_RECT.sid

and i would like to delete for instance the second row


texts: I am working on ArcGIS 10.0 and I want to SELECT (and not remove) all the spatial and address duplicates from the attribute table in a layer, either by code or in a table window or else with spatial operations.


texts: Using TauDem Tools in Python (like Arcpy)?

texts:
This may be a completely ignorant question so have some patience. But is it possible to use TauDem tools within an Python (Arcpy) script as per the other tools in ArcGIS? When you look at the help files of the TauDEM tools it gives a Syntax that looks similar to the syntax for the standard arcpy tools (this is what got me thinking).

In particular I am wanting to use the D-Infinity methodology for deriving flow direction and accummulation.

Feel free to shoot me down if I'm missing something obvious but I can't seem to find an answer anywhere.

I am using ArcGIS 10.


texts:
I am not a hydrologist but am providing GIS infrastructure to our hydrology department which actually doesn't use GIS very much yet. Therefore I am asking for some hydrology related insights.

Btw: We are using mainly ESRI ArcGIS but also other OpenSource GIS Software is possible to perform any of the described tasks.

Given some waterbodies alongside a large drainage system, I was asked if it is possible to aggregate (sum up) the water flow amount at every basin for the rivers upside the stream. That sounds like a pretty common task for me in hydrologial GIS Software. As I come mostly from transport networks, I was thinking of this as a directed graph, where water amounts act as edge weight, and aggregation is a simple dijkstra algorithm which sums up the weights. So I was thinking of scripting this task. But I wonder if there aren't ready to use tools for that.


texts: GDALWARP and OGR2OGR / CS2CS not agreeing on transformation?

texts:
I am trying to transform some data (both raster and vector) from WGS84 to a Croatian Coordinate system defined by:

+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs

This works great in CS2CS (Rel 4.8.0):

cs2cs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +to +proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs -r

And also great in OGR2OGR (version 1.9.2):

ogr2ogr -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -t_srs "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs" out.shp  in.shp

But in GDALWARP (version 1.9.2), the data comes out c.400m offset to the West (I think it's possibly a rotation error based upon it's appearance):

gdalwarp -s_srs "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" -t_srs "+proj=tmerc +lat_0=0 +lon_0=15 +k=0.9999 +x_0=5500000 +y_0=0 +ellps=bessel +towgs84=550.499,164.116,475.142,5.80967,2.07902,-11.62386,0.99999445824 +units=m +no_defs" in.tif out.tif

An example of this can be seen below, where I digitised a vector layer around existing raster data, transformed both using the above commands, and then loaded into QGIS:

I can't for the life of me figure out why this could be - it is worthy of note that I have tried this with multiple Raster datasets from a variety of sources, and they always seem to be offset to the West by the same amount.

I have tried upgrading to GDAL 1.10.1, and this has no effect upon the output.

Has anyone ever seen this before and have any idea why it could be?


texts:
Are there efficient ways of finding the union, intersection, and difference between two bounding boxes (box2df)?

I know that this can probably be done by first converting the bounding box (BOX2DF) to geometry type and then use ST_Union, ST_Intersection etc. But before doing this on a massive scale, I just want to check if that's naive and have any performance impacts with the conversion between bounding boxes and geometries.

Are there any more efficient alternatives like && is to ST_Intersects (over envelopes)?


texts: Set 'Load in QGIS' of all batch processing rows to 'no'?

texts:
I'm using 'Save Selected Features' algorithm from the Processing Toolbox.

I'm batch processing thousands of GeoJSON files to tab files but the same would be true for any conversion.  I've created a JSON input file to load in my thousands of records to convert, but QGIS opens each record created to the map view and therefore consumes all my system memory.

Is there a way to set all the rows to no without manually changing each, row by row?


texts: Plotting lat/long stored as text in Oracle DB using report builder/SSRS r2?

texts:
We use an Oracle 10g database.  However we write reports using SSRS or report builder, using the Oracle connection type to the datbase.  In oracle we have 2 text fields, one for lat and one for long.  (I am in the US and we use WGS84 coordinates).  Also, I only have SELECT rights on the oracle server.  I cannot manipulate the database.  I can create functions if that helps.

I need to plot these on a point layer in an SSRS map.  However, when I use the following code in my select statement, SSRS throws an error:

  Unsupported Oracle Data type USERDEFINED encountered

SDO_GEOMETRY('POINT(' || Longitude || ' ' || Latitude || ')',4326) as Coordinates

If I use the SQL code below I get an error that the FROM keyword is not found (makes sense since oracle doesn't use ::)

  geography::STPointFromText('POINT(' +
  CAST(longitude AS VARCHAR(20)) + ' ' +
  CAST(latitude AS VARCHAR(20)) + ')', 4326)

How can I use the lat and long stored in Oracle text fields as geo data types that SSRS can use to plot the point?
(NOT SQL Server... there is no SQL Server, only an SSRS server)


texts:
I have a rendering problem with the mapnik GEOS plugin.
I try to input a polygon WKT definition containing floats with around 16 digits.
Mapnik however does not render these. When I decrease the number of digits to 5 it does work. Does anybody have a solution?

mmap.background = mapnik.Color('steelblue')
s = mapnik.Style()
r = mapnik.Rule()
polygon_symbolizer = mapnik.PolygonSymbolizer(mapnik.Color('#f2eff9'))
r.symbols.append(polygon_symbolizer)
[ more style here ]
mmap.append_style('My Style',s)

wkt_geom = 'POLYGON ((5.12345 51.6, 6.8 52.6, 7.8 51.6, 5.12345 51.6))' #DOES WORK
wkt_geom = 'POLYGON ((5.123456 51.6, 6.8 52.6, 7.8 51.6, 5.123456 51.6))' #DOES NOT WORK

layer = mapnik.Layer('wkt_layer_bla','+proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs')
layer.datasource =  mapnik.Geos(wkt=wkt_geom)
layer.styles.append('My Style')
mmap.layers.append(layer)
mmap.zoom_all()
mapnik.render_to_file(mmap,'/foo/bar/bla.png', 'png')


texts: Anyone have any ideas why this tiny code snippet throws a COM exception?

texts:
System.Object obj = Activator.CreateInstance(Type.GetTypeFromProgID("esriLocation.LocatorManager"));
ILocatorManager2 locatorManager2 = obj as ILocatorManager2;
ILocatorWorkspace locatorWorkspace = locatorManager2.GetLocatorWorkspaceFromPath(@"Z:\ESRI\DataMaps10\streetmap_na");
locator = locatorWorkspace.GetLocator("Street_Addresses_US");
IAddressGeocoding addressGeocoding = (IAddressGeocoding)locator;  //this  fails

Inspecting the addressCoding object gives

'addressGeoCoding.matchFields' threw an exception of type Interop.COMException


texts: Interrupting ModelBuilder between two processes?

texts: I am using a model builder to create a tool that uses a DEM to create contour lines. Then I want to select, on screen, one of the contour lines to use it as an input parameter for the clip tool to clip the DEM. Is it possible to interrupt the model builder tool in run time so the user can select a specific contour lines and then the tool continues to the next step (clipping).


texts:
In ArcMap, it is very easy to reload table cache where it will call the CSV that was previously used to create the feature class to update the feature with any new data from the same CSV.
In ArcGIS Pro, it seems to me that once we create a feature layer from a CSV file, the attribute table of the feature layer is independent to the CSV file used. There doesn’t seem to be a way to reload table cache in ArcGIS Pro.
I’d like to be able to have the feature class updated using the same CSV (with new data added) instead of me recreating the feature every time we collect new data. This is so that I can preserve the symbology as well.
I’ve considered appending, but as some entries from previous import can change, I prefer to overwrite the attribute table in the feature class instead.
If this is possible using ArcGIS Pro or ArcPy, then how?
Edit: I tested with ArcMap and ArcGIS Pro, once the feature class is created, the attribute table is independent of the csv and reload table cache in ArcMap/ refresh in ArcGIS Pro will not update with the latest csv. So I guess there is not really away where I can update the attribute tables of a feature class without recreating the features using a new version of the csv?


texts: Update a postgis table with the resuly of SQL Executor

texts:
I have built a workflow in FME2015 in which I use SQLExecutor to run a query. I want to use the result of the query to UPDATE a column of table in my database.
I am having difficulties in passing the results into the table using a postgis writer with a writer mode in UPDATE.
I attach a screenshot of my workflow:

The column I want to update is called: nettype. And I need to update it based on the key: fr_gid. If I understand correctly I need to use the fme_where attribute but I am not sure in what way.


texts: How to get 3D GPS positions from geographic latitude, longitude and altitude?

texts: I have a list of GPS measurements for a point for a period of about 2 years. They are latitude, longitude and altitude. I want to generate the time series for the east, north and up components. Do I just need to convert the geographic coordinates to UTM coordinates or should I do something else?


texts: Adding custom projection to GDAL library?

texts:
Can i add my custom projection definition to GDAL library in CentOS?

For example, if i want to add custom projection def. in PostGIS i use simple SQL command for inserting projection. But i do not find how i add projection to GDAL library


texts: Consuming JSON directly to ArcGIS Online

texts:
I have a JSON API with a token information that I want to integrate with ArcGIS Online to create a heatmap. It's a news feed API that gives latitude, longitude, country, news in brief and news in detail.

I am unable to find a way to directly use the API to consume the JSON. Is there a way I can use the JSON directly or consume it in ArcGIS Server without using the Geoevent Server(used for real-time GIS)?


texts: SELECT geometry FROM spaces WHERE geometry ARE WITHIN OR OVERLAP polygons FROM OTHER table

texts:
I've got two tables, one containing LINESTRINGs and one containing POLYGONEs.

I would like to select all the LINESTRINGs that are containing within the outer boundary of all the POLYGONs combined.

I'm running PostgreSQL 9.3 and PostGIS 2.1.

Is that possible from within the WHERE statement?


texts: Making ArcGIS tool to generate ObjectIDs that start from Upper left rather than bottom left in Fishnet?

texts:
When I spatial join a fishnet with another fishnet, I see ObjectIDs start from bottom left of the feature rather than Upper left of features. How can I make ArcGIS to create feature class having records starting from Upper left of features?

Images:
What is happening:

What I need:

I know a tool "Spatial Sort" which does this. But I don't want to use that as I already have to apply many geo-processing tools on large data, and make it slower.

Any Idea?


texts: Fail to get heatmap from MULTIPOINT column of Spatialite layer in QGIS?

texts: I got a heatmap from a POINT column of spatialite layer. But, when I made another heatmap from MULTIPOINT column, it was blank. What is the possible reason of this problem?


texts: Rename feature classes in file geodatabase in bulk

texts:
I see may versions of this post but none quite get me to where I need to be.

I have 50,000 features classes in a file geodb.

They are a1_xxxxx_xx

I need to rename them all to a6_xxxxx_xxx

This is my code. Does not seem to quite work. Does create many a6 but a1 still remain when it completes.

import arcpy
from arcpy import *

env.workspace = r"C:\IND_Rivers\New\Outputs\A3\a3_03.gdb"

fcList = arcpy.ListFeatureClasses()

#attempts to replace feature classes with  "a1" with "a6"
for fc in fcList:
    Rename_management(fc, fc.replace('a1','a6'))


texts:
I have photo created by UAV, and I want to calculate bNDVI:

  (NIR-blue)/(NIR+blue)

I have these bands (filter is ok -100%). So using QGIS I have values -0,1-> 0,3 - but bNDVI must be from 0,65->0,80 (look at picture):

I have problems with the calculation; on QGIS, ArcGIS, and pix4D I get same result. I thought that I had to multiply val/256*2-1 but didn't work again. Color depth is 24. Can anyone help me?


texts: GeoServer Web Feature Service (WFS) in OpenLayers problem

texts:
I am running GeoServer 2.1.1, and have a datastore from ArcSDE populated with features projected in EPSG:26918. I can successfully add a WMS to an OpenLayers map, but am having issues with WFS. I have the map set to EPSG:900913 in preparation for a mashup with Google base layers.

I am under the impression that vector layers are reprojected on the fly, and that absent of any user specification, features will be rendered using the default style.

Here is the main chunk of my script:

    map = new OpenLayers.Map('map', {
        maxExtent: new OpenLayers.Bounds(
            -128 * 156543.0339,
            -128 * 156543.0339,
            128 * 156543.0339,
            128 * 156543.0339),
        maxResolution: 156543.0039,
        projection: "EPSG:900913",
        units: "m"
    });

    var bridgeLayer = new OpenLayers.Layer.WMS(
        "Bridge",
        "http://localhost/geoserver/geodatabase/wms", {
        layers: "geodatabase:Bridge",
        format: "image/png",
        transparent: true
    },{
        isBaseLayer: true
    });

    var wfsLayer = new OpenLayers.Layer.Vector(
        "Photo Points", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.WFS({
                version: "1.0.0",
                url: "http://localhost/geoserver/wfs",
                featureType: "Photo Point",
                srsName:"EPSG:26918",
                featureNS: "http://dev.geodatabase.org",
                geometryName: "SHAPE"
            }),
            visibility: true
    });

    map.addLayers([bridgeLayer, wfsLayer]);

And so, the WMS displays correctly (and when I add a Google base layer, it lines up correctly). However, the WFS layer fails to appear.

According to Firebug, however, my POST request is successful, and I receive an response. Furthermore, the layer object I defined looks to be a fully formed layer, although it does not appear to have any features in its "features" array. Not sure what that means. Screenshots appear below.

Something also that I don't understand, when I switch versions to "1.1.0", my POST request fails with a "NoApplicableCode" exception.

Additionally, I have tried reprojecting the data in GeoServer to spherical mercator EPSG:900918, and also redefining the srsName parameter, but both experiments have failed.


texts: QGIS expression for setting label visibility by scale

texts:
I am using QGIS 3.x and struggling with the labels and I would like to have some expression for setting the visibility of some layers based by scale.

How can I make some of them disappear?
I tried:
   if(@map_scale<2000) then ("fid" = 18) = NULL

where "fid" is my data attribute table column
NULL - label invisible when scaling smaller than 2000.
It doesn't work

According to the answer below, I used the following expression:
  @map_scale<2000 AND FID = 12

and I see, that the label I wanted to disappear is displayed, whereas all others I want to keep are gone.

How can I turn it around?


texts:
How can i convert a field 1,2,3,4... to 001,002,003,004...010...?

For example

1 change to 001
2 change to 002
3 change to 003

...

Is there any formula in field calculator?


texts: What is the information behind this geometry?

texts:
I am new to PostGIS. I have a geometry value. And I have used the following query to that geometry.

SELECT (ST_AsText(the_geom)) FROM (SELECT (ST_Dump(ST_Boundary(geom_original))).geom As the_geom FROM sample_table) As geom_text LIMIT 2

I have got the following result:

POINT Z (623566.485 6961791.436 96.825)
POINT Z (623599.006 6961805.281 94.874)

I want the explanation of this result. I have assumed something like that this geometry contains z(elevation) component also. 96.825 and 94.874 is showing the elevation. The first two values are boundary x and y values. But there should be two set of terminal points i.e. start and end and it is giving only one set of points. How to get start and point of this boundary? I have tried st_startpoint and st_endpoint but these are giving null value.

I have also executed st_astext query:

SELECT ST_AsText(geom_original) FROM public.sample_table As geom_text LIMIT 2

The result of the st_astext query is:

MULTILINESTRING ZM ((635205.589 7068587.674 226.674 0,635211.789 7068585.293 226.361 6.641,635222.143 7068578.926 226.15 18.796,635231.575 7068572.259 225.988 30.347,635233.497 7068566.092 226.011 36.806,635228.315 7068558.871 225.969 45.694,635220.565 7068558.979 226.061 53.445,635213.056 7068563.573 226.317 62.248,635208.348 7068574.825 226.304 74.445,635205.589 7068587.674 226.674 87.587))
MULTILINESTRING ZM ((623566.485 6961791.436 96.825 0,623576.041 6961796.754 95.624 10.936,623584.711 6961800.301 94.832 20.304,623599.006 6961805.281 94.874 35.441))

But the problem is that I am unable to correlate the st_astext and the st_astext of st_boundary result. And how to get exact coordinate value from these values?


texts: Remove/Fill holes in Polygons from a MultiPolygon using JTS Topology Suite programmatically

texts:
There are several questions/answers using ERIS,etc.
like:
Remove holes from polygon

But I am having a challenge applying the solution programmatically using the JTS Topology suite only.

Below is an example of MultiPolygon.. I want to remove/fill the holes.


texts:
I've written a python script to get tiles from, for example, starting here:
http://otile1.mqcdn.com/tiles/1.0.0/map/1/1/1.jpg

I've read that map servers will, understandably, limit the amount of downloading one can perform, though I haven't run into problems yet. I am interested here in OSM MapQuest, but will also be interested in USGS Topo and Satellite data.
How can discover what the limits are?


texts: Vertical datum query for coastal modeling

texts:
I have a USGS derived DEM with a vertical datum of EGM96 from NASA satellite data. My Bathymetry uses a lowest astronomical tide datum.

Has anybody come across this or similar issue before?

I need to correct both to a constant but the transformation I can access doesn't seem to recognize EGM96.

I'm using ArcMap 10.2


texts: I've looked through the PostGIS raster documentation, but can't find an indication of how to apply a convolution filter to raster data.  Does this capability even exist (yet?) There only appears to be preset operators such as min, max and std deviation.


texts: Synchronizing coordinates from GPS to central server

texts:
I have developed a tracking system using a GPS device connected to a Linux box. It collects the coordinates as it progress. It's done with a python script, which is listening to the GPS device and write the coordinates to a CSV file when it sees new reading by the GPS device.

My scenario is to use multiple of these Linux boxes each connected to a GPS device. I need all collected coordinates to be synchronized to a central server in realtime. The purpose of this is for, mapping all the devices in in realtime in the central server.

My question is about, how and what are the technologies can I use to synchronize the coordinates from multiple devices to the central server in realtime.

As of my research, I have came across of two primary way.

Using RSYNC - synchronize the collected data periodically to upstream with Linux's rsync
Push the data through HTTP API

There are 2 problem with above mentioned methods, one is, one is this types of synchronization can not be set to work on realtime. Always there will be a delay of sync interval. Second is, I am afraid, if the sync interval is very low, the central server will become resource hungry.

Am I missing a 3rd idea or one of these both will fit to my work?


texts: How to add habitat information from shapefile to an existing DF dependant upon grid reference?

texts:
I'm currently building a dataframe (DF) using R and I'm wanting to add habitat type as a variable. Within DF I have grid reference points as variables.

I've imported a shapefile containing habitat data using rgdal::readOGR.

I'm able to plot the shapefile using the base plot function but I've been unable to add the points from DF and unable to find information anywhere on how to then extract which habitat those points fall into.


texts: Incremental network loader Error - Geometric netwotk QA

texts:
Trying to merge 2 geometric networks with Incremental network loader I got this error -

How can I find this Closed polylines at my geometric network ?


texts: Where is list of Geographic Transformations in ArcGIS, and their areas of applicability?

texts:
I have some data in Provisional South American Datum 1956, and I need to project it to Web Mercator with WGS84 as the datum.

When I use the project tool in ArcGIS desktop, it requires a Geographic Transformation, and offers me 13 different ones: from PSAD_1956_to_WGS_1984_1 to PSAD_1956_to_WGS_1984_13

How do I find out which is the transformation meant for my area?

I remember finding a pdf on the Esri Forums  several years ago, which showed the various transformations between NAD 83/27 & WGS 84.

Is there something similar for all transformations in ArcGIS?


texts: Disconnect between raster symbology in Layers Panel and Layer Properties (QGIS 3)

texts:
I'm working on a Python plugin for QGIS that applies a customised colorRampItemList to a single band raster. The result displays correctly on the map and in the Layers Panel but incorrectly in the Symbology tab of Layer Properties. I can load the desired color ramp into that tab using setSourceColorRamp but the values and labels are different, and reflect what you get when clicking on Classify in Continuous mode. Is there any way of preventing the customised values/labels from being overridden? The fact that they can be manually edited and retained in Layer Properties implies that there should be.

The disconnect can be replicated in the Python Console by adapting lines 195-212 of test_qgsrasterlayer.py from the QGIS code on GitHub, e.g.:

myRasterLayer = iface.activeLayer()
myRasterShader = QgsRasterShader()
myColorRampShader = QgsColorRampShader()
myColorRampShader.setColorRampType(QgsColorRampShader.Interpolated)
myItems = []
myItem = QgsColorRampShader.ColorRampItem(10, QColor('#ffff00'), 'foo')
myItems.append(myItem)
myItem = QgsColorRampShader.ColorRampItem(100, QColor('#ff00ff'), 'bar')
myItems.append(myItem)
myItem = QgsColorRampShader.ColorRampItem(1000, QColor('#00ff00'), 'kazam')
myItems.append(myItem)
myColorRampShader.setColorRampItemList(myItems)
myRasterShader.setRasterShaderFunction(myColorRampShader)
myPseudoRenderer = QgsSingleBandPseudoColorRenderer(myRasterLayer.dataProvider(), 1, myRasterShader)
myRasterLayer.setRenderer(myPseudoRenderer)
myRasterLayer.triggerRepaint()

Currently working in QGIS 3.6.0 but have had this problem since migrating the plugin from 2.18 to 3.


texts: Ice type differentiation using satellite data

texts: I would like to classify glacial area by categories of snow (if there is) and ice, but what's most important: between old ice and fresh ice. They have different properties that are possible to recognize in field, but can you do this with satellite data? (preferably Landsat becasue of 30/15m spatial resolution)


texts: Labeling each side of polygon with its length in ArcGIS Pro?

texts: I am trying to figure out a way in ArcGIS Pro to label each side of a polygon with the length of the side without breaking it into individual polylines.


texts: Dividing map into features (whole squares)

texts: I'm using QGIS 1.8.0. I divided a map of a state with its municipalities into many features using a GRID (small squares), but instead of creating one ID per feature (one ID per square), it creates one ID per feature considering also the features divided by the municipalities lines, which I don´t want. I mean if a square is divided by the borderline between two municipalities, it considers two features.  I would like to have the map of the state with its municipalities divided into many "whole" squares, and each and every one of those squares should have one unique ID. Could anyone help me?


texts: Sentinel-1 missing data gap in adjacent acquisitions

texts:
Some months ago I posted Export each image from a collection in Google Earth Engine about exporting Sentinel 1 scenes from a collection and was successful with help from Rodrigo Principe.

I now have a new region of interest, which sits on the overlap between 2 scene footprints (sequential images in descending orbit). When I make a mosaic of the 2 scenes I get a small no data gap and lack of alignment between the 2 images:

// Create a geometry representing an export region.
var roi = ee.Geometry.Rectangle([-62.67, -8.31, -63.25, -8.85]);
var zoomPoint = ee.Geometry.Point([-62.35, -8.06])

// Load the Sentinel-1 ImageCollection.
var sentinel1 = ee.ImageCollection('COPERNICUS/S1_GRD')
.filterBounds(roi);

// Filter by metadata properties.
var IW_H = sentinel1

// Filter to get images with VV and VH dual polarization.
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
.filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))

// Filter to get images collected in interferometric wide swath mode.
.filter(ee.Filter.eq('instrumentMode', 'IW'))

// Filter IWs to get High res.
.filter(ee.Filter.eq('resolution', 'H'))

// Filter IW-Highs to get 10m res
.filter(ee.Filter.eq('resolution_meters', 10));

// Filter to get images from different look angles
var col_desc = IW_H.filter(ee.Filter.eq('orbitProperties_pass', 'DESCENDING'));

// Spatially mosaic the images in the collection and display
var col = col_desc.mosaic();
var colVV = col.select(0)
Map.addLayer(colVV, {min:-17, max:-3}, "VV");
Map.centerObject(zoomPoint, 20);

What am I missing?  Are these scenes not acquired sequentially?
https://code.earthengine.google.com/ca7a4b585f1600b771a2de07c8e75792


texts: Understanding map projection coordinate systems

texts:
I have been playing around with a variety of map projections for the purposes of transforming between various representations.

The thing that I still don't quite understand is how to interpret the (x,y) values of the projection. Take the Mercator projection, for example. I get that (x,y) represent locations on the tangent plane, except that they are still in spherical coordinates: e.g. x = lambda - lambda_0.

How do these (x,y) values translate to actual locations on the tangent plane (i.e. the map) if they're still in spherical coordinates? In my case, I ultimately need to identify actual pixel locations when I render the new projections, but generally speaking it is not clear to me what these spherical (x,y) values mean with regard to the mapping surface.

-----Update-----

Perhaps this additional info may help to explain my confusion. I found the image below on slide 60 of this powerpoint. I think it illustrates the step I am confused about. Here it's called the reference transformation. I've been reading through Snyder's "Map Projections: A Working Manual" and I think I've concluded that the rectangular coordinates output from the forward projections are arbitrary in scale. Where I am struggling is how to go from that arbitrary scale to my pixel grid.

-----Update 2-----

I hope it may be more clear if I explain how I have arrived at my confusion.

When resampling an image according to some function (e.g. "undistorting" an image, which is kind of like an inverse projection of some image distortion), the standard algorithm is:

(1) Apply the inverse projection to the normalized pixel coordinates of the image

(2) Find the min/max bounds of that projection and create a grid in that space

(3) Apply the forward projection to the normalized coordinates of that grid to determine where in the original image to sample

(4) Sample according to the output of step (3) in order to render the new image.

I have been trying to follow these steps for converting between an equirectangular map and various projections. It seems my mistake is that, in the algorithm I laid out above, the projected space is always a pixel coordinate system.

Am I correct to understand that each map projection has its own rectangular coordinate system? That is, Mercator rectangular coordinates mean something different than Gnomonic rectangular coordinates, etc. If so, is it based on the class of projection (planar, cylindrical, conic), or specific to the individual projection? Snyder lays out rectangular coordinates in tables in his manual, but they don't seem to correlate with each other.


texts: How to select points as geometry in QGIS when both multipolygons and points are stated within the wfs-server?

texts:
I'm trying to show data from the following server in QGIS: http://geodata.nationaalgeoregister.nl/bag/wfs?service=WFS&request=GetCapabilities&

I made a filter on the data I would like to appear in my map. This is how I stated the wfs-url in QGIS:

http://geodata.nationaalgeoregister.nl/bag/wfs?VERSION=1.0.0&cql_filter=postcode LIKE '6828%'

The following link is an example on how the data appears in xml-format (different filter):

http://geodata.nationaalgeoregister.nl/bag/wfs?SERVICE=WFS&version=2.0.0&request=GetFeature&typename=bagviewer:verblijfsobject&count=10&srsname=EPSG:28992&bbox=189869.57826037705,444557.7565065846,190175.32637351775,444863.5046197253

In the last link it show two kinds of geometry, one is called geometrie (points) and the other pandgeomotrie (multipolygon). The geometry appears in polygons on the map, though I would like the points to appear on the map.

Is there a way of doing this?


texts: Sextante in QGIS - downloaded, but does not show up

texts:
I downloaded Sextante, and I put it into the directory

C:\users[my_user_directory].qgis\python\plugins\sextante-1.0.9\sextant

I re-started the computer, but there is still no sign of Sextante in QGIS, neither in plug-ins manager, nor menu items.

Do I have one directory layer too many in the file path above? Is that why it doesn't work?

Thanks.


texts: Running classification on image collection using Google Earth Engine?

texts:
I have ImageCollection. composed of satellite images which shows Zebras on grasslands. The location of each zebra is stored within a point shapefile, which I uploaded as a "table" into Earth Engine.

I wish to extract the pixel values under each point within the shapefile (i.e. the pixel values for each zebra) and then run a supervised classification to detect the zebras.

When I try to run the code on my ImageCollection, this is the error that is returned:

   table2 = ee.FeatureCollection("users/example/abc..."),
    imageCollection = ee.ImageCollection("users/example/wsx...")

// Load training polygons from a Fusion Table.
        // The 'class' property stores known class labels.
        var polygons = ee.FeatureCollection(table2);

// Get the values for all pixels in each polygon in the training.
var training = imageCollection.sampleRegions({
  // Get the sample from the polygons FeatureCollection.
  collection: polygons,
  // Keep this list of properties from the polygons.
  properties: ['class'],
  // Set the scale to get Landsat pixels in the polygons.
  scale: 30
});

  imageCollection.sampleRegions is not a function


texts: Using field parameters in calculate field

texts:
I am attempting to make a generic tool to move existing points to an already defined XY.  What I want to do is have the user define the feature class, then the X field and Y field from that feature class, then run the Calculate Field with those X and Y fields.  Below is what I have, but I know I am going about it the wrong way (fairly new to ModelBuilder).  I can't find a way to associate those first Field Parameters with the feature class.

I know I can just use the same X and Y field names every time and it will work, but I want the tool to be as versatile as possible.


texts: Compute NDVI over each image in an image collection and export NDVI thresholded images

texts:
I want to compute NDVI for all images in Sentinel2 (less than 20% cloud cover) in an image collection comprising 5 years. Thereafter, I want to batch export the output as a binary TIF file where NDVI <0.5 is ignored. Once I can do this, I will do the rest of the analysis in QGIS through a model that is ready to go.
Here's what I've done so far in bits and pieces to the best of my abilities.
//add study area
var studyarea = ee.FeatureCollection('users/XXX/XXXXX/studyarea');

Map.addLayer(studyarea, {color: 'red'}, 'study area');
Map.centerObject(studyarea);

//
/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

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

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

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

// Load the Sentinel-2 ImageCollection and bound to study area
var collectionSentinel2 = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2015-07-01', '2020-10-15')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .filterBounds(studyarea)
                  .map(maskS2clouds)
                  .map(function(image){return image.clip(studyarea)});

//Mapping a Function over a Collection
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

//Applying NDVI to all images in the collection
var withNDVI = collectionSentinel2.map(addNDVI);

print (withNDVI);

So far so good apparently, but this is adding NDVI as a separate band, not sure if that's a good approach. Thereafter, here's the next bit of the code where I am thresholding NDVI images. I'm unable to use this as a function as with computing NDVI earlier.
//thresholding to extract healthy vegetation
//Identify all pixels below threshold and set them equal to 1. All other pixels set to 0
var ndvi_thrshld = ndvi.gt(0.50);

//Remove all pixels equal to 0
ndvi_thrshld = ndvi_thrshld.updateMask(ndvi_thrshld);

//Clip to study area
var ndvi_thrshld_clipped = ndvi_thrshld.clip(studyarea);

To export to Drive, I am using the following code snippet which works great for a few images but not sure how to iterate it over the image collection.
// Export the image, specifying scale and region.
Export.image.toDrive({
  image: ndvi_thrshld_clipped ,
  description: 'NDVI',
  scale: 10,
  region: studyarea,
  maxPixels: 18095956539
});


texts:
When re-saving a GeoJSON file to a shapefile, the .dbf significantly increases in size to the point where the file size becomes unmanageable.

For example, I re-saved a 250 MB GeoJSON file to shapefile and the resulting .dbf file in the shapefile grew to 2 GB.

Why is this?

I'm using QGIS 3.6.3.


texts: Routine for Transforming Raster Layer Projection

texts:
I am working with georeferenced raster images in pseudo mercator (3857). They are in TIFF with TFW- and JPG with JGW world files. The TIFFs are also available tiled.
I'd like to be able to transform either of these file types to OS/GB 27700 crs (preferably the JPG/JGW). While retaining the file format as is.

Is it possible to do this in JPG without changing to TIFF as an intermediate step?

If not, is there a way to transform the TFW file along with the new TIFF?

My experiments so far have been of mixed success. Black images, or images that are only a tiny part of the original. Through the "save as"- and "raster translate" routes.

I am now working with the JPGs. I can get the raster translate to work - creates images that have been warped. However they are produced with a new black outline, which isn't ideal (though can be cropped externally I would suppose).

Is there a way to carry-out the geoTIFF warp without producing the black cutline?

That said, main issue is I need a JGW file to accompany these renders.

Is there a routine for creating JGW file for an existing JPG?

QGIS example is here:
https://www.dropbox.com/s/zzqfw7hoi2v3tjv/Basic%20JPGs%200.1.qgs?dl=0
And JPG file folder (including failed renders), here:
https://www.dropbox.com/sh/bxxmcvb2wrb9gl4/AACi6S6GKGPpcLUWT0o_-nKAa?dl=0


texts: PyQgis how could I get the id from a selected feature

texts:
using self.iface.mapCanvas().selectionChanged.connect(self.getselectionid) how could I get the id of a manually selected feature on the Qgis map specifically on a vectorlayer? so it can be used to program a function.

I tried a couple of things to see if I could figure what selectionChanged is returning, these are my attempts and its results:

print("I selected {}".format(str(selection.feature.id())))
>>I selected <qgis._core.QgsVectorLayer object at 0x7f9d9f189c18>

print("I selected {}".format(str(selection[0].id())))
>>I selected <qgis._core.QgsVectorLayer object at 0x7f9d9f189c18>
>>I selected project_s_shapefile_f403eb35_14b9_4e61_b327_263dce57a8cb

print("I selected {}".format(str(selection)))
>>I selected <qgis._core.QgsVectorLayer object at 0x7f9d9f189c18>


texts: I am trying to get an idea of the average most recent imagery dates of Madagascar in Google Earth. Rather than zooming in and looking at the lower corner for the imagery date, what is the best way to find out how much of the country is covered by imagery from a given year?


texts: OSGEO Advanced Install Problems with Print Composer

texts:
I would like to have the ability to read/write file gdb data and export my data as a GEOpdf. I have always done the advanced install option so that I can read/write file gdb data.  However, now I need to get version 2.18, so that I can export geopdf from print composer.

My problem is that when installing the advanced option, I am forced to install the LTR version 2.14.16 due to some dependencies for different plug ins such as grass, etc.  This creates a problem in that when I go to create a map (in version 2.18), I do not have the ability to reference it in the Print Composer options. Or when I create in a map in 2.14.16, I do not have the ability to georeference.

Any Ideas out there


texts: Consuming gml generated by sql server As.GML()

texts:
Sql server has a method As.GML to produce GML file from geography/geometric type. Is there a good example to utilize this GML on other mapping platform such as QGIS? I have ask a related question here, but QGIS reports its format is not valid.

Thanks.


texts: Multiply raster with a variable ArcMap 10.1

texts: I want to multiply my raster with a value that would be a parameter in ModelBuilder.


texts: Adding vertex set distance from first or last vertex in polyline?

texts:
For ArcGIS 10.1, is there a python script or tool that will add only one vertex (node) a set distance from the first vertex, or alternatively, a set distance from the last vertex for multiple polylines?  The original vertices need to be retained for each polyline.

I have over 500 polylines and creating vertices individually in an edit session is taking forever!


texts: Shapefile of United Kingdom Constituencies of 1970 Elections

texts:
Does anyone know where to find shapefiles of older UK constituencies?

I am specifically looking for the shapefile of the constituencies used for the 1970 elections (of which the results can be found here: http://www.politicsresources.net/area/uk/ge70/ge70index.htm )

I found the shapefile of the current constituencies here: https://www.ordnancesurvey.co.uk/business-and-government/products/boundary-line.html


texts: QGIS WMS Layer returned legend text is distorted

texts:
QGIS version
3.16.6-Hannover
Mac OS installer.
On loading this (https://environment.data.gov.uk/spatialdata/survey-index-files/wms) WMS layer with the default settings I get a distorted legend for the possible layers. I'd like help debugging this. Is it a system setting in my OS, a setting within QGIS or something else?

EDIT: The space available to the legend seems to be squashed by this bit of window but I can't see anything on the GUI to resize it


texts: Identify parcels that have two types of land use

texts: I have two polygon layers in a geodatabase. One layer is land parcels. The other layer is land use. I want to identify the six largest parcels that contain both "upland forests" and "wetlands" (i.e. WHERE landuse = 4 AND 6 in a single parcels). How do I proceed???


texts: Is there any open source library for reading OCAD files?

texts:
I would like to read map files created by OCAD (.ocd). The file format specification is open and available on the web page, but it's binary and complex.

Is there an open source library that can be used for reading these .ocd files?


texts:
I have installed the OTB 5.10.1 version and when i want to use train images classifier for the random forest it doesn't create any output model file. and the library give me this error message:

ERROR: Waiting for at least one parameter...
This is the TrainImagesClassifier application, version 5.10.1
Train a classifier from multiple pairs of images and training vector data.

Complete documentation: http://www.orfeo-toolbox.org/Applications/TrainImagesClassifier.html

Parameters:
-progress Report progress
MISSING -io.il Input Image List (mandatory)
MISSING -io.vd Input Vector Data List (mandatory)
-io.valid Validation Vector Data List (optional, off by default)
-io.imstat Input XML image statistics file (optional, off by default)
-io.confmatout Output confusion matrix (optional, off by default)
MISSING -io.out Output model (mandatory)
-elev.dem DEM directory (optional, off by default)
-elev.geoid Geoid File (optional, off by default)
-elev.default Default elevation (mandatory, default value is 0)
-sample.mt Maximum training sample size per class (mandatory, default value is 1000)
-sample.mv Maximum validation sample size per class (mandatory, default value is 1000)
-sample.bm Bound sample number by minimum (mandatory, default value is 1)
-sample.vtr Training and validation sample ratio (mandatory, default value is 0.5)
-sample.vfn Field Name (mandatory, default value is )
-classifier Classifier to use for the training [libsvm] (mandatory, default value is libsvm)
-classifier.libsvm.k SVM Kernel Type [linear/rbf/poly/sigmoid] (mandatory, default value is linear)
-classifier.libsvm.m SVM Model Type [csvc/nusvc/oneclass] (mandatory, default value is csvc)
-classifier.libsvm.c Cost parameter C (mandatory, default value is 1)
-classifier.libsvm.opt Parameters optimization (optional, off by default)
-classifier.libsvm.prob Probability estimation (optional, off by default)
-rand set user defined seed (optional, off by default)
-ram Available RAM (Mb) (optional, off by default, default value is 128)
-cleanup Temporary files cleaning (optional, on by default)
-inxml Load otb application from xml file (optional, off by default)

It give me only the libsvm classification.


texts: How to display symbols on a Raster in QGIS?

texts:
I have a raster file in which each cell has a numeric code which indicates that the cell is directed towards one of its eight neighbours.

I am using QGIS to display each of the eight possible codes using a different colour right now, but wonder if there is a way to display symbols, such as little arrows, as well, especially at high zoom levels.


texts:
I'm wondering if it is possible to do pattern recognition with QGIS. I have 5 NDVI Landsat images and I try to find areas automatically (constant pattern). Is there a tool in QGIS?

The pic below shows three different raster layers from different years. You can recognize typical areas (e.g. bright or dark blue pattern) which I want to identify automatically.


texts: WFS works only in IE browser. Any work around is available for this issue?

texts:
WFS works only in IE browser for me.

It's not working in any other browser. I don't have much expertise in  checking browser compatibility, could someone suggest me how to resolve this issue or if any workaround is available for this issue to make it works in Mozilla and Chrome also.


texts: Data (type double) mismatch between ArcMap and ArcGIS Server 10.6.1

texts:
When I identify a feature which is having a double value as 0.00001 in ArcMap, it is showing what is expected. I have the layer published as a map service using ArcGIS Server. If I search for the same feature through Query Task of Map Service, I am getting the value as 9.999999999999999E-6. Below are some values where I identified such behavior.
Value in ArcMap : -123.0102,
Value in ArcGIS Server Rest Endpoint : -123.01020000000001
Value in ArcMap : 0.00001,
Value in ArcGIS Server Rest Endpoint : 9.999999999999999E-6
Value in ArcMap : -1.23456789,
Value in ArcGIS Server Rest Endpoint : -1.2345678900000001
Field Details from ArcMap: Number Format:  Numeric, Data Type: Double, Precision: 38, Scale:                    8
Products: ArcGIS Desktop 10.6.1, ArcGIS Server 10.6.1, Oracle DB


texts: QGIS Profile from Line tool requiring shapely.wkd?

texts: I am trying to install the Profile from Line tool for QGIS, and the system informed me I need shapely.wkd. A search on the QGIS site for shapely takes me to the Contour Plugin...can anyone help me out with this?


texts: Exporting BlenderGIS model maintaining georeferenced coordinates not possible

texts: Using Blender 2.83. Importing georeferenced data from other GIS System works well. Generated model with BlenderGIS. Problem: Exporting planes and bodies as xx.obj files does not maintain georeferenced coordinates. Importing the .obj files back to Blender 2.83 works but importing into other GIS does only preserve local x,y,z coordinates and not the georeferenced coordinates. Appears as shortcoming to use BlenderGIS. Only exporting georeferenced shape files possible.