texts: Greater control of period on ui.DateSlider in Earth Engine

texts:
I am using ui.DateSlider() on Google Earth Engine.
The period argument for that object seems to be the number of days (ref: https://developers.google.com/earth-engine/api_docs).

I would like to control the period to be monthly, or semi-monthly (say 1st-14th, and 15th-end-of-month). Are these two configurations possible?

I am unable to achieve it with the period being restricted to number of days as the months have different days.

Sample code:

// Use a DateSlider to create annual composites of this collection.
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1');
// Use the start of the collection and now to bound the slider.
var start = ee.Image(collection.first()).date().get('year').format();
var now = Date.now();
var end = ee.Date(now).format();

// Run this function on a change of the dateSlider.
var showMosaic = function(range) {
  var mosaic = ee.Algorithms.Landsat.simpleComposite({
    collection: collection.filterDate(range.start(), range.end())
  });
  // Asynchronously compute the name of the composite.  Display it.
  range.start().get('year').evaluate(function(name) {
    var visParams = {bands: ['B4', 'B3', 'B2'], max: 100};
    var layer = ui.Map.Layer(mosaic, visParams, name + ' composite');
    Map.layers().set(0, layer);
  });
};

// Asynchronously compute the date range and show the slider.
var dateRange = ee.DateRange(start, end).evaluate(function(range) {
  var dateSlider = ui.DateSlider({
    start: range['dates'][0],
    end: range['dates'][1],
    value: null,
    period: 365, //<= I would like to change this to monthly or semi-monthly
    onChange: showMosaic
  });
  Map.add(dateSlider.setValue(now));
});

How do I change the period to be monthly or semi-monthly?


texts: OpenLayers SVG sprite externalGraphic only works in Firefox

texts:
I am trying to set the externalGraphic property of a vector layer's style to an SVG sprite that has various viewboxes to position the icons dynamically by creating different views within my SVG file:

<view id="bikeview" viewBox="0 0 52 52"></view>
<view id="boatview" viewBox="52 0 52 52"></view>

...and then accessing them like so:

var defPoiStyle = {
    externalGraphic: 'images/activities/activitySprite.svg#boatview',
    pointRadius: 30,
    cursor: "pointer"
}

The problem is that this only works in Firefox, although with Chrome and Safari I can access each view by entering the URL directly into the browser like this, just not in the map.

I have also noticed that in Chrome's network log (and likely Safari's) that the #boatview part of my image path gets stripped out, but not so in Firefox. I assume it is some kind of XML request issue or something, but I feel like I got my declarations right:

<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
  "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
<svg
   xmlns:svg="http://www.w3.org/2000/svg"
   xmlns="http://www.w3.org/2000/svg"
   xmlns:xlink="http://www.w3.org/1999/xlink"
   x="0px"
   y="0px"
   width="52px"
   height="52px"
   id="svg3753"
   version="1.1">

This is too complicated for jsFiddle, but here is my map. Give it a go in Chrome and then in Firefox, and you should see what I mean (the boat icon is what is supposed to display).


texts: I'm using OpenScales in a web mapping project, and I would like to add a Google Maps layer. I know that there's some licence problems, but I don't care for these licence issues, I'm looking to add Google Maps layer to OpenScales like it has been made on OpenLayers.


texts: How to create Relationship Classes successfully

texts:
I'm new to ArcGIS and I'm making the first steps at this time. At the moment I've two different shapefiles. The first shapefile represents streets on a map with some detailed information in the attribute table. The 2nd shapefile represents crossroads of the same map including some detailed information in the attribute table too. Both attribute tables has IDs and logically it’s a tree where to node is a crossroad and the edge is a street.

Now the next step is to create a Relationship Class to "link" both shape files for future algorithms.  But when I want to create a Relationship Class I receive different errors. I've read this document before and tried to create such a class using an Add-in but I received a NotImplementedException.

To be sure that creating such a class is possible I tried to create a relationship class which ArcMap using ArcToolbox -> Data Management Tools -> Relationship Classes -> Create Relationship Class manually. I can choose both shapefiles but when I have to choose the Output Relationship Class I receive the error: “000177 : Relationship classes can only be created in the same database that contains the feature classes or tables.“

Can someone give me some notes what I've done wrong and tips to solve this problem?

Edit:

I've created a Relationship Class with ArcCatalog, a new personal database and the imported shapefiles. But I can't open this table and control if the table was successfully created and has content. Does anyone have had a similar problem?

Edit2:

What I've done:

1) Create a new Personal Geodatabase using ArcCatalog10

2) Import 2 Feature Classes (existing shapefiles filled with data)

3) Start creating a Relationship Class using ArcToolBox -> Create Relationship Class

4) For testing I created a one-to-one relation with this two imported shapefiles

5) The table was created successfully but:

5.1) If I want to import this in ArcMap I receive the message "Warning: A selected Item could not added to the map" and nothing happened

5.2) The preview function in ArcCatalog doesn't work for this Relationship Class

In summary I can't see if the Relationship Class is logically right and I can't see which data are selected and can be used for future algorithms.


texts: Handling ArcMap closing event in AddIn to save edits automatically

texts:
I'm trying to handle the closing event of ArcMap (10.0 SP5) within an ESRI Addin Extension using C# and Arcobjects. Furthermore I would like to handle 'a closing event', that is fired BEFORE ArcMap tries to stop a started edit session. As background: I want to create an extension that does the database handling automatically without asking the user if he wants to save the changes in the database when closing ArcMap (this should happen programmatically).

At first I tried to handle the ESRI.ArcGIS.ArcMapUI.IDocumentEvents.BeforeCloseDocument, but this event will be fired AFTER ArcMap tries to stop an active edit session. The situation is not different with the ESRI.ArcGIS.Editor.IEditEvents2_Event.BeforeStopEditingevent. In both cases the user is asked if he wants to save the changes before one can handle the event. And at last I tried the ESRI.ArcGIS.Desktop.AddIns.Extension.OnShutdown event, but this is fired much more later after the user is asked for saving the changes in the Geodatabase and saving the MXD. Now I am at my wit’s end.

Does anyone else has an idea or a workaround? Have I missed another event? Is it possible to ‘cast’ the ArcMap application window to System.Windows.Forms.Form to wire the FormClosing event?


texts: Multiple AND statements in QGIS Query Builder

texts:
In the QGIS Query Builder I am trying to filter by multiple fields. My filter expression is as follows:

"VENTILTYPE" IN ('Skydeventil', 'Trykreduktionsventil') AND  "STATUS" = 'I brug'  AND "AABNINGSGRAD"='100'

If I leave out one of the statements, everything works fine.
The following two expressions works as they should:

"VENTILTYPE" IN ('Skydeventil', 'Trykreduktionsventil') AND  "STATUS" = 'I brug'
"VENTILTYPE" IN ('Skydeventil', 'Trykreduktionsventil') AND "AABNINGSGRAD"='100'

However as soon as I add an extra QGIS gives me the following error:

  2019-07-31T16:13:28     CRITICAL    Layer Åben : Download of features for layer ugis:TL298 failed or partially failed: Download of features failed: Error transferring [Link for WFS] - server replied: Not Found. You may attempt reloading the layer with F5

What is causing this error?
Can the Query builder not handle 3 statements?


texts: DBF file type changed to .txt - how do i return to original source

texts: I have opened a dbf file as txt to edit some of the data, now I cannot find the original default programme to open the dbf file with. I tried changing the default programme for the dbf but the QGIS folder has hundreds of options. Can somebody point me to the right programme? Currently, when I open the file in QGIS it appears to have the attribute table beside it as opposed to the icon as usually.


texts:
I have created a web map using CARTODB API. I'm wanting to make an animation map with my points. The map points displays properly, but the time slider will not display as well as no animation of points.
Here's my HTML:

<!DOCTYPE html>

    <title>Easy example | CartoDB.js</title>
    <meta name="viewport" content="initial-scale=1.0, user-scalable=no" />
    <meta http-equiv="content-type" content="text/html; charset=UTF-8"/>
    <link rel="shortcut icon" href="http://cartodb.com/assets/favicon.ico" />

    <link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.15/themes/css/cartodb.css" />

I have a feeling that I am having the same issue as they did here:
Torque map slider formatting problem

Does anyone have a clue as to why my time slider won't display in my map?


texts: Export a raster from a WMS layer in QGIS Graphical Modeller

texts:
I would like to be able to export a section of a WMS to a georeferenced image using the graphical modeller.
I already have the extent setup for other things but I can't find an algorithm that will work.
I have tried "r.out.png" and "clip raster by extent" with no joy.


texts:
I'having trouble using writeOGR to write GeoJSON files, when there are MultiPolygon objects. Here is an example: I'm creating a SpatialPolygonsDataFrame object from GeoJSON, with a single feature made of two polygons. When I export to GeoJSON, the file that's written has a single feature with one polygon and a hole.

Why is this, and how can I persuade writeOGR to write two polygons?

# Create SpatialPolygonsDataFrame with one feature, two polygons
spdf <- readOGR(layer='OGRGeoJSON', verbose=FALSE, '{
"type": "FeatureCollection",
"features": [{
 "type": "Feature",
 "id": 0,
 "properties": {"label":"CENTRAL"},
 "geometry": {
     "type": "MultiPolygon",
     "coordinates": [[[[102.0, 2.0], [103.0, 2.0], [103.0, 3.0], [102.0, 3.0], [102.0, 2.0]]],
                     [[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 1.0], [100.0, 0.0]]]]
     } }] }')

# When I save and reload this object, coordinates has length 1,
# i.e. the feature is a polygon and a hole.
fn <- tempfile()
writeOGR(spdf, fn, layer='whatever', driver='GeoJSON')
length(fromJSON(file=fn)$features[[1]]$geometry$coordinates)

I'm using R 3.1.2 and rgdal 0.9.1 on OSX Yosemite.

UPDATE: Here is the GeoJSON tempfile that is created. As you can see, the coordinates array has been modified from the original.

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
"features": [
    { "type": "Feature",
      "id": 0,
      "properties": { "label": "CENTRAL" },
      "geometry": {
          "type": "MultiPolygon",
          "coordinates": [ [ [ [ 102.0, 2.0 ], [ 102.0, 3.0 ], [ 103.0, 3.0 ], [ 103.0, 2.0 ], [ 102.0, 2.0 ] ],
                             [ [ 100.0, 0.0 ], [ 100.0, 1.0 ], [ 101.0, 1.0 ], [ 101.0, 0.0 ], [ 100.0, 0.0 ] ] ] ] } }
    ]
}


texts: What does 'to the right' mean, regarding the >> operator in PostGIS?

texts:
I am reading the documentation on the >> operator in PostGIS, and I cannot understand the meaning of strictly to the right because on the plane, right or left depends on the direction at which one is looking.

Reference: https://postgis.net/docs/ST_Geometry_Right.html

What would be that direction?


texts: How can I see the queries (SQL statements) that are generated when I modify or edit the records of a PostGIS table, inside QGIS?

texts:
This is the DB manager, and it shows only the history of the statements that were used, but it does not show the records of the modified data in a table.
How do I see the history of the statements made in the database when modifying records?


texts: Setting pixels at edge of raster to NoData in FME

texts:
I am trying to generate GeoTiff rasters using FME Desktop 2015.1.  My input is an 8 bit raster that has areas with pixel values of 0 around its edges.  I want to set these to NoData.

I tried using the RasterBandNodataSetter to set the values to NoData. This was technically successful but pixels with a value of 0 in the middle of my raster were set to NoData as well which is not what I want.

I have also experimented with using the clipper to extract the areas I want hoping that the areas I do not extract would be set to NoData.  However the clipper transformer uses the rectangular extent of my clipper feature to clip while I need it to be tapered to the extent of my study area.

I have also tried to write an alpha8 band to the GeoTiff in addition to my main Gray8 band. I was hoping to set the alpha8 to 0 where there is no data so nothing would be displayed. However the GeoTiff writer seems to insist on turning by alpha8 band into a UINT8.

What is the best approach to achieve this?


texts: What is a good learning source for programmer who wants to understand the infrastrucre of map based web map?

texts:
I want to to write a simple web app with a map (tile cache because of high resolution raster) that once a user clicks the map he can do some queries on it (as in search in polygon, radius, nearest...) on the TAB (Not Google Maps or Bing, but my own tab I've got from a GIS friend).

From the little I've read I've encountered many new Terms as: Web Map Server / Web Map Service / Tile cache / WFS etc...

I'm looking for free (open source probably) solution that I can learn and use in .Net. I know client side will be OpenLayers but still need some reading material that will make some sense (and explanation) on all these terms - architecture. I know there is GeoServer, SharpMap, MapServer, TileCache... just not sure what each does / used for...


texts:
I have a pre-tiled image for some area, that

Tiles are in JPG format.
Tiles do not overlap.
Each tile has its own world file.
All tiles have the same projection.

Now, I want to consume those tiles as a tile service,

Is there any way to make such service without making a mosaic, i.e. tile service from the tile images immediately.
Should my tiles follow some scheme standard for naming and storage directory structure.
I prefer a solution based on ESRI/ArcGIS technologies, but any other solution will be useful.


texts: Accessing NOAA CO-OPS water levels using SOS with bounding box and time extent

texts:
I'm tring to access all water level data within a bounding box and temporal extent using the OGC SOS from the NOAA CO-OPS.  On the page here: http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/ there is a very useful SOS query builder tool, and I can formulate a query to get all the station data in a specific region:

http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?service=SOS&request=GetObservation&version=1.0.0&observedProperty=water_surface_height_above_reference_datum&offering=urn:ioos:network:NOAA.NOS.CO-OPS:WaterLevelActive&featureOfInterest=BBOX:-74.4751,40.389,-73.7432,40.9397&responseFormat=text%2Fcsv&eventTime=2012-10-26T00:00:00Z/2012-11-01T00:00:00Z&result=VerticalDatum%3D%3Durn:ogc:def:datum:epsg::5103&dataType=VerifiedSixMinute

But when I execute this query it tells me max of 1 day can be requested in a collection.

So I guess there are two approaches.

Approach 1: Do a query on collection with bounding box, but only the start time requested:

http://opendap.co-ops.nos.noaa.gov/ioos-dif-sos/SOS?service=SOS&request=GetObservation&version=1.0.0&observedProperty=water_surface_height_above_reference_datum&offering=urn:ioos:network:NOAA.NOS.CO-OPS:WaterLevelActive&featureOfInterest=BBOX:-74.4751,40.389,-73.7432,40.9397&responseFormat=text%2Fcsv&eventTime=2012-10-26T00:00:00Z&result=VerticalDatum%3D%3Durn:ogc:def:datum:epsg::5103&dataType=VerifiedSixMinute

which returns:

station_id,sensor_id,"latitude (degree)","longitude (degree)",date_time,"water_surface_height_above_reference_datum (m)",datum_id,"vertical_position (m)"
urn:ioos:station:NOAA.NOS.CO-OPS:8518750,urn:ioos:sensor:NOAA.NOS.CO-OPS:8518750:W1,40.7006,-74.0142,2012-10-26T00:00:00Z,0.446,urn:ogc:def:datum:epsg::5103,1.848
urn:ioos:station:NOAA.NOS.CO-OPS:8531680,urn:ioos:sensor:NOAA.NOS.CO-OPS:8531680:W1,40.4669,-74.0094,2012-10-26T00:00:00Z,0.319,urn:ogc:def:datum:epsg::5103,1.624

and then I could strip out the station_id and use them in single station queries.

Approach 2: Issue multiple collection request for periods of one day and then concatenate the responses together.

I guess I would use approach 1 if I had more days than stations, and approach 2 if I had more stations than days.

Is there a better way?


texts: Cannot see created SQL Server tables in ArcGIS

texts:
I have a MSSQL 2012 database that I can connect to from ArcGIS using 'Add Database Connection'. In this database I have a spatial table that I can see and open in ArcGIS, however when I create a table using 'SELECT * INTO newtable from..' I cannot see the new table in ArcGIS. I can see the new table in Management studio and it has all the columns and data. The only difference I can see from the old table is that there is no Primary key. If I set this it makes no difference and I still cannot see the table in ArcGIS. Here is the SQL script I'm using:

USE TrafficData;

IF object_id('TrafficData.dbo.a') IS NOT NULL
DROP TABLE TrafficData.dbo.a;

SELECT *

INTO TrafficData.dbo.a

FROM AMBASE_V15LINKS_POLYLINE;


texts: Area calculation in "Thomas Fire analysis"

texts:
I am studying this notebook: California wildfires 2017 - Thomas Fire analysis
I can follow the first part of the notebook. But I have got confused in this cell to calculate the burnt area:
ext = {"xmax": -13246079.10806628, "ymin": 4035733.9433013694, "xmin": -13438700.419344831, "ymax": 4158033.188557592,
       "spatialReference": {"wkid": 102100, "latestWkid": 3857}, "type": "extent"}
pixx = (ext['xmax'] - ext['xmin']) / 1200.0
pixy = (ext['ymax'] - ext['ymin']) / 450.0

res = burnt_areas.compute_histograms(ext, pixel_size={'x':pixx, 'y':pixy})

numpix = 0
histogram = res['histograms'][0]['counts'][1:]
for i in histogram:
    numpix += i

So where do the numbers 1200 and 450 in the pixx and pixy calculation come from? I have tried to investigate the metadata of the object and cannot find any hints. Are they some constant numbers?
Also, why the compute_histograms returns five classes while "burnt_areas" has only four "output_values" classes?


texts: AttributeError: DescribeData: Method FIDSet does not exist error

texts:
I am writing an add-in that in one step requires to find all layers that have selections. To achieve that I am using FIDSet property of Describe object. Unfortunately, I get an error: "AttributeError: DescribeData: Method FIDSet does not exist" (why "method" anyway?). While I perform test like this, in Python window in ArcMap:

desc = arcpy.Describe("O_PODST")
desc.FIDSet()

Everything works fine, but while using my add-in I get this error.
How to make it work?
Here is my code:

class ShowParcLU(object):
    """Implementation for Okienko_addin.button (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False

    def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")

        if arcpy.mapping.ListTableViews(mxd, "PARC_LU"):
            layers_with_selections = []

            for i in arcpy.mapping.ListLayers(mxd):
                if len([int(fid) for fid in arcpy.Describe(i).FIDSet.split(";") if fid != '']) > 0:
                    layers_with_selections.append(i)


texts: Showing labels just for a specific point type

texts: If I have a bunch of points all on the same layer and they are all categorised by different types in the Attribute Table, is there an easy way to just show labels for one type and not the others? Maybe through an expression?


texts: Calculating maximum distance between roads (Spatial Line Shapefile) in QGIS?

texts:
I’m writing my masters thesis and came across a problem I haven’t been able to find an answer to.

I’m working with QGIS and a Spatial Line Shapefile that represents all streets in Tibet, taken from Open Street Map. What I want to do is create a (raster?) file of Tibet from it with a resolution of, for example 1x1 m, that represents the distance of that 1x1 m area to the closest street. My endgoal is to calculate a correlation of the distance to the closest road with NDVI data that indicates the health of vegetation.

I was thinking about applying buffers to the lines but then I would have to have the buffers clip so that the border is at the exactle middle distance between the lines and it wouldn't show the distance along the buffer as far as I know.

Unfortunately I am at a loss of how to do this. I am also working with R. I don't know if it might be more useful for my problem.


texts:
I have code that creates either a UniqueValueRenderer or a ClassBreaksRenderer and applies it. I also have code that uses AnnotateLayerProperties to set labels.

When I create a renderer and then an annotation or an annotation and then a renderer the map goes haywire - sometimes one displays and the other one is missing and rarely (without pattern) both display.

Must I use the OnDraw event (and save the annotation instructions) to ensure that annotations are reproduced?  What can I do about the renderer?  Can I save that and restore it in the OnDraw? In what order should these be drawn?


texts:
What I have:
A large TIFF covering parts of several contiguous sheets of the British OS 6 inch first edition.  It has no georeferencing.  It's just a series of scans which someone has stitched together (actually gdal_transform extracted it from a larger .ecw)

A number of scans of original sheets.  These are also plain scans with no georeferencing but they do have latitudes & longtitudes along the borders so I can read off the corner coordinates in DMS.

I assume that the projection of all this data will be close to if not identical to the current OS projection.

What I want to achieve:
The large TIFF properly georeferenced so it can be used as a basis for georeferencing various tithe maps.

How I expected to go about it:
Set the existing CRS of one of the scanned original sheets to OSGB70 which, as far as I can see is the OS projection using lat/long, fire up the georeferencer, set GCPs at each corner, set the output to OSGB70 & run the transform.  Repeat for other sheets.  Using these georeferenced sheets pick off GCPs for features on the big map to use the georeferencer to georeference this.

What happened:
A grossly distorted map was produced.  I think I can see what the problem is; the original is in the OS projection but its CRS is purely pixels & the georeferencer treated it as degrees.  But what do I set as the initial CRS to set up the correct transform?

TIA


texts: Projecting data from 'EPSG:4326' to Lambert Azimuthal Equal Area with PyProj

texts:
I have a list of data points given in lon/lat coordinates, which I want to convert to cartesian ones using pyproj. My dataset refers to Italy, so I was considering using either the Lambert Azimuthal Equal Area projection or the Albers Conic Equal Area one. For now I am trying out the former, but since I am not very familiar with coordinate systems, I have run into a basic question with respect to the definition of the datum, ellipsoid and earth radius.

According to what I understand so far, each ellipsoid has a standard radius value. For instance, for the WGS84 ellipsoid the radius is R=6378137. However, when I specify both the ellipsoid and R value in pyproj.proj, I get different coordinates compared to when I only define the ellipsoid.

from pyproj import Proj, transform, Geod
import numpy as np

lambert_aea1 = {'proj': 'laea',
          'lat_0':43.2947,
          'lon_0':12.1295,
          'x_0':0.,
          'y_0':0.,
          'ellps': 'WGS84',
          'datum': 'WGS84',
          'R':6378137.0}

lambert_aea2 = {'proj': 'laea',
          'lat_0':43.2947,
          'lon_0':12.1295,
          'x_0':0.,
          'y_0':0.,
          'ellps': 'WGS84',
          'datum': 'WGS84'}

xi = [12.1295, 8.4555, 11.1193, 15.8035, 13.1496]
yi = [43.2947, 44.8834, 47.1653, 41.7059, 39.4241]
inProj = Proj(init = 'epsg:4326')
outProj1 = Proj(lambert_aea1)
outProj2 = Proj(lambert_aea2)
x1,y1 = np.array(transform(inProj,outProj1,xi,yi))
x2,y2 = np.array(transform(inProj,outProj2,xi,yi))

>> x1 = ([0., -290653.25071928, -76769.35932204, 306156.50224268, 88008.51279864])
>> x2 = ([0., -290185.4370178, -76641.501629, 305687.406043, 87878.46410642])
>> y1 = ([-21378.23435213, 161863.38128257, 409910.18924861, -191449.19444042, -451253.91439645])
>> y2 = ([0., 182924.98967511, 430536.02515148, -169835.95115785, -429281.93580092])

Now if I calculate the distances between the first point and all other points, and compare them to the distances I get by using the function pyproj.Geod.inv, I get the following differences:

geod = Geod(ellps='WGS84')
dist1 = []
for jj in range(1,len(xi)):
    dist1.append(geod.inv(xi[0],yi[0],xi[jj],yi[jj])[2])
dist2 = np.sqrt((x1[0]-x1[1:])**2 +(y1[0]-y1[1:])**2)
dist3 = np.sqrt((x2[0]-x2[1:])**2 +(y2[0]-y2[1:])**2)
print (dist2-dist1)
print (dist3-dist1)

>> [ 517.81025651  660.65516468  485.62015613  537.69336351]
>> [ -46.78192359 -102.50618725  -38.57430363  -70.06941143]

What am I missing here?

Also are these distance "errors" reasonable? This is the first time I deal with coordinate systems and projections, so even though this precision is probably fine for my purposes, I was surprised to find such errors for my relatively small region.

To rephrase my question, I am wondering why does my 'laea' projection using pyproj change, if in addition to the ellipsoid (WGS84), I specify the radius of the earth R = 6378137. From what I understand, this is the radius dictated by the WGS84 ellipsoid anyways. Is pyproj using some default radius despite me defining the ellipsoid, or am I doing something wrong?

I also cannot find any good documentation for pyproj and proj4 with information about the arguments needed for each specific projection type.


texts: Calculate area and perimeter of coordinate points in Python

texts:
I have a set of points that make up a polygon. I need to find the area and perimeter for the shape in Python. I have used the shapely library to calculate the area and perimeter but it is nonsensical because the points are latitude and longitude coordinates (in US States).
I need to convert the points to an equal-area projection system first. I was thinking of using Sinusoidal projection. I was able to do this and calculate area (code below), but I'm not sure how to interpret the result (also is it right?). What unit is it in? (m²?)
Next I need to calculate perimeter, but would the Sinusoidal projection distort it and throw it off? How would I deal with that?
I am a computer scientist and my work on a project has led me down this rabbit hole. How can I calculate the area/perimeter of this shape in the correct way?
def sinusoidal(lat, long):
    earth_radius = 6371009 # in meters
    lat_dist = pi * earth_radius / 180.0

    y = lat * lat_dist
    x = long * lat_dist * cos(radians(lat))

    return x, y


texts: linux gdal 3.1.0 ogr2ogr ERROR 1: Unable to find driver PostgreSQL

texts:
I have seen a number of posts similar to this but I have another issue. I try to complie gdal-3.1.0 to import osm.pbf to PostgreSQL. My os is Red Hat Enterprise Linux Server release 7.6. I need Oracle and PostgreSQL driver. Oracle is working but PostgreSQL throw error Unable to find driver PostgreSQL.
I install PostgreSQL client libpqxx-devel-5.0.1-2.rhel7.1.x86_64 and set enviromen variables:
LD_LIBRARY_PATH=/opt/sas/shared/devdi94/thirdparty/oracle/12.2.0.1/client/lib:/usr/pgsql-11/lib:$LD_LIBRARY_PATH
ORACLE_HOME=/opt/sas/shared/devdi94/thirdparty/oracle/12.2.0.1/client
PKG_CONFIG_PATH=/usr/pgsql-11/lib/pkgconfig

export LD_LIBRARY_PATH ORACLE_HOME PKG_CONFIG_PATH

Then, ./configure --with-threads --with-pg=yes --disable-static
results looks good
OCI support:               yes
PostgreSQL support:        yes

make & make install Log is clear without errors.
Ogr2ogr is working for oracle but not for PostgreSQL
ogr2ogr -f PostgreSQL PG:"dbname='xxx' host='xxx' port='xxx' user='xxx' ACTIVE_SCHEMA='xxx' password='xxx'" "dolnoslaskie-latest.osm.pbf"

>>> ERROR 1: Unable to find driver `PostgreSQL'.

I don't know where the problem is. have I some missing library or enviroment variable?

I try install version gdal-3.0.1 and both drivers working. Maybe gdal-3.1.0 have bug.


texts:
Everytime I get a new scene, I have to build a specific number of band combinations. It's almost 50 band combinations. I have to do it for almost every scene, over and over.

Is there any automatic process that allows me to build those band combinations automatically?

I am using ENVI+IDL 5.1


texts: Getting R to recognize buffer around las file

texts:
I'm trying to clip from a catalog using a shapefile. The file is huge, so I've divided it into two (ctgw1 and ctgw2). The final radius should be 47 (the ctg is in feet) and I'm trying to set a buffer of 328ft (100m). This all runs fine and creates valid las files, but when I run it, I receive the error code:

ERROR buffer greater than 0 is needed

ctg <- catalog(folder = "D:/All LAS")
shp.w <-shapefile("C:/Users/elliot59/Downloads/plot4w/plot4w.shp")

bufferw <- lasclip(ctg, shp.w, radius = 328.084)

plots <- seq(1:244)
for(n in 1:length(plots)){
  plot.n <- plots[n]
  if(bufferw[[plot.n]]@header@PHB$`Number of point records` == 0) {
    next
  }
  writeLAS(bufferw[[plot.n]], file = paste('plot', plot.n, ".las", sep = ""))
}
ctgw1 <- catalog(folder = "C:/Users/elliot59/Desktop/buffer.w1")
ctgw2 <- catalog(folder = "C:/Users/elliot59/Desktop/buffer.w2")

opt_output_files(ctgw1) <- ("C:/Users/elliot59/Desktop/Out")
opt_laz_compression(ctgw1) <- TRUE

thinnedctgw <- lasfilterdecimate(ctgw1, random(2))

dtm <- grid_terrain(thinnedctgw, 1, kriging(k = 10L))


texts: Difference between Geoserver stores: ImageMosaic, WorldImage, ImagePyramid

texts: What are the differences between the Geoserver stores ImageMosaic, WorldImage, and ImagePyramid? I noticed that ImagePyramid is no longer included in the latest Geoserver. Is one of these more efficient than the others?


texts: Cannot federate the server at this time. 10.3.1

texts:
When I go to federate the ArcGIS Server I get the following error message:

  "Unable to federate at this time"

When I look into both the Portal and AGS logs I get the same error message:

  "Failed to update the security configuration. Cannot update the
  security configuration as server is not registered with Portal".

When I set the AGS logs to DEBUG I get this message:

  " Testing to see if this server has been registered with Portal as a
  federated server. Failed to test if the server is registered with
  Portal"

Then it states: "File or directory not found" in relation to http://www.w3.org/1999/xhtml and this of course cannot be accessed as it is a completely disconnected environment.

I have also checked the following at 10.3.1:

I was able to add store credentials of a secured ArcGIS Server service from the AGS I am trying to federate against so that proves that the certificates work as expected.
ROOTs and INTERMEDIATEs have been imported into the Portal Java keystore.

The main confusion stems from the fact it was federated earlier this morning, the URLs I used were similar to these (with the actual names removed for security reasons)

SERVICES URL: https://machinename.domain/arcgis

ADMIN URL:    https://machinename.domain:6443/arcgis

The correct PSA details were used in the username and password.


texts: QGIS label calculation always renders zero

texts:
I'm trying to use the area of a paddock as its label. The layer is a shapefile containing a single polygon. In Layer Properties/Labels/Label this layer with I define an expression:

concat(tostring(round($area/10000,2)),'Ha')

The preview correctly yields 4.67Ha but canvas renders 0Ha.

The shapefile .prj file:

PROJCS["NAD83_NSRS2007_UTM_zone_59N",GEOGCS["GCS_NAD83(NSRS2007)",DATUM["D_NAD83_National_Spatial_Reference_System_2007",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",171],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]

The measure tool ellipsoid is Clarke 1866

QGIS 2.8.3/Windows.


texts: How to change a marker color based on a rule in Qgis

texts: I am styling point markers in a layer that I imported using "Add Delimited Text Layer". I want the style markers based on the color marker. So each color relates to a different species (I have two different species).


texts: Field Calculator text field to new shapefile

texts:
I have point features being selected by a polygon feature and copied to their own shapefile.  I also have certain values of the polygon feature being applied to the newly created point file.  The point file already has the required fields that need to be populated.  The integer/floating fields are copying correctly, but the lone string field (SfcType)/row[4] value isn't being applied.  I realize this has to do with the fact that it is a string, but I can't seem to find the right syntax.

Can anyone point out my (obvious) mistake?

from arcpy import * import os
arcpy.env.overwriteOutput = True

#Make feature layer of feature classes to allow selection
MakeFeatureLayer_management (inPolyfc, "Polylyr")
MakeFeatureLayer_management (inPointfc, "Pointlyr")

#Get name of OID field for sql expression
polyOIDfld = Describe ("Polylyr").OIDFieldName

#iterate through polygons with cursor
cursor = da.SearchCursor ("Polylyr", ["OID@", PolyIDfld, "Slope", "endelev", "SfcType", "NearID"])
for row in cursor:
#sql to select individual polygon
sql = '"{0}" = {1}'.format (polyOIDfld, str(row[0]))
    #select individual polygon
SelectLayerByAttribute_management ("Polylyr", "", sql)

#select intersecting points
SelectLayerByLocation_management ("Pointlyr", "", "Polylyr")

#Check to make sure points are selected, ie at least one point actually intersects polygon
if Describe ("pointlyr").FIDSet:
    FeatureClassToFeatureClass_conversion ("Pointlyr", outPath, str(row[1]))
    print "Created Shapefile Named: " + str(row[1])
    Fc = outPath + "/"+ str(row[1]) + ".shp"
    arcpy.CalculateField_management (Fc, "Slope", str(row[2]))
    print "Added SLope to " + Fc
    Fc = outPath+ "/" + str(row[1]) + ".shp"
    arcpy.CalculateField_management (Fc, "endelev", str(row[3]))
    print "Added Rwy End Elev to " + Fc
    Fc = outPath+ "/" + str(row[1]) + ".shp"
    arcpy.CalculateField_management (Fc, "SfcType", str(row[4]))
    print "Added Sfc Type to " + Fc
    Fc = outPath+ "/" + str(row[1]) + ".shp"
    arcpy.CalculateField_management (Fc, "NearID", str(row[5]))
    print "Added NearID to " + Fc
    #Create empty feature class if no points intersect
else:
    print "No points were under the " + str(row[1]) + " surface"


texts: Creating Strip Map Index in Portrait view for Data Driven Pages?

texts: I'm trying to generate a strip map with a north to south orientation. Every time I try to create a strip map for this polyline I get pages in landscape orientation. I would like to have a portrait orientation as my layout for the final product. I've tried changing the page orientation to vertical with no luck as well.Any ideas on how I can fix this issue.


texts:
I am trying to find an easy way to have private tables, yet exposing them to certain clients if they are authenticated.

This authentication can be achieved by a token system that allows the client to perform queries only if the token or the credentials are valid.

The approach here was to have a function with the authentication logic and whatever queries or operations you need to expose, and call that function from the client along with his credentials.

Problem with this is that you end up with a pg function for every operation, need to precisely define the input and output of the function, no polimorphy, so every time anything changes (like a data type) all functions need to be recreated. Not very flexible.

So another approach would be to have just one function that acts as a wrapper for whatever SQL you send it, checks credentials, sanitizes the sql and sends back the response.

CREATE OR REPLACE FUNCTION runSQL(thesql text,client text, hash text) RETURNS SETOF RECORD as $$
DECLARE
sql text;
expiry time;
client_info RECORD;
m   varchar;
--here go the blacklisted strings
blacklist varchar[] := array['UPDATE ','INSERT ', 'ALTER ', 'DROP ', 'TRUNCATE', 'DELETE '];
BEGIN
    expiry:=INTERVAL '1 hour';
    sql:= 'SELECT * FROM auth WHERE client_id = $1 AND client_hash = $2 AND created_at > (NOW() - $3)';
        EXECUTE sql using client, hash, expiry INTO client_info ;
            -- check if it's a valid user and has an up to date token
        IF client_info IS NULL THEN
            RAISE EXCEPTION 'Authorization failed for client %', client;
        END IF;

    sql:= thesql;
            -- sanitize the sql
    FOREACH m  IN ARRAY blacklist
            LOOP
            IF position(m in upper(thesql)) = 1 OR position(concat(' ', m) in upper(thesql))>0 THEN
                    RAISE EXCEPTION 'Not enough privilleges';
            END IF;
            END LOOP;
            --send back results
        RETURN QUERY EXECUTE sql ;
END;
$$ LANGUAGE plpgsql
SECURITY DEFINER;

And usage would be:

SELECT * FROM runSQL('SELECT cartodb_id, the_geom FROM sometable', 'client_id', 'token') as f(cartodb_id int, geom geometry)

My concern here is not opening a huge security hole in my database. Do you think it would be enough just to check for potentially dangerous strings in the query (like INSERT, DELETE, DROP, ALTER, GRANT, REVOKE...) or are there other security issues with this approach?

Any other ideas on how to have a SECURITY DEFINER function that only allows certain operations like read, or access to just some tables?

Thanks!

David


texts: I am trying to open the scripter or scripting tool in Whitebox GAT so that I can edit Export GeoTIFF tool to make Whitebox export a TIFF file from dep/tas format, which was suggested by Whitebox developer in another post. But nothing pops up when click the tool. I am in Windows 10 and I have Java the newest version.


texts: Splitting street in Address table which participates in Relationship Class (1-M)?

texts:
I have a table "Addresses" and the Feature "streets".
The feature "Streets" are lines.

I need to do a relation with this table with this feature. What I do?

I tried built a "relationship class". The type is "1-M". It is the line of a street with the table of a address. But happens an error. I can´t split the line of a street. Anyone knows how to fix it?

That is the error:
"

  The feature could not be split. The field is not nullable.

"

I found the error "FDO_E_FIELD_NOT_NULLABLE" description on a ESRI site. But they don't told how to fix it. This is the link:
http://edndoc.esri.com/arcobjects/8.3/ComponentHelp/esriCore/fdoError.htm

I need the relationship and to split the streets!
How to solve/fix this error?

The real problem is: I have a Feature of streets, called "Eixos". And a Table of Addresses, called "Logradouros". I tried with Domains. But when I run a script python that make "Table to Domain" give me the error of "

  ERROR 000464: Cannot get exclusive schema lock

". But I cannot kill the other users. Do you know a good solution to do a relationship Feature Street vs Table Address?


texts: Trying to use bringToFront to order combined tile layers layergroup for leaflet basemap

texts:
I've combined a from a geoserver WMS tilelayer (trails) seperately with a mapBox streets-baseic and mapbox sateline using LayerGroup to form two basemaps, "Streets" and "Imagery" available for a leaflet map. I'm wanting the WMS trails to become part of these base maps and not be an overlay users can control. I'm new to leaflet and am having trouble with ordering of the TileLayers within the streets layer group.

The map starts with the trails on top of the streets basemap fine. The trails will still display over the sateline imagery when switching to the Imagery base map.  However, when the base map is switched back over to the Streets base map, the trails disappear.

I've tried to implement trails.addTo(map).bringToFront(); after the layer control, thinking it is an ordering issue and even played around with zIndex, without success. I've looked at the follwoing two posts. I could try to remove and re-add the layers each time the base map is switched, but it seems the answer should be simpler.

How to change layer order in leaflet.js

https://stackoverflow.com/questions/12848812/layer-ordering-in-leaflet-js

Can anyone tell me how I can order the combined tile layers within the base map?

Here is my code

"use strict";

//map
var map = L.map('map', {
        center: [IVY.default_latitude, IVY.default_longitude],
        zoom: 14
    });

//default basemap layer

//street layer using mapbox basic
var mapboxbasic = L.tileLayer('https://api.mapbox.com/v4/mapbox.streets-basic/{z}/{x}/{y}.png?access_token=' + IVY.mapbox_api_key, {
        maxZoom: 18
    });

//City Trails
var trails = L.tileLayer.wms('https://bloomington.in.gov/geoserver/publicgis/wms', {
        layers: 'publicgis:BikePedTrails',
        format: 'image/png',
        transparent: true,
        version: '1.1.1'
    });

//aerial layer using mapbox satellite
var mapboxsat = L.tileLayer('https://api.mapbox.com/v4/mapbox.streets-satellite/{z}/{x}/{y}.png?access_token=' + IVY.mapbox_api_key, {
        maxZoom: 18
    })
;

//Base Map Layer Groups

var streets = new L.layerGroup([mapboxbasic, trails]).addTo(map);
var imagery = new L.layerGroup([mapboxsat, trails]);

/*
//experiment use bringToFront inside layergroup def - no difference
var streets = new L.layerGroup();
    mapboxbasic.addTo(streets).addTo(map);
    trails.addTo(streets).addTo(map).bringToFront();

var imagery = new L.layerGroup();
    mapboxsat.addTo(imagery);
    trails.addTo(imagery).addTo(map).bringToFront();
*/

//layer controls

var basemap = { 'Streets': streets, 'Aerials': imagery}
;

L.control.layers(basemap).addTo(map);

// bring trails to front

//base layer change method

//only works when switching to aerials - not when switching to streets
map.on('baselayerchange', function(layer) {
trails.bringToFront();
});

//these don't make a difference
//trails.addTo(map).bringToFront();
//trails.addTo(basemap).bringToFront();


texts: Snap lines (map matched trajectories) to road network in PostGIS

texts:
Using OSRM I have map matched a series of GPS points to a custom road network. My goal is to provide insight in where but more importantly how many times the vehicle has visited each road segment in the road network by counting the number of overlapping map matched trajectories per road segment. However, since OSRM only stores 6 decimal places of fixed precision for coordinates (about 10cm precision at the equator) and my roadnetwork uses a 7-decimal places fixed precision format (about 1cm at the equator), there's up to 9cm of error introduced by that truncation. Therefore, I first need to snap the map matched trajectories (red) to the nearest road segment (green).

My desired output is a table snapped_lines in which the map matched trajectories are snapped to the road network (i.e. the red lines lay exactly on top of the green lines). Since the map matched trajectories are not too far away from the road network, I was hoping ST_SNAP should work. Executing this query takes ages however, so I was wondering whether something is wrong.
DROP TABLE IF EXISTS snapped_lines;
CREATE TABLE IF NOT EXISTS snapped_lines AS (
  SELECT
    ST_Snap(
      a.line_geom,
      b.geom,
      1.1) as geom
  FROM
    mapmachted_trajectory AS a
    roadnetwork AS b);


texts: Using python string wild cards within the field calculator ArcGIS 10.3.1

texts:
I have two fields. One field has an ID. The other field has a list of ID's. Using the Field Calculator I with to calculate 1 or 0 if either field contain the same ID number. I was successful at writing an arcpy row by row query, but there's millions or records and a linear test of rows is not efficient. I need to calculate it in one hit.

I need to take a field and make it's value a wild card and test the other field for that value. I have tried the "*" and "in" methods. But it's not returning the correct value.

Can someone assist me in why this is wrong?

Code One:

def myCal(reportnumber,spatialjoin):
 reportnumbercheck = ('*'+(str(reportnumber))+'*')
 if (reportnumbercheck in spatialjoin):
  return 1
 else:
  return 0


texts: Looking for a places API that can return residential and business addresses for a city

texts:
I'm working on a web application for a UK-based client that has two location functions -

1) report a location using keyword or postcode

2) search for a specific place using keyword or postcode

The results need to be limited to the client's home city, but should include residential addresses, so users can report their own homes if needed. I have not been able to find an API that does exactly what I want.

I have tried Google Places Search and Autocomplete, but it only returns businesses and POIs - residential addresses are excluded.

I have also tried OS Places, but Find returns results for the whole country, and Nearest is limited to 1000m radius, so cannot cover a whole city.

So, the question is, am I using the above APIs correctly, or have I missed something that means I can actually search residential addresses within a city (~50km radius).

Or, is there another API I could be using that does what I want?


texts: Which D3 microlibraries are required to load TopoJSON/GeoJSON files and translate into path definitions?

texts:
D3 v4 is modular, based on around 30 "microlibraries", intended to allow discrete functions to be broken out and used on their own.

I'm trying to assemble a minimal set of D3 microlibraries to do the following:

Asynchronously load a TopoJSON file from the server using .json function
Process this data
Using a projection, turn it into a SVG polygon or path definition strings (not SVG elements)

...and that's it. I'm using something else (Raphael.js) to process the path strings, so I don't need the rest of D3.

I thought I'd only need topojson and d3.geo but I'm either not understanding something or I'm lacking the capability to load JSON files and the capability to convert data to path definitions, and I'm struggling to find those features in the D3 micro-libraries.


texts:
I would like to move the cursor to a given x,y coordinates of the MapCanvas.

Is there an easy way to do this using PyQGIS?


texts: Online source of polygon nodes for individual newly created zipcodes

texts: When a new zipcode is created it would be very convenient to obtain its polygon nodes from an online source, without having to download a national dataset. Does an online web-page or web-API exist, where you can provide a zipcode as a search-term and get back its shapefile nodes expressed as a series of latitude/longitude pairs?


texts:
I just installed QGIS 3.4 on my computer (macOS Sierra 10.12.6) and I'm trying to install QuickMapServices, but I'm unable to. My plugins window looks like this:

and the settings tab says "The settings on this tab are only applicable for Python Plugins. No Python support detected, thus no settings available":

I downloaded QuickMapServices from the QGIS Python Plugin Repository, but I can't access it. I have Python 3.6.0 on my computer.


texts: No FID in attribute table when creating a fishnet

texts:
I have a vector polygon and I want to create a 3 x 3 km fishnet over the polygon.

However, when I create the fishnet and check the attribute table of the new layer there is no FID column in the table.

Can anyone explain how to create a fishnet and retain the FID column?


texts: How to bind attributes to a Qt custom form?

texts:
I have followed this helpful post on creating a custom form. I have named the Qt objects the same as the fields in my attribute table. I have also set the edit widget properties to text edit for the line edit widget etc.

The form pops up when I click on a feature, but none of the data appears in the custom form.

Can anyone suggest what might be the problem?


texts: Projected Coordinate System issue using GDAL and qgis when defining x and y

texts:
I have a numpy.array which I would like to save as GeoTiff using GDAL in python. The array was created using the following PCS in mind:

  +proj=lcc +lat_1=48 +lat_2=48 +lat_0=48 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs

where lat_0=48 and lon_0=0 are the standard latitude and longitude where there should be no distortion. I looked at this guide but with no success. The problem is that when I define the rasterOrigin in degrees and define pixelWidth and pixelHeight in metres this does not seem to work and when I load the raster in qgis the raster is misplaced. See the example below taken from the guide I pointed out but tweaked.
The

  outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))

does not seem to work and the image is placed at the bottom left corner taken from the lat/lon details in the PCS and not from the SetGeoTransform which are different.

import gdal, ogr, os, osr
import numpy as np

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

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.SetProjCS( "myproj" )
    outRasterSRS.SetWellKnownGeogCS( "WGS84")
    outRasterSRS.SetLCC(53.5,53.5,53.5,0,0,0)
    outRasterSRS.SetLinearUnits("Meter", pixelWidth)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()

def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    reversed_arr = array[::-1] # reverse array so the tif looks like the array
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster

if __name__ == "__main__":
rasterOrigin = (-20,40)
pixelWidth = 10
pixelHeight = 10
newRasterfn = 'test.tif'
array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                  [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                  [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)


texts: Vector boundary generation -similar to Image Boundary

texts:
I want to automate the creation of boundary polygons from vector files. This is not footprint or extent creation (I already have scripts to do this). I want to create a tool similar to the QGIS, Image Boundary plugin which creates a single polygon which describes only the areas which have data. It should work on points/lines and polygons. For the points/lines, I guess it would need to create a buffer and then use that in generating a multipart boundary polygon (so it won't be exact but based on the buffer you set).

The purpose to use the output in a metadata catalogue where when you click on a location you only see a list of all files which have data in that spot and not just that the extent/footprint is within.

I found a tool from a post here that claims to do this but it is circa 2008 and isn't python based- http://blog.shoutis.org/2008/10/announcement-boundary-generator-v-01.html

An opensource QGIS/GDAL solution would be ideal but arcpy is fine as well.


texts: Extracting raster values at points using Open Source GIS?

texts:
How can I extract values from a raster by points?

I prefer not in Arcgis.

I prefer in Qgis or Mapwindow or other open source gis.


texts: Using Season instead of Month for Time Display in Time Manager?

texts:
Is there a way to make the season appear next to the year in Time Manager instead of the month? In the settings, I currently have DateTime selected with %Y-%m in the box. This makes the displayed time read 2017-10, for example.

I'm trying to make my frames into a GIF when I'm done and I'm hoping to not have to edit each individual frame in GIMP. Alternatively, if there's a way to edit the frames outside of QGIS in another program, I'm open to that, but I realize that a different stack exchange is dedicated to that kind of stuff.


texts:
MXD located in below path. I could able to open MXD from Windows explorer from dev box.
I am running the script in development box and accessing mxd from the server
isFile always fails.

\\server1\d$\arcgisserver\directories\arcgissystem\arcgisinput\Unknown_Unknown\poly_test.MapServer\extracted\v101\POLY_TEST.mxd

if os.path.isfile(mxdPath):
    print "Entered into successful"
    mxd = arcpy.mapping.MapDocument(mxdPath)
else:
    print "sorry"


texts: I wanted to know if there are any free available to download informative GIS/Geospatial audio or video podcast other than ESRI Instructional Series Podcasts  and Speaker Series Podcasts in the market. I am interested in general news/advancements in GIS.


texts:
I can't get my ArcPy script to delete a terrain. If it exists, it should be deleted. However it throws an

  ERROR 000317: {path} cannot be deleted

My script is a lot longer than this, but this is what it boils down to. fds is a feature dataset in a file geodatabase. The feature dataset contains a terrain with 2 associated feature classes:

a points feature class used as masspoints
a polygon feature class used as boundary/softclip

Here's my code:

import os
import arcpy

fds = r'D:\Dumfries & Borders.gdb\DEM'

arcpy.env.workspace = fds
l = arcpy.ListDatasets()
for d in l:
    if arcpy.Exists(os.path.join(fds, d)):
        arcpy.Delete_management(os.path.join(fds, d))

Earlier on I tried to have the script delete the whole feature dataset at once. Didn't work either. I then decided to delete the datasets in it one by one. The ListDatasets() lists only one: the terrain, which sparks the error mentioned above.

It is possible to delete the terrain or the feature dataset by hand in ArcCatalog, or using the Python window in ArcCatalog. However, even after closing all other programs (to preclude any locks from obstructing) except my Python interpreter (PyCharm) it still throws the same error.

So, any ideas?
Will post this to the ESRI forums as well.

edit follow-up:
removing the spaces and & from the gdb's name does not help, as ListDatasets() won't list anything now anymore. Also scriptwise trying to delete the feature dataset now throws an

  ERROR 000732: Input Data Element: Dataset {path} does not exist or is
  not supported.


texts:
In the sequence of my previous question (Can't use Raster Calculator with two different rasters), I am stuck again in my project.

With my source rasters, I created 3 different outputs using the following code:

calc100 = QgsRasterCalculator('aa@1*(1/28) + ao@1*(2/28) + am@1*(1/28) + ct@1*(6/28) + ln@1*(6/28) + ea@1*(2/28) + gl@1*(6/28) + gm@1*(2/28) + rs@1*(2/28)', 'D:/maps/z100.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100.processCalculation()

calc100to5000 = QgsRasterCalculator('ln@1 * (6/28+0.665*(d@1-100)/4900) + ct@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gl@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gm@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + rs@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ao@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ea@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + am@1 * (1/28-0.665*(d@1-100)/(4900*22)) + aa@1 * (1/28-0.665*(d@1-100)/(4900*22))', 'D:/maps/z100to5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100to5000.processCalculation()

calc5000 = QgsRasterCalculator('ln@1 * 0.879 + ct@1 * 0.033 + gl@1 * 0.033 + gm@1 * 0.011 + rs@1 * 0.011 + ao@1 * 0.011 + ea@1 * 0.011 + am@1 * 0.005 + aa@1 * 0.005', 'D:/maps/z5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc5000.processCalculation()

However, I just want to create one single output raster. I have one raster with the distance to shoreline at each point, and I was thinking if it's possible to insert an if statement to identify the value in that raster and use the respective formula to the pixel in question, something like this:

if d@1 < 100:
calc100 = QgsRasterCalculator('aa@1*(1/28) + ao@1*(2/28) + am@1*(1/28) + ct@1*(6/28) + ln@1*(6/28) + ea@1*(2/28) + gl@1*(6/28) + gm@1*(2/28) + rs@1*(2/28)', 'D:/maps/z100.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100.processCalculation()

if d@1 >= 100 and d@1 <= 5000:
calc100to5000 = QgsRasterCalculator('ln@1 * (6/28+0.665*(d@1-100)/4900) + ct@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gl@1 * (6/28-0.665*(d@1-100)*6/(4900*22)) + gm@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + rs@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ao@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + ea@1 * (2/28-0.665*(d@1-100)*2/(4900*22)) + am@1 * (1/28-0.665*(d@1-100)/(4900*22)) + aa@1 * (1/28-0.665*(d@1-100)/(4900*22))', 'D:/maps/z100to5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc100to5000.processCalculation()

if d@1 > 5000:
calc5000 = QgsRasterCalculator('ln@1 * 0.879 + ct@1 * 0.033 + gl@1 * 0.033 + gm@1 * 0.011 + rs@1 * 0.011 + ao@1 * 0.011 + ea@1 * 0.011 + am@1 * 0.005 + aa@1 * 0.005', 'D:/maps/z5000.tif', 'GTiff', aco_ant.extent(), aco_ant.width(), aco_ant.height(), entries)
calc5000.processCalculation()

If it's not possible, could you recommend a viable way?


texts: Python IDEs (PyScripter,PythonWin) crashing when installed with ArcGIS for Desktop?

texts:
I am learning Python, so here is a real newbie question!

Everything writen in ArcGis Python's windows is working fine. The same lines, correctly writen in an IDE, when run, are crashing it.

Working on ArcGis 10.2 using Python 2.7.5, 32 bits, Windows 7

IDE:

PythonWin32 (pywin32-218.win32-py2.7)
PyScripter (PyScripter-v2.5.3-Setup.exe)

For PyScripter I get a message Windows Error "python.exe is not working" followed by aPyScripter message: "EOFError: [Errno 10054] An existing connection was forcibly closed by the remote host" . Internet searches didn't get an answer I could understand.

For PythonWin, it just stops and crashs.

Installations path are as followed (by default):
Python 27: C:\

PythonWin: C:\Python27\ArcGIS10.2\Lib\site-packages\pythonwin

PyScripter: C:\Program Files (x86)\PyScripter

I suppose these issues have been met before and solve... So what I am missing here? some other packages or libraries? What would you advice me to use to script in Python?

For info my lines are:

import arcpy

from arcpy import env

env.workspace = "D:/_PYTHON/Python/Data/Exercise07"

fc = "Results/airports.shp"

delimfield = arcpy.AddFieldDelimiters(fc, "STATE")

cursor = arcpy.da.UpdateCursor(fc, ["STATE"], delimfield + " = 'AK'")

for row in cursor:

    row[0] = "AK2"

cursor.updateRow(row)

del row

del cursor


texts: Efficient clipping of polygons by other polygons using R

texts:
I have a dataset with overlapping polygons. I would like the clip the overlapping polygons (which I need to retain as overlapping) to a boundary (another polygon). I am using sf and st_intersection but due to the size of my data (60,000 polygons) this is not a quick process (I stopped it after an hour). For comparison I've tried using QGIS (I also stopped that after an hour with it still processing) and ArcMap (completed the task in 3 minutes). Is there a way I can complete the task using R in a reasonable time frame?
The code below recreates my data (I am after Points_Buffer being clipped to Boundary not Boundary_Clip):
library(sf)

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/England_gor_2011_clipped.zip",
              destfile = "England_gor_2011_clipped.zip" , mode='wb')
unzip("England_gor_2011_clipped.zip", exdir = ".")
file.remove("England_gor_2011_clipped.zip")

Boundary <- st_read("England_gor_2011_clipped.shp") %>% st_transform(., crs=27700)

Boundary_Crop <- st_crop(Boundary,st_bbox(c(xmin = 401000, xmax = 500000, ymax = 172000, ymin = 74000), crs = st_crs(27700)))

Points_Buffer <- st_sample(Boundary_Crop, 60000, exact = TRUE) %>% st_buffer(., 2000)

Points_Buffer_Clip <- st_intersection(Points_Buffer, Boundary)


texts: Extremely Laggy whenever I start georeferencing TIFF Image on ArcMap 10.6

texts:
I have been running my ArcMap on my new work computer for a month or so, it went fine. But since last week, ArcMap becomes slower whenever I load the map/georeferencing, these days when I georeference an image, it keeps loading and sometimes ArcMap will even stop responding for a while.

My daily tasks involve in making temporary selection from layers, therefore I have many temp layers (I'm not sure if that's the reasons why it lags)

I have tried to clear cache from a
ArcMap options and turn off snapping.


texts: how can I use OpenLayers2 Image layers as base layer?

texts:
I have created one code to display an image layer on my map using OpenLayers, I am using following code:

tiled = new OpenLayers.Layer.WMS(
                "layername", "domainname/geoserver/workspace/wms",
                {
                    LAYERS: 'layername',
                    STYLES: '',
                    format: 'image/png',
                    tiled: true,
                    transparent: "true",
                    tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
                },
                {
                    buffer: 0,
                    displayOutsideMaxExtent: true,
                    isBaseLayer: false,
                    yx : {'EPSG:900913' : true}
                }
            );
            map.addLayer(tiled);

My complete code is:

function init() {
    map = new OpenLayers.Map('mapdiv');
    var mapnik = new OpenLayers.Layer.OSM();
    var zoom = 15;
    map.addLayer(mapnik);
    map.setCenter(new OpenLayers.LonLat(151.1522771, -33.7288868).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:900913")), zoom);
    tiled = new OpenLayers.Layer.WMS("layername", "url", {
            LAYERS : 'layername',
            STYLES : '',
            format : 'image/png',
            tiled : true,
            transparent : "true",
            tilesOrigin : map.maxExtent.left + ',' + map.maxExtent.bottom
        }, {
            buffer : 0,
            displayOutsideMaxExtent : true,
            isBaseLayer : true,
            yx : {
                'EPSG:900913' : true
            }
        });
    map.addLayer(tiled);
}


texts: Why does the SRID of an ogr2ogr-imported MSSQLSpatial table change when using UID/PWD vs. trusted_connection?

texts:
I'm using python to run ogr2ogr commands to load SHP to MSSQLSpatial.

The command works fine when using trusted_connection, but when I use UID and PWD, the SRID is randomly changed to SRID:846 (while the data loads correctly, and the geometry previews correctly, it won't display in QGIS, which tipped me off to the problem)

Here are the commands, with line breaks are they are currently working in Python:

with Trusted Connection:

cmd = 'call "C:\Program Files\QGIS 2.18\OSGeo4W.bat" ' \
      'ogr2ogr -skipfailures -overwrite -progress -nln ' \
      '"plangis.temp_pyogr2ogr" ' \
      '-f MSSQLSpatial "MSSQL:driver={SQL Server};server=servername;database=databasename;trusted_connection=yes" "shp\statistical_neighborhoods.shp" -s_srs EPSG:4326 -t_srs EPSG:2877 -lco geom_name=shape -lco UPLOAD_GEOM_FORMAT=wkt'

with UID/PWD:

cmd2 = 'call "C:\Program Files\QGIS 2.18\OSGeo4W.bat" ' \
      'ogr2ogr -skipfailures -overwrite -progress -nln ' \
      '"plangis.temp_pyogr2ogr" ' \
      '-f MSSQLSpatial "MSSQL:driver={SQL Server};server=servername;database=databasename;UID=username;PWD=password" "shp\statistical_neighborhoods.shp" -s_srs EPSG:4326 -t_srs EPSG:2877 -lco geom_name=shape -lco UPLOAD_GEOM_FORMAT=wkt'

As you can see, both commands are identical except for the connection information.

When I print the assembled command, both are also identical except for the connection information.

I can run these at the OSGeo4W command, and the results are the same.

I'm testing the SRID output in MSSQL Server using shape.STSRID.


texts: R cannot read *.shp shapefile with rgdal

texts:
I'm trying to open ne_10m_urban_areas.shp with rgdal::readOGR, but it gives me the bellow error

library(rgdal)
a <- readOGR(dsn = path.expand("C:/.../Desktop/"),
layer = "ne_10m_urban_areas.shp")

  > Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  :

  > Cannot open layer

a <- shapefile("C:/.../Desktop/ne_10m_urban_areas.shp")

  > Error: file.exists(extension(x, ".shx")) is not TRUE

file.exists("C:/.../Desktop/ne_10m_urban_areas.shp")

  > [1] TRUE


texts: How to split multiple lines of the same feature using a line?

texts:
So obviously there is a split tool to split a single selected line. In my instance I would like to split 2+ lines of the same feature using a tool that is similar to the split tool.

In this image I have two road lines that are identical features. I want to split these lines using another line (the red line). I have tried a variety of different tools and the closest one to making it work is the


texts: Vector to Raster for Smoothing purposes using ArcGIS Desktop?

texts:
I have a vector file that represents transit scores by Census Tract (CT) of a city. It shows harsh boundaries where one CT has great access to transit and right beside that CT is one that has extremely poor access. I want to smooth these boundaries so as to blur these hard differences.

Is this appropriate? If so should I convert vector to raster, and resample to a larger grid size using bilinear method? I want to maintain the score value for future analysis.

Would the Smoothing Edge Zones (http://webhelp.esri.com/arcgisdesktop/9.3/index.cfm?TopicName=Smoothing%20zone%20edges) be a more appropriate method to use?


texts: How can I georeference a second image without losing georeferencing for first image in ArcGIS 10?

texts:
What happens:

1) Take a map with an image that was georeferenced earlier, which works fine.

2) Add a second image to the TOC and fire up the georeferencing toolbar.

3) Select the new, ungeoreferenced image from the drop-down list...

4) ... georeferencing on first image vanishes.

How do I georef the second image without losing the first georeferencing?


texts: PostGIS intersect query: Raster footprints of multiple high resolution GeoTIFF's or not?

texts:
I have multiple (381) high resolution GeoTIFF's (2x2 meter resolution, 1 band, 700 MB each) that I somehow want to upload to a PostGIS-server to be able to run ST_intersects or ST_3Dintersects queries to. I need the resolution preserved since every 2x2 meter square holds a height value in meters that new buildings can have before they are visible from a number of cultural heritage places. And to be extra clear the desired intersect operation I want to run is in the z-dimension, not only based on the coordinates. Ideally I would like to be able to create a polygon shapefile with a z-attribute for a planned building and run a st_intersect query to see from which of the 381 points the new building would be visible from.

When asking about this at Stackexchange before, the question got closed because it was to broad(which it indeed was). One user questioned the need to load the raster's to postgis when the footprint of the image was all I needed. Since then I've been trying to understand what a raster footprint is, but I'm not sure I've figured it out yet:

I've tried to Polygonize a raster with QGIS (Raster>Conversion>Polygonize) which worked and stores the z-values I want to use BUT the vector version of the 700 MB GeoTIFF takes up 4,33 GB of storage. So I'm guessing this approach is even worse then uploading the raster straight to PostGIS with raster2pgsql? So I'm guessing this is not what the user had in mind when suggesting footprints as a way forward?

I've also tried a plugin for QGIS called Image Footprint which from the plugin description sound like just the plugin I need:

  This plugin create a catalog layer from directories of images, how a memory layer, where each feature has metadata and footprint of a image.
  See presentation in http://pt.slideshare.net/LuizMotta3/qgis-ibama-imagefootprint

The plugin gives me two options: 1. Bounding Box or, 2. Valid Pixels.

The bounding box shapefile is only 176 KB big which is handy, has two columns "filename" and "meta_html", but no values from the raster is present in the shapefile which only contains one feature and some meta data about the area of valid pixels, number of holes, projection etc.
The Valid pixels shapefile is also 176 KB and has almost the same meta data except for more info about pixel size and extent.

But none of the above results from this plugin gives me any actual info from the 2x2 meter cells in the original raster. Although the attributes contain a path to where the actual raster is stored on the hard drive. So this can not be the way forward either?

My question is how would you store these high resolution GeoTIFF's on a PostGIS-server to be able to run the above mentioned st_intersects query with the fastest possible performance?

Server runs Ubuntu 17.10 with Postgresql 9.6.6/PostGIS 2.4.1 on a Intel i5-6600k CPU with 4 cores and 4 logical processors, 16 GB om physical RAM and a 1 TB SSD harddrive.

Metadata from one of the rasters:


texts:
I would like to take geometries from a vector dataset and reduce them to a hash.  This hash would then be used to verify the integrity of that data and also identify identical geometries.

Do any appropriate algorithms exist that could be used?  What pitfalls could I encounter?


texts: Calculating accumulated walking pace from starting location outwards using R?

texts:
Goal

In R, I would like to calculate an accumulated cost surface around a starting location, with the cost representing the time (say, hours) to walk from the location outwards in all directions (say, 8 directions in raster's terms). The walking time (pace) is devised via the Tobler's hiking function (https://en.wikipedia.org/wiki/Tobler%27s_hiking_function).

I have indeed used ArcGIS to achieve that in the past, using the 'Path Distance' tool (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-path-distance-tools-work.htm), which implements the method devised by Cixiang Zhan, Sudhakar Menon, Peng Gao (https://link.springer.com/chapter/10.1007/3-540-57207-4_29). I would like to reproduce something similar in R.

Earlier threads

I have made a search in this Forum, but at the best of my knowledge the issue has been always tackled from an ArcGIS perspective.

Issue

In theory, the implementation of the methodology is simple, but I have hard time to devise the correct workflow in R. Also, I do believe that the gdistance package is the key to the calculation of an accumulated cost surface, but I cannot wrap my head around how to properly use the transition() function.

What I did so far

# Define the Tober's function for speed in kmh
tobler.kmh <- function(x){ 6 * exp(-3.5 * abs(tan(x*pi/180) + 0.05)) }

# calculate the slope from an input DTM
SLOPE <- raster::terrain(dtm, opt='slope', unit='degrees', neighbors=8)

# calculate the speed in Kmh
v.kmh <- calc(SLOPE, tobler.kmh)

# time in h to cross 1 km (p stands for 'pace')
p.hkm <- 1/v.kmh

# express the speed in meters per hour
v.mh <- v.kmh * 1000

# time in hours to cross 1 m
p.hm <- 1 / v.mh

Where I am stuck

Now, with the ultimate goal to devise isochrones, I think I would need to accumulate the pace (in hours per meters = p.hm) outwards from a given location. I guess that the gdistance's accCost() function would make the trick, but I cannot understand how to tackle the intermediate step that entails creating a transition matrix via the transition() function. The latter (at the best of my understanding) would allow to define a traversing cost other than distance in which one incurrs when moving from one cell to an adjacent cell.

I do not know: (a) if the above formula to calculate the pace (in hours per meters) has to be somehow implemented withing the transition() function; (b) if I am approacing the whole issue from the wrong perspective.


texts: How to connect and access Mapd Database by GeoServer

texts:
In my project requirement to access table data from Mapd Database on GeoServer

Mapd(GPU Based database) Database documentation link:
https://www.mapd.com/docs/latest/mapd-core-guide/users-and-databases/

But in GeoServer provided some extension to connect and access particular database

GeoServer Extension link:

http://docs.geoserver.org/latest/en/user/data/database/index.html#data-database

How can I connect MapD databse with geoserver is there any plugin.


texts:
First, let's set up the data. Download adjoining Seoul, Incheon, Gyeonggi-do region shapefiles from this website (in Korean). I attach the files, as they are not accessible outside Korea. Then read into R, name, and combine (credit):

library(sf)

#Read, name and combine regions
sca_nsdi <- rbind(cbind(st_read("LARD_ADM_SECT_SGG_11.shp"), Name="Seoul"),
                  cbind(st_read("LARD_ADM_SECT_SGG_28.shp"), Name="Incheon"),
                  cbind(st_read("LARD_ADM_SECT_SGG_41.shp"), Name="Gyeonggi-do"))

#Plotting reveals sub-regions (the legend was natively clipped)
plot(sca_nsdi[, "Name"])

I wish to merge the sub-regions of these 3 regions. I have tried applying this solution, but without success:

library(dplyr)

sca_nsdi <- sca_nsdi %>% group_by(Name) %>% summarise(geometry=st_union(geometry)) %>%
  ungroup()

plot(sca_nsdi[, "Name"])

The sub-regions are imperfectly merged. (For some reason, the bounding box is also shrunk to the size of Seoul.)


texts: Modelbuilder: nested iterators for intersect

texts:
I have two geodatabases: one with line layers (pipelines), and one with polygon layers (stripmap). I would like to use model builder to take the first line layer and intersect it with the first polygon layer. Then take the second line layer and intersect it with the second polygon layer, and so on. In my testing each geodatabase has 3 layers.

I have used two iterators to do this, by creating a submodel (see image), but when I run it, the first two output shapeiles do not contain any features. The final layer works correctly. I have run intersect independently on the first line layer and first polygon layer and it does generate a shapefile containing the overlaps.

I have read through similar questions about nesting iterators, but can't seem to find the answer.


texts: Openlayers + WMS layer. Use scale or resolution + DPI

texts:
I have a problem with displaying map with WMS layer in Openlayers. I'm able to display it but the image of the map is not sharp especially at higher zooms.

I have a documentation from the map provider and below is part of description:
The data scale set is:

  1:20 mil., 1:6 200 000, 1:3 000 000, 1:1mil., 1:500 000, 1:200 000,
  1:75 000, 1:25 000, 1:12 500 , 1:6 250

For the best quality of displayed data while zooming is recommended to count with recalculation
96dpi/200dpi (Map DPI = 200), it means 0,48 x default scale set.

I'm not sure if this DPI setting that is wrong in my code. I found out that Openlayers standard DPI is 72 DPI.

Can you suggest anything how can I fix it?

WMS getcapabilities:

<SRS>EPSG:32633</SRS>
<LatLonBoundingBox minx="-8.060897" miny="32.143274" maxx="56.352488" maxy="66.324824"/>
<BoundingBox SRS="EPSG:32633" minx="-1700000" miny="3800000" maxx="2300000" maxy="8000000"/>

My code:

<!DOCTYPE html>
<html>
  <head>
    <title>WMS map</title>
    <script src="http://openlayers.org/api/OpenLayers.js"></script>

<script type='text/javascript'>

var map; function init() {

{
   var mapoptions = {
      maxExtent: new OpenLayers.Bounds(-1700000, 3800000, 2300000, 8000000),
      numZoomLevels:14,
      scale: [20000000, 6200000, 3000000, 1000000, 500000, 200000, 75000,
      25000, 12500, 6250],
      maxResolution: "auto",
      projection: "EPSG:32633",
      displayProjection: new OpenLayers.Projection("EPSG:32633")

    }};
    map = new OpenLayers.Map( 'map',mapoptions );

   var wms = new OpenLayers.Layer.WMS( "layer_name" ,'https://wms.dot.com',
   {layers: 'layer_name',
   format: "image/png"},
   {}
);
map.addLayer(wms);
if(!map.getCenter()){ map.zoomToMaxExtent(); } }
</script>
</head>

<body onload='init();'>
    <div id='map' style='width: 1000px; height: 1000px;'>
</div>
</body>
</html>


texts: Uploading maps made using ArcGIS for Desktop to Google Maps?

texts:
I've made a few maps using ArcGIS for Desktop. I wanted to know if it was possible to upload these maps onto Google Maps.

Is there a way to do this?


texts: Raster incorrectly reprojected to OSGB(27700)

texts:
Ubuntu 12.04, QGIS 1.8.0, GDAL/OGR 1.7.3, GEOS 3.2.2

Open a Spatialite db with OS Mastermap poly layer (EPSG:27700);
Add OpenLayers Plugin (1.0) from QGIS offical repository;
Add a Google Satellite layer (EPSG:3857) and make OS poly fills transparent;
Zoom & drag map a bit till it 'settles down' - then both aligned nicely;
Make poly layer not visible and save map as image (jpg + jpgw);
Add saved jpg back in as a layer (EPSG:3857) - still aligned;
Reproject (warp) jpg layer from EPSG:3857 to a EPSG:27700 Gtiff;
Add this Gtiff layer - no longer aligned but shifted about 6m North.
???!!!???

(On the fly CRS transformation is enabled)

I've had this issue before with the unstable flavour of UbuntuGIS QGIS, so went back to the stable & all was ok. This time I'm getting the 'shift' in the stable UbuntuGIS QGIS and no amount of wiping out the libraries (with Synaptic) & reinstalling has so far got rid of this 'shift'.

Any help or ideas gratefully received...

No problem reproducing the error (thanks underdark), but some more digging might narrow the issue a bit more.

The proj4 CRS projection string for OSGB/27700 in GDAL1.9.1 is:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=375,-111,431,0,0,0,0 +units=m +no_defs

The same thing in GDAL1.7.3 is:

+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs

Even though QGIS is reporting the CRS as 27700 with the same proj4 string as the earlier GDAL1.7.3 it's probably using the later GDAL1.9.1 parameters for the reprojection behind the scenes. That would be fine if the image lined up after reprojection - but it doesn't.

To test a step further I created a QGIS user-defined CRS with the same projection parameters as the earlier GDAL version, and then reprojected with this - everything now lines up fine.

Are the EPSG:27700 projection parameters the likely suspect in the later GDAL (which QGIS is using) - or am I barking up the wrong tree?


texts: Creating an arc with set angle and range using QGIS?

texts:
I'm currently using QGIS 3.6.

For a while now to map the arcs of fire for historical fixed weapon emplacements. Creating an arc to a set range to a set angle is extremely important for my analysis work.

I currently use the buffer tool to create a circular polygon around a point to the desired range and direction. I then have to essentially delete points until the desired angle is achieved. This isn't that precise, however as I have to estimate the correct angle. Doing this for every point is very time consuming and I was wondering if there was a quick and simple way to achieve the same result in QGIS, without manually setting the angle by deleting points and without plugins (if possible).

I am not a technical person and have no programming experience.

See below for an example of what I am hoping to achieve.


texts: Avoiding Raster to Shapefile (Polygon) errors 000863 and 000308?

texts:
Have you experienced converting a polygon from raster? The raster dataset is a satellite image created from zonal statistics? The problem is when I was converting the raster to polygon errors came out as:

  ERROR 000863: The input is not within the defined domain.
  ERROR 000308: Invalid field type

I do not know how to get around these errors.


texts: Clip Raster by Mask Layer to jpg not tiff

texts:
Is it possible to set your output file for the Clip Raster by Mask Layer tool as a jpg and not a tiff?  Are there extra parameters or commands I can use to do this?

I can export as tiff and convert to jpg with another program but I am doing this for a large amount of features and it's extremely inefficient.


texts: Merging overlapping buildings but not adjacent buildings using ArcMap?

texts:
I have a dataset that has buildings generated from satellites which leads to imperfect building outlines. The buildings circled in red are one building but the computer drew them as having multiple pieces.

How can I merge them into one building while keeping the other buildings that are adjacent but not overlapping from merging. I tried every union, merge, dissolve, and nothing was able to distinguish between overlapping and adjacent. Or I wasn't able to make it distinguish between the two. I am using ArcMap.


texts: Custom error message with 'Error Code -1' in ArcGIS script tool

texts:
I have a script tool in ArcGIS 10.2 and I have altered the tool's validation code to show my custom error messages using setErrorMessage method.
It works but when an error message is displayed it shows error code -1 (ERROR -1) together with the message text, and since it doesn't link to anything I don't want to show it at all.

Custom message looks like this:

I have been looking for this in the documentation but I can't find it. How can I hide it or prevent it to show?

Here's my validation code:

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

    if self.params[0].value:
      if not arcpy.Exists(self.params[0].value):
        self.params[0].setErrorMessage("Invalid geodatabase: path typed incorrectly or it does not exist.")
      else:
        desc1 = arcpy.Describe(self.params[0].value)
        if desc1.workspaceFactoryProgID != 'esriDataSourcesGDB.AccessWorkspaceFactory.1':
          self.params[0].setErrorMessage("This workspace is not a personal geodatabase")

    if self.params[1].value:
      if not arcpy.Exists(self.params[1].value):
        self.params[1].setErrorMessage("Invalid database connection: path typed incorrectly or it does not exist.")
      else:
        desc2 = arcpy.Describe(self.params[1].value)
        if desc2.workspaceFactoryProgID != 'esriDataSourcesGDB.SdeWorkspaceFactory.1':
          self.params[1].setErrorMessage("This workspace is not an SDE database connection.")

    return


texts: How to merge two point feature classes based on spatial location using ArcPy?

texts:
I have two point feature classes, each containing about fifty points.
These 2 point feature classes should be merged, or the points of feature class 2 added to feature class 1, as feature class 2 does not contain any information, but only the points itself.

But if the distance between a point of feature class 1 and feature class 2 is smaller than 50 meters, they should be merged to one point.

The figure should illustrate the problem

How can I do this programmatically (with an ArcPy script)?

One way I could think of is:

run buffer around feature class 1,
delete points from feature class 2 that fall into the buffer,
add remaining points of feature class 2 to feature class 1

How can I achieve this?


texts: Calculating albedo from Landsat in ArcGIS Desktop?

texts: I need to calculate albedo from Landsat images using ArcGIS Desktop. I am fairly new in RS and GIS and I know that I have to change DNs to TOA and surface reflectance but I cannot find a procedure on how to do this in ArcGIS Desktop.


texts:
Instead of two shapefiles that are saved in fcList in this loop I would like to do KernelDensity on the unique values of one shapefile.

Does anyone have an idea how this might work?

import arcpy,os
from arcpy import env
env.workspace = r'C:\Users\lhark\Desktop\SHK\ArcGis_Proj_Canada\shapefiles'
outworkspace = r'C:\Users\lhark\Desktop\SHK\ArcGis_Proj_Canada\KernelDensity'

fcList = arcpy.ListFeatureClasses ()
for fc in fcList:
outKDens = KernelDensity(fc,"NONE",4670)
outKDens.save(outworkspace + "/" + os.path.splitext(fc)[0] + "_kd.tif"


texts: Publishing linear referencing capability in ArcGIS Server?

texts:
I created a linear referencing event and want to publish as linear referencing service in ArcGIS 10.2 for Server . The problem is when i want to publish, High warning shows the linear referencing event table not recognized .in other hand the linear referencing item in capabilities not found .

How can i publish the linear referencing service?


texts: Difference between TIFF files exported by ArcGIS and GDAL

texts: I am doing some testing were I have created a radom raster in ArcGIS in TIFF format. It created a .tif and a .tfw file (where it stores the coordinates). Then I imported this TIFF file into GDAL, did some operations as a numpy array with it and then exported again the array to TIFF file with GDAL. The result is a .tif file, but no .tfw file. However, this resulting TIFF file has all the right coordinates and CRS defined in ArcGIS. How come the one exported from GDAL doesn't need a .tfw file and the one exported from ArcGIS does? Are they different types?


texts:
I am interested to demarcate the boundary of a specific city and then work on the data collected from that city in QGIS. But I am unable to specify boundaries of the required city from the vast google satellite imagery.
Can you please tell me how can I demarcate a whole specific city in QGIS?


texts: Flex Geoprocessing Widget output parameter invalid

texts: I want to publish a ArcToolbox with a script on a Flex application as a widget. I have ArcGIS server from which I will publish it as a geoprocessing service. The tool has an input parameter which is a text file and the output parameter has to be a folder which will route the user to choose a folder on computer. When I try to publish this tool, it says that the output parameter - folder is invalid. Is it possible that the folder parameter is not supported by flex output data? Any way to do that?


texts: Counting points in polygons using ArcGIS Desktop?

texts:
I have a table as follows:

I have points and polygons. Every polygon has a variety of points inside. Some have 6, some have 12, some have 0.
FID_2* is polygon ID.
I would like to populate Count( 9th column above) automatically according to FID_2*, that is, polygon ID.
For example, when FID_2* says 46( first row) I want Count to be populated as 6 because I would like to match each point with the polygon they are in. (Polygon 46 has 6 points)

This is another picture showing my points and polygons:

Could you help me solve this?


texts: How can I export raster data from Arcmap to an R  *.grd file format?

texts:
I want to work with raster data in R, but I have problems exporting my raster dataset from ArcMap. By "export data"/GRID I get a file with an *.aux extension. What I need is *.grd.

Can anybody help me out?


texts: Configuring proxy.cgi in Windows using Apache and GeoServer?

texts:
I want to dispaly my layers on a server. I'am using geoserver 2.3.2 and apache 2.4

Till now i Create a proxy file in referring in this link : http://trac.osgeo.org/openlayers/browser/trunk/openlayers/examples/proxy.cgi  and i Put it in cgi-bin folder.  I Installed python and in the starting line of proxy.cgi i put the location of python.  #!C:/Python27/python.exe -u

I found in some tutorials that i should configure the file httpd.conf  and i have no idea how to configure it and where should i add this code :

" ProxyRequests Off
ProxyPreserveHost On

Order deny,allow

Allow from all

ProxyPass /geoserver //localhost:8080/geoserver

ProxyPassReverse /geoserver //localhost:8080/geoserver "


texts: Export shp geometry into Wkb format using PyQGIS

texts:
I try to get features geometry from a polygon layer. I will insert this geometry into one of my postgis table so I think it would be easier if I export the geometry into WKB format. But sadly I just recognise that QgsGeometry doesn't support function to export geometry into wkb format (it's just for wkt and geojson).

Here is my script, but it's using exportToWkt function:
from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon
from qgis.core import QgsMapLayer
from qgis.core import QgsFeature, QgsGeometry
selectedLayer =iface.mapCanvas().layers()[0]
for feature in selectedLayer.getFeatures():
    geom = feature.geometry()
    print geom.exportToWkt()

Any ideas to export this geometry into wkb format ?


texts:
I want to edit object on WFS layer but have one problem. I use OpenLayers.Control.ModifyFeature(myVecLayer) to editing.By the logic i do this.

Press button and create WFS layer with chosen object
Add OpenLayers.Control.ModifyFeature to map controls  and activate
this
Edit object
Press button again and remove WFS layer with
my object, deactivate OpenLayers.Control.ModifyFeature control.
Select next object (or same object) and do "1", "2"

And object not editing.
But if I do not remove layer after editing, but just add new WFS layer to map OpenLayers.Control.ModifyFeature control works

My code:

function ff(){
if(counter==0){
    //WFS
        var myStyle = OpenLayers.Util.extend({},
                                OpenLayers.Feature.Vector.style['default']);
                            myStyle.strokeWidth = 1.5;
                            myStyle.strokeColor = "#ff0000";
                            myStyle.fillOpacity = 0.1;
                            OpenLayers.ProxyHost = "proxy.cgi?url=";
                            var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
                            renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;
                            myVecLayer = new OpenLayers.Layer.Vector("Редактируемый участок", {
                                strategies: [new OpenLayers.Strategy.BBOX(), new OpenLayers.Strategy.Save()],
                                style: myStyle,
                                protocol: new OpenLayers.Protocol.WFS({
                                    version: "1.0.0",
                                    srsName:"EPSG:900913",
                                    url: "http://localhost:8080/geoserver/wfs",
                                    featureType: "filedata",
                                    featureNS: "http://www.opengeospatial.net/cite",
                                    defaultFilter : new OpenLayers.Filter.Comparison({
                                        type: OpenLayers.Filter.Comparison.EQUAL_TO,
                                        property: 'id',
                                        value: "162",
                                    })
                                }),
                                renderers: renderer
                            })
                            map.addLayers([myVecLayer]);

                            if(ppp==null){
                                if (edit==null){
                                    edit = new OpenLayers.Control.ModifyFeature(myVecLayer);
                                    map.addControl(edit);
                                    ppp=map.getControlsByClass("OpenLayers.Control.ModifyFeature");
                                }

                                alert(ppp);
                            }
                        //  ppp.autoActivate=true;
                            for(var i=0;i<ppp.length;i++){
                                ppp[i].activate();
                                alert("2 "+ppp[i].active);
                            }

                            counter=1;

                }

Its will nice if u help.


texts: Distance from rivers, taking into account terrain elevation

texts:
I am interested in building a raster of the distance from rivers in a territory under study. I have digitized the course of the rivers in ArcGis 10.1.
I have used the 'Euclidean Distance' tool to create a raster of the distance from rivers. Now, since we are moving in a 3D world, not in a 2D one, I would like to 'weight' the calculated distances in order to take the elevation into account.

So, I would like to know how can I create a raster depicting the distance of any point from rivers, also taking into account the elevation variable.


texts: Intersection of two shapely polygon geometries

texts:
I have Shapely polygon geometries, A circle(c) and list of polygons within the circle. I want to find the union of the circle with the intersections. (circle + intersections)

I was trying to make a MultiPolygon of difference and intersection:

poly_list = [poly_N, poly_E]

# finding the intersection of circle and polys
intersec_list = []
for poly in poly_list:
    intersec_list.append(circle.intersection(poly))

MP_intersec = MultiPolygon(intersec_list)

# finding the difference of circle and polygons
diff_list = []
for inter in intersec_list:
    diffi = c.difference(inter)

MP_diff = MultiPolygon(diff_list)

# find the union of circle and intersection
# output: a MultiPolygon or a geometric collection with 4 polygons and the difference of circle.

MP_out = MultiPolygon([MP_diff, MP_intersec])

but MP_out gives invalid geometry probably because I have overlapping polygon as well.


texts: Do Sentinel-2 L1C data downloaded from USGS Earth Explorer consist of reflectance values or digital numbers?

texts:
I was wondering if Sentinel-2 L1C data downloaded from USGS Earth Explorer consist of reflectance values.
And is it necessary to use reflectance values for time series analysis of Sentinel-2 especially in yield determination and phenology extraction?
I created some vegetation indexes (LAI, NDVI, EVI) from Sentinel-2 L1C data after directly downloaded it without any processing. How accurate is what i did? Should i have done something different?


texts: Error attempting to remove layer using arcpy in ArcGIS Pro?

texts:
I am trying to remove layers from an ArcGIS Pro 1.1 project, but I receive a ValueError with the traceback of

File "C:\Program Files\ArcGIS\Pro\Resources\ArcPy\arcpy\_mp.py", line 1280, in removeLayer
    return convertArcObjectToPythonObject(self._arc_object.removeLayer(*gp_fixargs((remove_layer,), True)))

Here is code to add raster layers from a geodatabase

rasters = arcpy.ListRasters("*",None)
aprx = arcpy.mp.ArcGISProject("Path/to/aprx")
aprxMap = aprx.listMaps("Map")[0]
for raster in rasters:
     raster = env.workspace + "\\" + raster
     aprxMap.addDataFromPath(raster)
     aprx.save()

This successfully adds the raster objects as layers to the current map, but when I try to remove using the following code, I get an error.

for raster in rasters:
         aprxMap.removeLayer(raster)
         aprx.save()

Any ideas why this is happening?

Raster is a string, and when I have tried to write as an input string, I still receive the error.


texts: Determining how interactive web map was created?

texts:
I found this on github, http://mbostock.github.io/d3/talk/20111116/airports.html, with no map credits.

Does anyone know how this interactive web map was created?

I need to make a map in this style and cannot find out which program to use or if it requires coding to create.


texts: Selecting from Point column in PostgreSQL?

texts:
How to select from Point column data type(native) in PostgreSQL?

select a.* from indsolv.address As a where a.location = '93.2321, 21.0321'

Getting error for the above query!

ERROR:  operator does not exist: point = unknown
LINE 1: ...ct a.* from indsolv.address As a where a.location = '93.2321...
                                                         ^
HINT:  No operator matches the given name and argument type(s). You might need to add explicit type casts.
SQL state: 42883
Character: 55


texts: Channel Easements for County: Need help with Query or Python

texts:
I have a layer for easements for a whole county and need to find each easement either MISSING or having a PART of what the county calls "STW_ARCS." These ARCS are what the county have data on for how much linear footage in each (or not each) easement.
.
                   I'm having a very hard time using tools to find these easements as described above.