texts: All WMS layers not showing in OpenLayers?

texts:
I'm trying to display information on OpenLayers from different sources, I guess the problem is that some of these layers has information for the same areas as another layer already has. And OpenLayers doesn't allow overlapping of layers?

I found something about setting "allOverlays=true" should do the trick, but nothing changes. When display the layers one by one everything works fine. I just can't combine them.

Here's my code:

var layers = [

    new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: 'http://localhost:80/WMSServer1',
        params: {'LAYERS': 'buildings,restrictedarea'},
        ratio: 1
      })
    }),

            new ol.layer.Tile({
      source: new ol.source.OSM()
    }),
    new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: 'http://localhost:80/WMSServer2',
        params: {'LAYERS': 'areas,parks,buildings'},
        ratio: 1
      })
    }),

    new ol.layer.Image({
      source: new ol.source.ImageWMS({
        url: 'http://localhost:80/WMSServer3',
        params: {'LAYERS': 'areas'},
        ratio: 1
      })
    })

  ];
  var userPos = "";
    if ("geolocation" in navigator) {
      navigator.geolocation.getCurrentPosition(function(position) {
          do_something(position.coords.latitude, position.coords.longitude);
        });
    } else {
      /* geolocation IS NOT available */
    }

  var map = new ol.Map({
    layers: layers,
    target: 'map',
    allOverlays:true,
    view: new ol.View({
      center: [-10997148, 4569099],
      zoom: 4
    })
  });

  var lonLat = new ol.LonLat( -0.1279688 ,51.5077286 )
      .transform(
        new ol.Projection("EPSG:4326"), // transform from WGS 1984
        map.getProjectionObject() // to Spherical Mercator Projection
      );

var zoom=16;
var markers = new ol.Layer.Markers( "Markers" );
map.addLayer(markers);
markers.addMarker(new ol.Marker(lonLat));
map.setCenter (lonLat, zoom);

    if ("geolocation" in navigator) {
      navigator.geolocation.getCurrentPosition(function(position) {
          do_something(position.coords.latitude, position.coords.longitude);
        });
    } else {
      /* geolocation IS NOT available */
    }
    function do_something(latUser,longUser){
        var newCoord = ol.proj.transform([longUser, latUser], 'EPSG:4326', 'EPSG:3857');
        console.log(newCoord);
        map.setView(new ol.View({
                center: newCoord,
                zoom: 10
         }));

    }


texts: How to extend polylines to the nearest point

texts:
Could anyone tell me a way of getting my polylines to search for a point within a certain number of meters of their start/end points and extend to it?

There was some problem uploading a picture so I'll explain further. I have contours with gaps like here and am able to close them, but have found that other gaps, caused by edges to the area and faults are also being filled, the lines also travel completely straight when this is not quite the case. I also have a point feature class of where annotations used to be. The gives one more point when the gap should be slightly curved, and identifies the gaps that I actually want to fill. So I want to do something similar to closeing rings, but instead of extending the lines to meet eachother I want to extend them to meet the point.


texts:
I set up a script to print some statistics about line segments. The "logfile" is basically the same as "print", but to a textfile instead of on screen only.

  LogFile("Feature {0} has a length of {1} meters at coordinates
  {2}".format(row[0], row[1], row[2]))

This will give me an output like:

  Feature 1 has a length of 346.239783695 meters at coordinates
  (470765.69285000017, 6313169.6829)

I want to drop the decimals (i.e., turn the values into integers for reporting purposes only) and just have

  Feature 1 has a length of 346 meters at coordinates (470765, 6313169)

How do I alter the format within the script to accomplish that? I assume it will be super easy and something I should know, but I've been out of the Arc/Python world for a few years.


texts:
Are there any known hacks to get subsurface data into Google Earth?

Sub-ocean-surface works with the release of Google Ocean, but you cannot go below the sea floor.

The other major limitations of Google Earth (no vertical data, single placemark balloon) have hacks to get around them (thin buildings, custom icons). I'm hoping someone has a subsurface hack.

Edit:

Possible approaches might be

Find the default image and replace it with a transparent PNG
Capture the network requests to download tiles and respond with transparent PNGs
Convince GE to set the Level of Detail (LOD) for the primary layer to some negative or infinite value.

On OS X, I think the following folder/code looks promising:

cd Google Earth.app/Contents/MacOS/shaders
$ ls atmosphere_ground_sun_off*
atmosphere_ground_sun_off.arbfp1
atmosphere_ground_sun_off.asd
atmosphere_ground_sun_off.ps_2_0
atmosphere_ground_sun_off.arbvp1
atmosphere_ground_sun_off.cfg
atmosphere_ground_sun_off.vs_2_0

Those files appear to contain assembly code instructions to the graphics card, and there is mention of "groundTexture", "groundRayleighMap", and other possible hack-able things.


texts:
What tool does:

Simply creates tbf file from raster extent.

What is the problem:

I ran toolbox in ArcMap 10.5.1 and ArcGIS Pro 2.2.4.

When I run it as standalone script it is working just fine.

I am not new to Python nor Python Toolboxes, but today:

1) I received an error which I never encountered before (error shown below)

2) when I declare rasters = arcpy.ListRasters('tif'), list of rasters is empty.

Error code:

  Running script Tool...

  Traceback (most recent call last):
    File "", line 54, in execute
    File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing_base.py", line 527, in set_
      self[env] = val
    File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing_base.py", line 587, in setitem
      ret_ = setattr(self._gp, item, value)
      enter code here

  RuntimeError: Object: Error in accessing environment
   Failed to execute (Tool).

Tool code below:

import arcpy
import os

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

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

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

    def getParameterInfo(self):
        """Define parameter definitions"""
        # Pierwszy parametr
        inside = arcpy.Parameter(
            displayName="In workspace",
            name="in_gdb",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input")

        p = [inside]
        return p

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

    def updateParameters(self, parameters):
        """Modify the values and properties of parameters before internal
        validation is performed.  This method is called whenever a parameter
        has been changed."""
        return

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

    def execute(self, parameters, messages):
        """The source code of the tool."""
        arcpy.env.overwriteOutput = True

        arcpy.env.workspace = parameters[0].valueAsText
        output = parameters[0].valueAsText
        rasters = arcpy.ListRasters('*.tif')

        i=1
        for ras in rasters:

            arcpy.AddMessage('Working on raster nr:  {0} of {1} '.format(i, len(rasters)))
            f = arcpy.Raster(ras).extent

            north = f.YMax
            west = f.XMin
            south = f.YMin
            east = f.XMax

            i+= 1
            ras_split=ras.split('.')[0]

            with open(os.path.join(output,'{0}.tfb'.format(ras_split)), "w+") as text_file:
                text_file.write('north={0}\n'.format(north))
                text_file.write('west={0}\n'.format(west))
                text_file.write('south={0}\n'.format(south))
                text_file.write('east={0}'.format(east))

        return

Fixed all mistakes - tool fully working.


texts: Add/display all images of mycollection in google earth engine

texts:
I started with GoogleEngine recently and I was wondering how I can display all images of my collection as layers. I guess there is a better way than manually copy past / typing:

 var img1 = ee.Image('COPERNICUS/.....);
 var img2 = ....;
 var imgx = ....

 Map.addLayer (img1);
 Map.addLayer(img2);
 Map.addLayer(imgx)....

This questions was already asked in 2017 but I was unable to retrace the link/shared google group...


texts: Pan and Zoom via mouse not working in New Data Frame?

texts:
Inserted a new data frame for a small overview map in the corner of my display. Left click and selected properties, set the same coordinate system as the larger display, set the Extent Indicators to the Layers data frame, and left the extents as Automatic.

I select a GIS basemap as the background in the new data frame and can see the little red box showing the location of the larger display. My problem is, I cannot use the mouse or mouse wheel to pan or zoom in this New Data Frame box. I can set the scale in the scale bar and also use the magnifying-glass buttons but not my mouse.

Suggestion commonly found on google searches: make sure I have the Tools checked and not just Layout Tools. I do. Left click on the gray space near the toolbars, I have Tools and Layout Tools checked.

A lot of people seem to have their movement buttons grayed out. I do not! None of the zoom/pan buttons are ever grayed out. I can click them and my mouse arrow changes into the little hand, I just cannot pan or zoom using my mouse or mouse wheel if I double click into the New Data Frame, I have to use the magnifying glasses.

I did find a work around but I don't know how it's working: when I save and close out of the MXD and then open it up again, I can double click into the New Data Frame and the pan and zoom respond to my mouse and mouse wheel.

I am using ArcMAP 10.5.1 for desktop. This has been a nuisance but since I know a work around it is not impeding my work....I am just hoping to finally get an answer once and for all.

Edit: Here is a step-by-step of what I am doing and where I think the program is not working as it should.

Using ArcMap 10.5.1 for desktop
Open new MXD and save to local location with simple name
Set coordinate system and add GIS basemap layer
zoom in to desired location and set scale (panning and zooming with mouse and mouse wheel works fine here)
check that both tools and layout tools are checked
insert new data frame and set active
set to same coordinate system and add GIS basemap
Open properties and set extent indicators to layers, close properties and little red square is visible in new data frame
double click inside new data frame. Edge of data frame is striped grey showing that I am inside of it but mouse and mouse wheel still move/zoom the layout view, not the data frame. This is where I think the program is not working correctly and if I want to adjust the scale of my new data frame I have to use the magnifying glass tools + the scale bar.
However, if I save and exit and immediately open the same MXD, I can double click into the new data frame, select the pan button (arrow changes to the little hand), and use the mouse/mouse wheel to pan and zoom. A good work around but I feel like the program should allow me to do this whenever I first insert a New Data Frame.
Please not that I have not added any features or shapefiles to this example MXD but I have had this problem with inserting new data frames in MXD's that also have callouts, shapefiles, etc displayed.


texts: Expression String Builder window always opens when I start QGIS

texts: When I open QGIS, the Expression String Builder window always opens and I have to cancel it in order for my project to open. How can I stop this behavior?


texts:
In this map I have added a line of code to make it zoom to the vector layer:

map.zoomToExtent(vector_Layer.getDataExtent());

This zooms in to the maximum zoom at the point in the Gulf of Guinea where the equator meets the Greenwich meridian. It is not supposed to zoom here, it is should show Scotland.

I've passed the same bit of code in the console in Chrome and that works fine.
Does anyone know why it is not working properly when the page loads and does anyone know what I can do to fix it?

Thanks.


texts: Correcting scale of heatmap data when plotting GeoTIFF file as map in R?

texts:
I'm trying to use R (and R studio) to examine a GeoTIFF file that I've gotten from NASA's Socioeconomic Data and Applications Center (SEDAC).

I'm a complete beginner in GIS, so as a first step, I was trying to get R to load in the tif and plot the same map that I've seen rendered as a PNG on their site:

I downloaded the data set from their page (though at a lower resolution, just to cut down on file size), and saved as pop-count-2020.tif in my R Studio workspace.

Then I ran this program:

install.packages(c("rgdal", "raster"))
library(raster)
library(rgdal)

pop_count_2020 <- raster(x = "pop-count-2020.tif")
plot(pop_count_2020)

And I'm seeing a plot like this:

I am getting some data rendered, which seems to line up with the darkest areas of the SEDAC original, so it feels like my scale is working at too large an increment, so most things are getting rounded down into the smallest bucket.

Is this actually what's happening or have I misunderstood?

If so, how can I correct this?

I'm not sure what to try next as I'm having a hard time trying to find the right terminology for what I'm trying to do.


texts:
my layers are in wgs84 and I need to zoomToExtend to  specific Bound's  but it seems I can not make it....
The  code seems ok but nothing happens.

map = new OpenLayers.Map({
            div: "map",
            fallThrough:false,
            units: "m",
            projection: googleMercator,
            displayProjection: wgs84
       });
.
.
.
.
        var minLat = 38.79;
        var minLng = 24.09;
        var maxLat = 42.08;
        var maxLng = 25.61;

        map.zoomToExtent(new OpenLayers.Bounds(minLng,minLat,maxLng,maxLat).transform(wgs84, googleMercator));


texts: Choropleth map turns one color when filtered, sql error

texts:
My data is generating a choropleth map perfectly but when I try to filter (both in the filter tab and by writing my own SQL) it turns all one color and tells me I have an SQL error.
This happens even when the SQL is generated by CardoDB's filter widget.
I have a big file - about 56,000 rows. Is this a file size issue or something else?


texts: How to overlay layers with different SRID using QGIS and PostGIS?

texts:
I have a map of the Netherlands imported from a SHP-file into my PostGIS environment with SRID 4289 (Rijksdriehoek Amersfoort). When I connect using QGIS to the PostGIS database and add the layer with municipality information, all goes well en looks great.
Also, I have a database table with a column with latitude and longitude information (I guess this is standard WGS84). This information I convert to a point in PostGIS SQL using: ST_SETSRID(ST_MakePoint(longitude, latitude), 4289).
But when I load a query on this table into the QGIS environment it plots the points at the wrong location (far below the municipality object and if I zoom in very much, I do see the shape of the Netherlands appearing). I have tried switching between 4289 and 4326 but no differences occur.

Is this a conversion issue between the different SRID's or is there something else?
What I can do to properly overlay the map with the points properly?


texts: Calling a custom tool from Python script

texts:
I have a script based gp tool published on web (via Flex App). The script goes to many different mxds to write a report.

Then I created another model to  do some geoprocessing in a different toolbox and mxd. Now from the first published tool I want to combine one, so that the whole process will execute once called from Flex app.


texts: Loop to increase value on Parameter of Model?

texts:
I'm working on the ArcGIS Model Builder to develop a model that makes a Raster Calculation using a Value Parameter and analyzes the maximum pixel value of the result, comparing to a Reference Value.

If the maximum pixel value is lower then the reference value, I need to increase "1" at the Variant Parameter and do the Raster Calculation again, and after to do it many times, until the maximum pixel value is higher or equal to the Reference Value. I'll continue to other calculations that I'm not working on yet.

That is a picture of the model, and I need to connect the "lower" box to the "Variant".


texts: Why is "define projection" in ArcGIS not working on my GeoTIFFs? Problem with the extent?

texts:
I have about 800 geoTIFF files that have an unspecified projection, but since I know the projection (WGS 84 UTM Zone 32N), I have used the “define projection” function to batch define the projections of these files.
Now I need to clip all of these files to a new extent using a shapefile with the same projection. However, when I checked the files I saw that they do not overlap. My shapefile shows up all the way at the top of the screen and the geoTIFFs are tiny dots in the middle of the screen. Actually, they are so tiny that I can’t see them in global view, but I know that they are there. The data has almost the same extent, so they should be the same size.

The spatial reference data and extent data is shown below for each file type. The difference I see is in the extents, but I don’t know how to fix this problem, or even if this is the source of my problem at all. Does anyone know what might be wrong and how to make my files overlap?

geoTIFF file
Spatial reference
XY Coordinate System
WGS_1984_UTM_Zone_32N
Linear Unit Meter (1.000000)
Angular Unit Degree (0.0174532925199433)
False_Easting 500000
False_Northing 0
Central Meridian 9
Scale_Factor 0.9996
Latitutde of Origin 0
Datum D_WGS_1984

Extent Top 6720.20271636
Left 370.036258324
Right 461.540519838
Bottom 6651.1691216

Shapefile
WGS_1984_UTM_Zone_32N
WKID: 32632
Authority: EPSG
Projection: Transverse_Mercator
False_Easting: 500000.0
False_Northing: 0.0
Central_Meridian: 9.0
Scale_Factor: 0.9996
Latitude_Of_Origin: 0.0
Linear Unit: Meter (1.0)

Geographic Coordinate System:
GCS_WGS_1984
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_WGS_1984
Spheroid: WGS_1984
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314245179
    Inverse Flattening: 298.257223563

Extent
Top 6733114.02831481
Left 364349.954284994
Right 501090.000365895
Bottom 6620131.02208354


texts:
I am having a vehicle route with a set of latitude and longitude values. User need travel on the specified route, If user deviates from the route (greater than 100 meters), we need to identify in which point the user moved away from the expected route. The user can again join on existing route after deviation. I need to calculate,

How far is the user travelling on the same expected route?
How far the user deviates from the expected route?

Is there any approach?


texts: In Geotools, why does reprojecting the origin of a GridGeometry2D results in off-the-grid coordinates?

texts:
I have a simple GridGeometry2D in Geotools 13.1 that is created like this:

GridEnvelope grid = new GridEnvelope2D(0, 0, 100, 100);
Envelope2D envelope = new Envelope2D(CRS.decode("EPSG:3857"), 10.0, 10.0, 1.0, 1.0);
GridGeometry2D geom = new GridGeometry2D(grid, envelope);

For a quick test, I'm trying to find the grid coordinates for point (10.0, 10.0):

DirectPosition out = geom.getCRSToGrid2D().transform(
                (DirectPosition) new DirectPosition2D(10.0, 10.0),
                (DirectPosition) new DirectPosition2D());

I expected the result for the x field to be 0.5 (since the default pixel anchoring mode is center), but instead this happened:

DirectPosition2D[-0.5000000000001137, 99.5]

I understand that the 99.5 is probably because the GridGeometry2D considers that y coordinates increase towards the north, but the GridEnvelope2D ycoordinates apparently increase towards the south (as it is common in image processing).

What I don't understand is why the x coordinate is negative if the minimal x for my grid is 0. Why is this happening?


texts: Library for minimum bounding box "by width" in JavaScript or Python

texts:
I'm looking for a library (or at least a code snippet) to calculate the minimum bounding box "by width" in JavaScript or Python. Below is an illustration from the ArcGIS documentation of that functionality. Ideally the function would give me also the angle, width and height of the bounding box. Any pointers?

Image courtesy of the ArcGIS documentation


texts: Is it possible for ArcGIS Portal to modify security settings of the source zip file?

texts:
Is it possible that during the upload of a zipped file geodatabase to ArcGIS Portal, that portal is able to change the security settings of the source zip file?
To walk through an example.
I have an existing Hosted Feature Service that is published and shared to different groups on ArcGIS Portal. I created this from a zipped File geodatabase.
The data needs updating.
So I use the 'Overwrite' tool inside of ArcGIS Portal admin site.
After I select the correct zipped filegeodatabase, the process runs for just under a minute and the hosted service is now back up and running with the latest data.
BUT
The source zip file - which was taken from a file share on a server, has had its security properties changed, so much so that even the Owner of the zip file is no longer myself (even though I created the file literally 5 mins before). I cannot delete the file or move it and its completely locked out from myself.
I'm at a complete loss to understand how ArcGIS Portal would modify the security settings as part of an upload.
I have experienced this 3 time now, but - I cannot successfully replicate the issue at will, so hence the question.
Is this even possible that ArcGIS Portal has the ability to change the source file security permissions as part of an Overwrite action?
Some additional details regarding Process Monitoring. I ran ProcessMonitor whilst uploading a zip file to ArcGIS Portal. This upload went successfully. the Process Operation that i was looking for, was a 'FileSystemControl' operation. However this process never occured on the zip file (which kinda makes sense, because the zip file permissions haven't changed).
So my only guess is that - its possible this operation only occurs when the actual 'error' is occurring.


texts: Register and display real time tracking of my smartphone (gps) on a map

texts:
I have 2 smartphones (android galaxy and iphone) and i would like to display the position, according to the gps coordinates, on a map.

Is there any app to download or any kind of service where i could register? Every time I start the app, the coordinates of the smartphone would be sent to some web service and few seconds later that actual position would be displayed on a map. The map should be accessed/viewed by a user from a desktop pc through a normal Internet Browser.

edit: any known open source projects where you could setup your own "tracking system" and expand/modify it?


texts: Copy all shapefiles and feature classes in a directory and its sub-directories to a geodatabase

texts:
The script I have put together should copy all shapefiles, dwg's, and feature classes (GDB) in an existing network directory to an existing geodatabase with a feature dataset regardless of the file. The script only needs to throw an exception or tell me that the file exists because it only has a copy function (arcpy.copy_management(indata,outdata)).

Here's the code:

import arcpy, os, shutil, time

arcpy.env.workspace = os.path.normpath('//svcfs01/datavolume2/Standards/Water Resources/GIS/_Water Resources GIS Database/New Folder')
targetGDB = os.path.normpath('//svcfs01/datavolume2/Standards/Water Resources/GIS/_Water Resources GIS Database/_Unprocessed_Raw_Data.gdb/WRSpatial/')
featureClasses = [arcpy.Describe(a).name for a in arcpy.ListFeatureClasses()]
feature_classes = []
targetPath = []
for dirpath, dirnames, filenames in arcpy.da.Walk(os.path.normpath('//svcfs01/datavolume2/Standards/Water Resources/GIS/_Water Resources GIS Database/New Folder'),datatype="Any",type="Any"):
    for filename in filenames:
        feature_classes.append(os.path.join(dirpath,filename))
        file, filext = os.path.splitext(filename)
        targetPath.append(os.path.join(targetGDB,file))
        if arcpy.Exists(targetPath)==False:
            try:
                arcpy.Copy_management(feature_classes, targetPath)
            except Exception as e:
                print "Unable to copy %s to %s" %(feature_classes, targetPath)
        else:
            print "%s already exists....skipping....." %(feature_classes)

The error I get is:

Unable to copy [u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\New Folder\\Drainage_Area.shp'] to [u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\_Unprocessed_Raw_Data.gdb\\WRSpatial\\Drainage_Area']
Unable to copy [u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\New Folder\\Drainage_Area.shp', u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\New Folder\\TC.shp'] to [u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\_Unprocessed_Raw_Data.gdb\\WRSpatial\\Drainage_Area', u'\\\\svcfs01\\datavolume2\\Standards\\Water Resources\\GIS\\_Water Resources GIS Database\\_Unprocessed_Raw_Data.gdb\\WRSpatial\\TC']

print e
Object: Error in executing tool


texts: Cropping attributes to only show those inside a layer

texts: I am currently learning to use ArcMap and I am trying to crop the lines (roads) in the below image to only those inside the state polygon. Is there a way to do this. I have tried via the attribute table, but this is as far as I could get.


texts: 'outdata = driver.Create' raster only outputs to byte even though set to Float64

texts:
I am trying to produce a raster in QGIS using the random forest classification predict.proba function. When I print these results in the python console it appears as decimals (between 0 and 1), however, when I open the raster in QGIS, there are only values of 0 and 1.
I have specified outdata = driver.Create(ProbSlides, ImageSize[1], ImageSize[0], gdal.GDT_Float64), but when I open the layer properties, the raster is in Byte format. Is there any other command that would get around this issue?
## part 13 predict probabilities
Probability_Full_Prediction = Estimator.predict_proba(FullFeatures)
print('Displaying Probability Prediction')
print(Probability_Full_Prediction)
Slides_Probability= Estimator.predict_proba(FullFeatures)[:,0]
print('Displaying Probability of Slides Prediction')
print(Slides_Probability)

##part 14 reshape and produce class raster
Full_Prob_Raster = (Probability_Full_Prediction.reshape(2,2968,3830))
print("Full Prob Raster")
print(Full_Prob_Raster)
Probab_Slides_Raster=(Slides_Probability.reshape(ImageSize))
print("Prob Slides Raster")
print(Probab_Slides_Raster)

#part 15 make into georeferenced image
fileName = os.path.join(MyFolder + MyImageRoot + BandList[0] + extension)
ImageFile = gdal.Open(fileName)
ProbFull = os.path.join(MyFolder + MyFullProbabilityRaster + '.tif')
ProbSlides = os.path.join(MyFolder + MyProbabilitySlidesRaster + '.tif')
driver = gdal.GetDriverByName("GTiff")
outdata = driver.Create(ProbFull, ImageSize[1], ImageSize[0], gdal.GDT_Float64)
outdata = driver.Create(ProbSlides, ImageSize[1], ImageSize[0], gdal.GDT_Float64)
outdata.SetGeoTransform(ImageFile.GetGeoTransform())
outdata.SetProjection(ImageFile.GetProjection())
outdata.GetRasterBand(1).WriteArray(Probab_Slides_Raster)
outdata.FlushCache() ##saves to disk!!
outdata = None

#part 16 open new Prob raster data in QGIS
print('Displaying Slides Probability')
fileInfoSlides = QFileInfo(ProbSlides)
baseNameSlides = fileInfoSlides.baseName()
rlayerSlides = QgsRasterLayer(ProbSlides, baseNameSlides)
if not rlayerSlides.isValid():
    print ('Slides Prob Layer failed to load!')
else:
    QgsProject.instance().addMapLayer(rlayerSlides)

print('Displaying Full Probability')
outdata.GetRasterBand(1).WriteArray(Full_Prob_Raster)
outdata.FlushCache() ##saves to disk!!
outdata = None
ImageFile = None
fileInfoFull = QFileInfo(ProbFull)
baseNameFull = fileInfoFull.baseName()
rlayerFull = QgsRasterLayer(ProbFull, baseNameFull)
if not rlayerFull.isValid():
    print ('Full Probability Layer failed to load!')
else:
        QgsProject.instance().addMapLayer(rlayerFull)


texts: Using Python to iterate through directory to create event layer and perform kernel Density for each file?

texts:
I am new to Python for ArcGIS, but have been using ArcGIS a while.
I have a large directory of tables with point data. For each file of each folder I want to create a kernel density layer.
I have managed to write a script that finds all the files (using os.walk) and creates an event layer however the resulting kernel density raster is blank, i.e. all cells have the same value.

When I use the GUI to create the kernel density, using the event layers that I created with the script the kernel density layer comes out beautifully.

I think it is the file formats that I'm having issues?? I tried different types but only .tif creates anything at all.
Here is my function to run the kernel density after the file has been found and the event layer created.

def RunKernel(feature):
inFeatures = feature
populationField = "NONE"
cellSize = 1.5E-3
searchRadius = 5000

exten = ".tif" # extention of output file
# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
# Execute KernelDensity
try:
    outKernelDensity = KernelDensity(inFeatures, populationField, cellSize,
                              searchRadius, "SQUARE_MAP_UNITS")

    featureName = feature.replace(".lyr", exten)

    print(featureName)
# Save the output
    outKernelDensity.save(featureName)

    print (feature, "done")
except Exception as err:
    print(err)


texts: Seeking dataset to create reproducible ArcPy example for ArcGIS Pro

texts:
I would like to create a reproducible example for ArcPy (for ArcGIS Pro) and cannot locate any example feature classes or rasters that come with installation.

I am looking for something like the TemplateData.gdb that came with ArcGIS for Desktop as answered at Seeking datasets to construct minimal reproducible examples for ArcMap.


texts: get lat/lon values of every start and end point of polyline

texts: I have a polyline shapefile with 166 different line features.  How would one populate the attribute table with x and y values of the start and end points of each line?


texts: Transform coordinates from custom projection to lat long PyProj

texts:
I am trying to convert a pair of coordinates specified in US survey feet on a custom transverse mercator projection to lat/long using PyProj. I must be doing something wrong because I cannot get the correct decimal degree values:

import pyproj

# Coordinates in ft and false easting/northing in ft
x1,y1 = [1383566.534, 1627687.047]
sourceProjection =  "+proj=tmerc +lat_0=42 +lon_0=-115.5 +k=0.9996 +x_0=1640416 +y_0=328083 +datum=NAD83 +units=us-ft +no_defs"

x2,y2 = pyproj.transform(pyproj.Proj(sourceProjection), pyproj.Proj(init='epsg:4326'), x1, y1)

print "Original:    {0}, {1}".format(x1, y1)
print "Transformed: {0}, {1}".format(x2, y2)

The output I get is:

Original:    1383566.534, 1627687.047
Transformed: -119.384930994, 53.6300677515

but the correct result is -116.503170  45.562201

What am I doing wrong please?


texts: Batch model stores incorrect model parameter values

texts:
I am trying to batch process about ten files using a model I created. There are five model parameters, which I can fill in for the first feature class without a problem. However, when I try to fill in the output file information for the subsequent feature classes, the model "remembers" the values for the first feature class and changes it. For example, for the output for the intersect of fc1, I specify fc1_intersect. Likewise with fc2, I specify fc2_intersect. As soon as I click away to fill in the next values, the output of fc2's intersect will be changed to the same as fc1's intersect.

This causes the model to fail. I have not specified any default output locations or anything. I used to have a similar problem in ArcMap, where the tools would "remember" the previous layers used, but closing and reopening it would solve it. This is not working now with ArcCatalog, and the time I was hoping to save by batch processing is now being wasted trying to figure out why it keeps changing my input.


texts:
Re-posting the question with updates.
I have a following piece of code that is allowing me to change the symbol properties.

layer = iface.activeLayer()
iter = layer.getFeatures()
layer.rendererV2().symbol().symbolLayer(0).setName('triangle')
layer.rendererV2().symbol().setColor(QColor.fromRgb(255,255,255))
layer.triggerRepaint

How do i get the following expression incorporated into the python script

u'size_dd_expression': u'CASE WHEN  "BCCH"  < 200 THEN 25\r\nWHEN "BCCH">     200 THEN 12 \r\nEND'

OR

u'color_dd_expression':u'CASE WHEN  "STATUS"  = \'ON\' THEN  color_rgba(0,255,0,255)\r\n  WHEN  "STATUS"  = \'OFF\' THEN  color_rgba(255,0,0,100)\r\nEND'

OR

 u'angle_dd_expression': u'"AZIMUTH" +180'


texts: Snapping doesn't work - MapInfo Pro 17 (64 bit)

texts: I've experienced a strange behaviour during drawing in MapInfo. Snap function doesn't work. Pressing the S key has no effect. I can't turn off snap. When I start drawing in snap mode, the program freezes. Any ideas what the reason is?


texts: What are Pros and Cons to ArcGIS Online Organizational Account?

texts:
Esri is making a big push towards ArcGIS Online as a means for map applications. They pushed it really hard at the recent Southeast User Conference and did a good job at highlighting its advantages.

Primarily: Configuration over Customization

I'll admit, this seems to be a significant ROI savings with as simple as some of the maps are to create with little or no programming needed.

What I'm wondering is, for any of you that might have experience with ArcGIS Online, what the advantages and disadvantages are to moving towards this product?

My organization is looking towards a 30 day trial of the product and want to be prepared if we decide to try and sell it to the decision makers.


texts: Debugging FMW memory usage near Group-based Transformers

texts:
I have a large FME FMW file which is using a lot of memory at runtime. I have been tasked with trying to find ways to reduce the memory footprint.

Early on in the workbench there is (what I assume is) a group-based transformer - LineOnAreaOverlayer, which looks for lines that intersect polygons. I suspect that the group-based nature of this transformer is pushing up memory usage, as the thousands of features which are fed into it need to be held in memory and processed all in one go (at least that's the theory). My aim is to halt processing when features are first output from this transformer to see how much memory is being used at this point. I'm using inspection points on the transformer's outputs.

However, it seems I either don't understand group-based transformers or inspection points. I break processing on the first feature that is output from the transformer, and the connector feature counts in the rest of the FMW all show features having already gone right through the transformation and being written. I thought that all input features were backed-up in front of the LineOnAreaOverlayer before anything else happened?

Can anyone help me clarify this please? Also, any other useful information on tracking memory usage in an FMW would be much appreciated - right now my only options are the log's peak memory usage numbers or trying to match a procmon log with the fme log, which will be a tedious task.


texts: Summing attribute values of overlapping polyline portions in ArcGIS Desktop?

texts:
I have polylines of 80+ bus routes operating in a city. The polylines have attribute values for the number of passengers, number of trips, etc. Many of the polylines overlap at certain areas of the road network. What I want to do is to save the overlapping portions into separate individual shapefiles, that have attribute values that correspond to the sums of those original values. I know how to do this the long way (find the overlaps, cut them and save the separate shapes, manually calculate the attribute values, enter and save the data) but is there a way to do it easier/faster without using scripts?

For example in the image above, the Green Route has 1,000 passengers and travels 20 trips a day while the Orange Route has 500 passengers and 10 trips a day. I'd like to cut out the overlapping area and save it as a separate shapefile with data on 1,500 passengers and 30 trips. This is just a simple example with two routes but I have 80+ in all that overlap in many locations with each other.


texts: ArcGIS Foreground progress bar sitting at 99%

texts:
Running the Interpolate Polygon to Multipatch tool and watched it work from 1% through 99% in 1% increments. Got to 99% and has been sitting there for a while.

Is it still working or did something break? Is this normal?


texts:
I have been tasked with getting OS Mastermap topography data into QGIS system to compliment our other map layers. The data comes in .gz format. Which I am told needs to be decompressed before I can do anything further? After decompressing it using WINRAR it turns it into a file which has no file type? I assume it is supposed to be .gml?

I have read I then need to convert it into shape file format so that QGIS can read it, so I use OS GML to OGR python plug in. Which seems to do nothing...or says "processing stage failed with code 1".

Can anyone give me any advice or a brief walkthrough on how to get OS Mastermap data into QGIS?


texts: QGIS: Rendering WMS/WMTS raster at proper resolution?

texts:
QGIS 2.18.24 – which I use for work and cannot update – seems to struggle with some WMS and WMTS layers. They seem to render at an incorrect resolution, although they appear aligned properly.
This is Mapbox (WMS):

And this is esri World Topo (WMTS):

As you can see (click to enlarge), both look really crappy. Seems like the resolution is way off. Other layers do not have such issues. Am I missing something obvious here?


texts:
I am georectifying scanned versions of Survey of India maps, which are in
Everest spheroid and india_nepal datum. The exact parameters are:

+proj=longlat +a=6377301.243 +b=6356100.2284 +towgs84=295,736,257,0,0,0,0
+no_defs

I encountered two issues, one minor and one more significant.

I define a custom CRS using these parameters. Give it a name
"Everest_India_Nepal_custom". Save the CRS. Use it when georectifying a
toposheet. Save the toposheet with this new CRS. (gtif file). But when I
reopen that gtif file in a new project, its properties show the CRS
parameters, but do not show the custom name I had given for ease of
recognition. When saving it said :User 100001: Everest...
But when opeining it says "User 100009: *Generated CRS....". So it does not
remember the name I gave. Although the name still shows up in CRS options.
bigger problem: if I try to validate this custom CRS by using the
"calculate" button at the bottom of the custom CRS window, i get rubbish.
But the CRS definition is correct, and it works--when my rectified topo is
placed on google earth, it sits perfectly (error <10m).


texts: Error when intersecting large spatial polygons data in R : rgeos_binpredfunc_prepared: maximum returned dense matrix size exceeded

texts:
I am trying to intersect two polygons that are SpatialPolygonsDataFrame. I am finally getting output after intersecting a small subset of the data frames. But when i run this on the entire dataset which is about 100,000 polygons in each of the data frame, i get an error message "rgeos_binpredfunc_prepared: maximum returned dense matrix size exceeded"
I tried splitting up the data frames using split and a group called district. I then ran intersect on a loop but that is so super slow.

sample dataset - https://drive.google.com/open?id=1Ha47ZqwG7bCuxExJk6zOs8KoqDiqD5A4

rm(list=ls())
##LOADING DATASETS REQUIRED FOR THIS
load("sampledata.Rdata")

#booths data - dataframe - upbooths_raw
#booths data - spatial file - up_booths_sample
#village data - up_villages_sample
#loading libraries that i might use
library(rgeos)
library(raster)

#new CRS which is good for Uttar Pradesh India
newcrs <- CRS("+proj=utm +zone=42N +datum=WGS84 +units=km")

#transforms booths and villages to new crs
up_booths_transform <- spTransform(up_booths, newcrs)
up_villages_transform<- spTransform(up_villages_sample,newcrs )

#create booths with 1km buffer
up_booths_transform_buffer1km <- gBuffer(up_booths_transform, byid = TRUE, width = 1, capStyle="ROUND")

up_booths_transform_buffer1km <- gBuffer(up_booths_transform_buffer1km, byid = TRUE, width = 0)
up_villages_transform <- gBuffer(up_villages_transform, byid = TRUE, width = 0)

#above is done as per suggestion on stackexchange because of the topology error that i get
#if i do not include this
#Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
#TopologyException: Input geom 1 is invalid: Self-intersection at or near point
#1463.9798370253977 3002.5185235944564 at 1463.9798370253977 3002.5185235944564

#Error in RGEOSBinPredFunc(spgeom1, spgeom2, byid, func) :
#rgeos_binpredfunc_prepared: maximum returned dense matrix size exceeded

#this runs perfectly on the small sample data i have attached but i get above message when
#i run this on the entire dataset
pi <- intersect(up_booths_transform_buffer1km, up_villages_transform)

My super slow approach to solving this error is painful...

#district is a higher geographic area that is meaningul and present in both dataframes, i really i would like to avoid this because district names differ in the data frames and i would have to clean the names to match before i use this strategy below.
boothlist<-split(up_booths_transform_buffer1km,up_booths_transform_buffer1km$district_name_17)
villagelist<-split(up_villages_transform,up_villages_transform$DISTRICT)

listofintersection<-list()
y<-1
#THIS IS SUPER SLOW
for (i in 1:length(boothlist))
{
  for (z in 1:length(villagelist))
  {
    if(names(boothlist)[i]==names(villagelist)[z])
    {
      temp1<-boothlist[[i]]
      temp2<-villagelist[[z]]
      #temp1 <- gBuffer(temp1, byid=TRUE, width=0)
      #temp2 <- gBuffer(temp2, byid=TRUE, width=0)
      listofintersection[[y]] <- intersect(temp1, temp2)
      names(listofintersection)[y]<-names(boothlist)[i]
      print(names(listofintersection)[y])
      y<-y+1
      print(y)

    }
  }
}


texts: Leaflet WMS layer not working until zoomed

texts:
I created a toggle button to add and remove wms layer from the map. In the default zoom level, it is working fine. But when I remove layer and zoom in or out, then do "Add Layer", the layer does not show up until map is zoomed in or out.
I tried layer.refresh(), but no improvement. Below is my sample code.

var map = L.mapbox.map('map', 'examples.map-4l7djmvo').setView([24, -89], 5);
var mylayer =      L.tileLayer.wms("http://nowcoast.noaa.gov/wms/com.esri.wms.Esrimap/analyses", {
format: 'image/png',
transparent: true,
layers: 'NSSL_RAS_QPE_72HR'
});

$('#Check').click(function() {
   if($(this).text() == "Add Layer")            {
         $(this).text("Remove Layer");
         map.addLayer(mylayer);
         mylayer.refresh();
       }
       else
           {
             $(this).text("Add Layer");
             map.removeLayer(mylayer);
              mylayer.refresh();
           }
});


texts: How to use GeoDjango PostGIS ForceRHR/ForcePolygonCW

texts:
I have a model with a field called 'geom' that is a MultiPolygonField. When I query for instances of the model, I'd like the data in this field to follow the right hand rule. Because I am using PostGIS on the back end, I have support for ForceRHR (I am using Django 2.06), but I can't figure out how to use it.

In views.py, I get a URL request with name value pairs, and I convert that dictionary into a dictionary suitable for use in a filter:

def field_datasets(request):
    client_name = 'MSFF'
    grower = Growers.objects.filter(name__exact=client_name).get()

    fields = serialize('geojson', Fields.objects.filter(grower__exact=grower.id).filter(
    Q(*[(k,v)for k,vs in request.GET.lists() for v in vs], _connector=Q.OR)
    ))

    return HttpResponse(fields, content_type='json'

I've tried this:

...
    test = Fields.objects.filter(grower__exact=grower.id).filter(
        Q(*[(k,v)for k,vs in request.GET.lists() for v in vs], _connector=Q.OR)
    )

    test.update(geom = ForceRHR('geom'))
    fields = serialize('geojson', test)
...

and I've tried this:

...
    test = Fields.objects.filter(grower__exact=grower.id).filter(
        Q(*[(k,v)for k,vs in request.GET.lists() for v in vs], _connector=Q.OR)
    )

    for item in test:
        rhr_geom = ForceRHR(item.geom)
        test[item].geom = rhr_geom

    fields = serialize('geojson', test)
...

The documentation says:
Accepts a single geographic field or expression and returns a modified version of the polygon/multipolygon in which all exterior rings are oriented clockwise and all interior rings are oriented counterclockwise. Non-polygonal geometries are returned unchanged.

So how else do you pass a geographic field?


texts: Wrapping vector tiles generated by geojson-vt

texts:
I've followed the example in geojson-vt, and am able to generate vector tiles from my geoJson source.  However, I can't read these tiles with mapbox-gl-js or Leaflet.MapboxVectorTile because I don't seem to be wrapping the tiles in a protocol buffer correctly.

Is there an example or package that automate this protobuf encoding of vector tiles?  I've looked at geobuf, but that does not appear to be able to handle the tiles from geojson-vt.


texts: Retrieve data from MongoDB in OpenLayers

texts:
I wonder if it's possible to retrieve data "directly" from MongoDB to display them in an OpenLayers Map? Maybe using some sort of API?

If you know more about this, or if you think it's a stupid question, let me know. Any kind of help will be appreciated.


texts: Is GDAL/OGR (specifically CoordinateTransformation) thread-safe?

texts:
I'm using the Java bindings if it matters. I have a heavily multi-threaded program and this is one of the bottlenecks that I have in a synchronization block. If I could remove the synchronization, it would improve speed a bit I believe, but I can't find any information whatsoever whether it's thread-safe or not, or how best to personally test for that.

The specific line would look something like: FROM_WGS84.TransformPoint(coordinates[0], coordinates[1]); where FROM_WGS84 is a CoordinateTransformation class created using CoordinateTransformation.CreateCoordinateTransformation(), but with arguments supplied obviously.


texts: Saving output files with unique names to avoid ERROR 000725 from ArcPy?

texts:
I'm a GIS Analyst who is VERY new to Python.  I have successfully run the script below.  However, when I want to run it again with a different input file it gives me the 000725 error...output file already exists.

How do I automatically create a unique output filename every time I run this script?

>>> #ADD NETCDF FILE POINTS TO THE MAP
... arcpy.MakeNetCDFFeatureLayer_md(in_netCDF_file="Z:/ENGINEERING/Stormwater/Stormwater Files/Rainfall/NWS NetCDF Files/nws_precip_1day_20180625_netcdf/nws_precip_1day_20180625_conus.nc", variable="observation", x_variable="x", y_variable="y", out_feature_layer="observation_Layer", row_dimension="x;y", z_variable="", m_variable="", dimension_values="", value_selection_method="BY_VALUE")
...
Runtime error  Traceback (most recent call last):   File "<string>", line 2, in <module>   File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\md.py", line 125, in MakeNetCDFFeatureLayer     raise e ExecuteError: ERROR 000725: Output Feature Layer: Dataset observation_Layer already exists.
>>>


texts:
recently I have installed the Geoserver with vector tile layer extension, getting the data in pbf (Mapbox vector tile format).

http://localhost:8080/geoserver/gwc/service/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&LAYER=VectorTileLayer:tbl__analysis&STYLE=&TILEMATRIX=EPSG:4326:{z}&TILEMATRIXSET=EPSG:4326&FORMAT=application/json;type=geojson&TILECOL={x}&TILEROW={y}

I want to convert this request URL to RESTfull like this.
http://localhost:8080/geoserver/gwc/service/wmts/1.0.0/blocks/{z}/{x}/{-y}.pbf

I have tried many methods but none of them working
My goal is to get the data in RESTful link for vector tiles .pbf layer
Update..
Logical reply in the question.
Geoserver currently doesn't support any restful service for URL mapping for wmts.


texts: How to ceate a point layer from a table in PostGIS?

texts:
I have a table in PostgreSQL(postgis). This table has latitude and longitude values in it.
Now I want to plot these point in qgis. The only thing I could get done was connecting to the database from QGIS.

Could anyone please let me know on how to go about it.Thank You.

Here is the screenshot below

Below is the connection error I get while connection to the database from QGIS


texts: Aggregating numerous polygons into one with PostGIS

texts:
I have this PostGIS query which aggregates many polygons according to a variable ("region"):

SELECT ST_AsGeoJSON((ST_Collect(polygons)),3) as geom FROM table GROUP BY region;

The ST_Collect returns a collection of polygons as it should but, since my goal in aggregating polygons is to reduce the number of vertices, how can I turn this collection of multipolygons into a single multi polygon?


texts: QGIS 2.8 Composer - extent overview not shown

texts:
I have a print composition with three maps:

a general map (1),
a medium-scale map (2)
and my detailed map (0).

I want a box showing the extend of map0, in map2.

For doing that I block the layer in map0 and then go to map2 and create a new Overview. Then I select map0 as Map Frame, but nothing happens; no extent frame is shown.

I've tried redrawing the map, also tried with different blending modes. I did it a couple of months ago, with QGIS2.6.

Could it be a bug of the new version?


texts: Run-time error -2147467259 (80004005) in VBA when clicking and saving polygon to shapefile?

texts:
I'm fairly new to this project and to VB and got some of old code to fix. The task is to click a polygon on the screen with the mouse and then save it to a shapefile.

When executing the code, it breaks with the error message:

  Run-time error -2147467259 (80004005) Automation error

The debugger turns to this line where the destination shapefile is opened:

Set pWorkspace = pWorkspaceFactory.OpenFromFile("I:\Projekte\LBDatenBank\Karten\Bildrahmen.shp", 0)

I: is a network drive.

I did some research and it seems that the WorkspaceFactory does not usually open shapefiles, yet this code used to work without changes a few months ago.

If I change the string to a path to a directory or something else the error stays the same. Even if I enter no path at all and just random characters.

I think the problem might be within the OS?


texts: Implementing Spatial Autocorrelation using QGIS or PostgreSQL or any other free application?

texts:
I new to GIS.

I have two layers in a PostgreSQL database.
The first layer contains points where events happened in Europe.
  The second contains a shapefile of Europe.

I have been able to map the points on the map of Europe using QGIS.
I would now like to implement Spatial Autocorrealation (Moran i).

I have seen a few demos of this being done using GeoDA but it appears that it only handles a single shapefile.

Can anyone point me in the right direction on how to implement this using QGIS or PostgreSQL or any other free application?


texts: Alternative quicker workflows to the standard Buffer tool in ArcGIS for Desktop 10.3

texts:
I have a point layer representing oil and gas exploration wells containing 500,000 features and have created a 50km dissolved buffer polygon to use in another workflow:

The completion time for the buffer was 16 hours, of which is unacceptable as the process would form part of a script.

Are there any quicker more efficient alternative workaround workflows that can be applied to reduce the processing time yet create the same end product?


texts: How to create a surface with a slope of 2%?

texts: I am working on a project where I have a polygon of a basin that I should create his base whose elevation begins at -1m and go down in every direction with a slope of 2%. Does anyone know how to do it automatically using ArcGIS?


texts: Why azimuth result is wrong for this pyqgis code?

texts:
In this code I want to obtain azimuth from a single point to each pixel center, but when run it, in raster created azimuth result for each pixel is wrong, what is correct code?

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

pntLayer = QgsVectorLayer("E:/Sample_Dataset_Lorestan/Point.shp", "Point", "ogr")
feat = pntLayer.getFeatures().next()
pntxy = feat.geometry().asPoint()

g_NoData_value = -9999

# Open tif file
ds = QgsRasterLayer("E:/Sample_Dataset_Lorestan/subset_a1.tif","subset_a1")
global g_pixelWidth, g_pixelHeight
g_pixelWidth = ds.rasterUnitsPerPixelX()
g_pixelHeight = ds.rasterUnitsPerPixelY()

# extent of the layer
ext = ds.extent()
global g_originX ,g_originY
g_originX ,g_originY = (ext.xMinimum(),ext.yMinimum())
cols = ds.width()
rows = ds.height()

# Get Center pixel x y
def pixel2coord(x, y):

  xp = (g_pixelWidth * x) + g_originX + (g_pixelWidth/2)
  yp = (g_pixelHeight * y) + g_originY + (g_pixelHeight /2)
  return(xp, yp)

# Open tif file
ds = gdal.Open("E:/Sample_Dataset_Lorestan/subset_a1.tif")

band1 = ds.GetRasterBand(1)

driver = ds.GetDriver()
outDs = driver.Create("C:/Temp_NoiseOnePoint/Azimuth_lyr.tif", cols, rows, 1, gdal.GDT_Float32)

outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.float32)

for i in range(0, rows):
  for j in range(0, cols):
    rspnt = pixel2coord(i, j)
    az = pntxy.azimuth(QgsPoint(rspnt[0], rspnt[1]))
    if az < 0.0:
        azim = 360 + az
    else:
        azim = az

    outData[i,j] = azim

# write the data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats

outBand.SetNoDataValue(g_NoData_value)

# georeference the image and set the projection
outDs.SetGeoTransform(ds.GetGeoTransform())
outDs.SetProjection(ds.GetProjection())
outBand.FlushCache()

#outBand = None
outDs = None

This picture shows point in green and azimuth raster result in grayscale, min value in white color is 0.5622 and max value in black color in 359.251,It seems that north side is incorrect!


texts: LSA-SAF data reprojection with GDAL and R

texts:
The Land Surface Analysis Satellite Applications Facility
(LSA SAF) distributes several products
derived from the EUMETSAT satellites data. I am working with R to read and process the data (the
file of this example is available from the
Static Auxiliary Data page
of the LSA-SAF service.)

library(rgdal)
library(raster)

r <- raster("HDF5_LSASAF_USGS_DEM_Euro.hdf")

Although the file is read correctly, the projection info is missing.

> r
class       : RasterLayer
dimensions  : 651, 1701, 1107351  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 1701, 0, 651  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : /home/oscar/temp/HDF5_LSASAF_USGS_DEM_Euro.hdf
names       : HDF5_LSASAF_USGS_DEM_Euro
values      : -32768, 32767  (min, max)

If I use gdalsrsinfo to parse the file and obtain the
projection, I get ERROR 1: ERROR - failed to load SRS definition
from ... which is more or less the same if I use GDALspatialRef
in R:

> GDALSpatialRef('HDF5_LSASAF_USGS_DEM_Euro.hdf')
[1] ""

As far as I know these files use the
GEOS projection.
However, if I set this PROJ.4 string and the file is transformed
to a latitude-longitude projection, the results are erroneous
(look at the extent of the projectExtent output):

projection(r) <- CRS("+proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs ")

> projectExtent(r, CRS(projLL))
class       : RasterLayer
dimensions  : 651, 1701, 1107351  (nrow, ncol, ncell)
resolution  : 8.983153e-06, 9.043695e-06  (x, y)
extent      : 0, 0.01528034, 0, 0.005887445  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0

It seems that the LSA-SAF data
are not yet geo-corrected
and the images are expressed in the "raw" Column-Row system, so
this could be the root of the problem. There are
Latitude and Longitude static files
which could be used for remapping the rasters. However, I would
like to use a more direct approach with a correct PROJ.4 string
and a projection transform with gdalwarp or rgdal::spTransform
in R.


texts: QGIS 3.10.2 - Install PostgreSQL and PostGIS

texts:
QGIS 3.10.2 in Ubuntu 18.04.

I'm having problems installing PostgreSQL and PostGIS following the instructions here: Installing PostGIS under Ubuntu.

This doesn't work, 'postgresql-9.1-postgis' isn't in the repository any more and has apparently been replaced with 'postgresql-10-postgis-2.4'.

Would it be safe to attempt an install as follows?

$ sudo apt install postgis

$ sudo apt install postgresql-10-postgis-2.4


texts: How to compare a classified raster dataset to another polygon dataset?

texts:
I am going to perform an unsupervised classification on my dataset and will get a map as a result. The objective is to rebuild an already existing map (of different soil types) with my new dataset. Unfortunately I don't have any underlying data information to the existing map - just the polygon features and their different classes.

My problem is: How can I validate which unsupervised classification result has a better correlation with the existing map?

The main issue will be in comparing the spatial patterns of the polygons and thereby derive which classification result represents the existing map better.
(I don't want to make just a visual comparison, but an analytic- mathematical one.)

My dataset from where I will perfom the unsupervised classification is raster based. The other (already existing) map consists of polygons (in shapefiles) with nothing more than the information of which type of class has been assigned.

I'd be happy if you could help me out with a proper methodology or process. (I am working either in R, QGIS and ARCGis).

After some research I came across a very promissing statistical approach: the Rand index.
In R there are various implementations (package mclust & flexclust) - I will work them through in the next days and hopefully can deliver more information on the topic in a while.


texts:
I have an SQL table containing ID, longitude and latitude (WGS84) of several points. I would like to add a new column containing the new latitude of every point projected 5Km north.

I see that PostGIS offers the tool ST_Project, which appears to do exactly what I would like to do. The issue is that I am not familiar with its syntax. I don't know how to apply ST_Project to all the points at the same time and how to add its output to a new column in the original table.


texts:
My goal is to be able to query cached tiles from a WMTS server using the EPSG4326 projection.
I would like to be able to use (1) the geographic coordinates of a point in degrees and a (2) zoom level (plus if required (3) the meta-information provided by the WMTS server), and obtain the (x,y) coordinates of the tile that contains my point.
For example, in EPSG3857, I can translate the coordinate (32.839833, -117.244648) at zoom level 13, corresponding to Mt Soledad in San Diego, to obtain the "tile coordinate" (3304,1428) using code provided on the internet (https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python). On online services I can retrieve the tile and verify that it contains the orginal point I was looking for:
https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/13/3304/1428
https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/WMTS/tile/1.0.0/World_Imagery/default/default028mm/13/3304/1428.jpg
I attempted the same thing on a EPSG4326 service. I published an image server on on a local ArcGIS server with cached tiles containing the same data, but this time using the EPSG 4326 projection. After, some poking around, I manually found the coordinates of my tile:
http://192.168.73.100/server/rest/services/imagery_sd_4326/ImageServer/tile/13/19750/15210
The scale and footprint of that tile are visibly different from the corresponding tile on arcgisonline. My questions are the following:

What is the meaning of the scale parameters that I see in some of the descriptions of the EPSG coordinate system (for instance: https://www.linz.govt.nz/data/linz-data-service/guides-and-documentation/wgs-84-web-mercator-tile-scale-set-definition)
In WMTS service pages, TileMatrixSet sections have attributes defined. Where can I find a complete list and description of these attributes. (for instance: ScaleDenominator, TopLeftCorner, MatrixWidth, MatrixHeight). I assume that these parameters have to do with tile layout.
In general, how can I use the information in the WMTS service page to calculate a tile coordinate containing a given GPS point.


texts: How to define an event for features selected in mapcanvas in qgis using python

texts: I have created a push button, when it is clicked it calls the function 'selectfeatures'. This function will call the Select feature by free hand tool of QGIS and allows the user to select features in the mapcanvas. What I need is my event to recognize the features selected in mapcanvas and call another function.  How should this be accomplished?


texts: Combine multiple dynamic Leaflet filters

texts:
I'm very new to Leaflet and made a map using a GeoJSON Layer. I'd like to filter points from the GeoJSON Layer by a checkbox input. If I consider only one Parameter (Par1) with levels A, B and C I can do it: I  different subsets of my layer by pre-defining all possible filters:

function filterA(feature, layer) {
            if (feature.properties.Par1 === "A") return true;
        };

The same I'd do for B and C, each resulting in a different Layer:

var dataA = L.geoJSON(oli_data.responseJSON, {
            pointToLayer: MakeMarker,
            filter: filterA,
            onEachFeature: MakePopUp
        }).addTo(map);

The same way I can create a subset for levels B and C using their filters. I let the user select the layer by:

            controler = L.control.layers(MapsTypes, {
            "A": dataA,
            "B": dataB
        }, {
            collapsed: true
        }).addTo(map);

Eventually I want to let the user filter the data by a combination of levels of Parameter1 and Parameter2. For that I would like to use Checkboxes, one with Parameter1 as title and all levels of Parameter1 as possible choices, and the same for Parameter2.

How can I achieve that without having to make extra layers for each possible combination of levels but rather by filtering by the actual user input?


texts: Making a popup (which is a DIV containing a panel consist of textbox and a button) to float within the map DIV only

texts:
The HTML for the popup is:

       <input type="text" value="" id="annotationText" />
       <button type="button" ng-click="setComment()">OK</button>

CSS:
.commentBox {
    position: absolute;
    z-index: 9999;
    left: 72px;
    top: 58px;
    display:none;
}

#commentBox:before {
    position: absolute;
    top: -7px;
    left: 1px;
    display: inline-block;
    border-right: 7px solid transparent;
    border-bottom: 7px solid #ccc;
    border-left: 7px solid transparent;
    border-bottom-color: rgba(0, 0, 0, 0.2);
    content: '';
}

#commentBox:after {
    position: absolute;
    top: -6px;
    left: 0px;
    display: inline-block;
    border-right: 6px solid transparent;
    border-bottom: 6px solid #ffffff;
    border-left: 6px solid transparent;
    content: '';
}

The commentBox:after and commentBox:before are the small arrow-tip and it's shadow on top-right corner of commentBox.
Issue: when I click on map on extreme right-side the popup-div shows outside of the main Map Div. Please refer below screenshot.

When clicked somewhere else:

How do I make this DIV to display within the limits of MAP DIV when clicked extreme right on MAP?


texts: Cutting area from polygon in Google Maps API?

texts:
I found this example here.

var map;
function initialize() {
  map = new google.maps.Map(
    document.getElementById("map"), {
      center: new google.maps.LatLng(37.4419, -122.1419),
      zoom: 13,
      mapTypeId: google.maps.MapTypeId.ROADMAP
    });
  var bounds = new google.maps.LatLngBounds();
  var triangleCoords1 = [
    new google.maps.LatLng(25.774252, -80.190262),
    new google.maps.LatLng(18.466465, -66.118292),
    new google.maps.LatLng(32.321384, -64.75737),
    new google.maps.LatLng(25.774252, -80.190262)
  ];
  var triangleCoords2 = [
    new google.maps.LatLng(29, -69),
    new google.maps.LatLng(25, -70),
    new google.maps.LatLng(23, -72),
    new google.maps.LatLng(29, -69)

  ];
  for (var i = 0; i < triangleCoords1.length; i++) {
    bounds.extend(triangleCoords1[i]);
  }
  for (var i = 0; i < triangleCoords2.length; i++) {
    bounds.extend(triangleCoords2[i]);
  }
  map.fitBounds(bounds);

  bermudaTriangle = new google.maps.Polygon({
    paths: [triangleCoords1, triangleCoords2],
    strokeColor: '#FF0000',
    strokeOpacity: 0.8,
    strokeWeight: 2,
    fillColor: '#FF0000',
    fillOpacity: 0.35
  });
  bermudaTriangle.setMap(map);
}
google.maps.event.addDomListener(window, "load", initialize);

html,
body,
#map {
  height: 100%;
  width: 100%;
  margin: 0px;
  padding: 0px
}

But if I change that to my needs, it will not work.

var map;
function initialize() {
  map = new google.maps.Map(
    document.getElementById("map"), {
      center: new google.maps.LatLng(37.4419, -122.1419),
      zoom: 13,
      mapTypeId: google.maps.MapTypeId.ROADMAP
    });
  var bounds = new google.maps.LatLngBounds();
  var triangleCoords1 = [
     new google.maps.LatLng(51.7339654, 6.4478929),
     new google.maps.LatLng(51.4992925, 6.4666535),
     new google.maps.LatLng(51.5837765, 6.1284897),
     new google.maps.LatLng(51.7339654, 6.4478929)
  ];
  var triangleCoords2 = [
     new google.maps.LatLng(51.6374564, 6.2919939),
     new google.maps.LatLng(51.6474904, 6.3974619),
     new google.maps.LatLng(51.5722125, 6.3693991),
     new google.maps.LatLng(51.6374564, 6.2919939)
 ];
 for (var i = 0; i < triangleCoords1.length; i++) {
    bounds.extend(triangleCoords1[i]);
  }
  for (var i = 0; i < triangleCoords2.length; i++) {
    bounds.extend(triangleCoords2[i]);
  }
  map.fitBounds(bounds);

  bermudaTriangle = new google.maps.Polygon({
    paths: [triangleCoords1, triangleCoords2],
    strokeColor: '#FF0000',
    strokeOpacity: 0.8,
    strokeWeight: 2,
    fillColor: '#FF0000',
    fillOpacity: 0.35
  });
  bermudaTriangle.setMap(map);
}
google.maps.event.addDomListener(window, "load", initialize);

html,
body,
#map {
  height: 100%;
  width: 100%;
  margin: 0px;
  padding: 0px
}

What am I doing wrong?


texts: Geoserver / Geowebcache diskquota.index.db H2 quota database file is 120GB, how to reorg / shrink or delete/rebuild?

texts:
I have Geoserver 2.8.2 running on a Windows server in Tomcat.  I have a drive specifically setup just for GeoWebCache.  All has been fine until today, I restarted the machine, and discover that drive is out of disk space.  It is 126GB in size, and was using about 6GB.  I discovered that the "diskquota.index.db" file is 120GB in size.  The diskquota.data.db is only 37MB. Clearly something has gone wrong.

I think I need to somehow reset / delete / recreate this database, but I'm not entirely sure what the proper procedures would be to do so.

Here is the directory listing for that folder:

05/22/2017  12:03 PM    <DIR>          .
05/22/2017  12:03 PM    <DIR>          ..
05/22/2017  12:03 PM        38,584,368 diskquota.data.db
05/22/2017  12:03 PM   123,511,767,088 diskquota.index.db
05/22/2017  12:03 PM           705,672 diskquota.trace.db
           3 File(s) 123,551,057,128 bytes
           2 Dir(s)         565,248 bytes free


texts:
I am having issues getting arcpy to run correctly on Linux Red Hat.  I solved some earlier problems which appeared to be related to permissions on the init_server file.  Now that si resolved, when I 'import arcpy' in Python, it just hangs, not even timing out after 20 minutes.

Does anyone have any experience of using arcpy on Linux and did you face any issues setting it up?

Thanks in advance for any help/suggestions.

EDIT

AN update.  I have rebooted the server and edited the LD_LIBRARY_PATH in init_server.sh and also changed the call to set_pyton_anv.sh (as it didn' appear to be calling it ) and it failed with Cannot open X display "(not specified)"  I then started an XWindows session on the box in Hummingbird and then, and only then, did import arcpy run correctly.

Not really sure what to say other than What?  Why?  How?

Does anybody know any reason why this should be?


texts: Mask clouds in an ImageCollection google earth engine

texts: Is It possible to create .median() Surface Refletance image (landSat 8 RS) without clouds before calculating the NDVI in goolge earth engine? I need to work with a year (eg. 2016-01-01, 2016-12-31), so I need to remove the clouds to create a single image free of clouds and then calculate the NDVI.


texts: GeoServer WMS Request - How to use alternate style?

texts:
I have a GeoServer map request for which I want to use an alternate style...

https://xyz.com/geoserver/wms?LAYERS=earth:WAYS&CQL_FILTER=[(USAGE IN ('X'))]&STYLES=STYLE2&SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image%2Fpng&SRS=EPSG:4326&TRANSPARENT=TRUE&_OLSALT=0.89383&BBOX=0,0,45,45&WIDTH=256&HEIGHT=256

I get back a  No such style: STYLE2

If I run the request at the GeoServer layer preview and change the style in the URL, I get back the proper result. It's just that it doesn't work when I try it through my network call. I can change the style in the network call URL to the built-in 'point' style and it works without problem. The only thing I can think of is that Geoserver is not looking in my workspace for the style, but rather the built-in, generic workspace where 'point' is stored. Perhaps there is a prefix I can place on the styles parameter, such as: ...&styles=earth:STYLE2... Tried that, but it didn't work. I think I also tried defining STYLE2 as a secondary style in the layer, but that didn't work either. The style is definitely created as I can see it in the GUI list of styles and it works in the layer preview.

Any ideas for the proper syntax or what I might be doing wrong?


texts: QGIS2Web cannot click url as popup disappears

texts:
I have successfully created an interactive map using QGIS2Web. I can hover over the polygons and see the url in the popup, but cannot click on it or access the url as the popup disappears when I maneuver off the icon and onto the popup.

I followed the tutorial to do this. There is an example map which you could see has this problem.

If you cannot click on the link but just having the url has no use.

Is there a way to fix this?


texts: Creating a Graph from a lines AND a points shapefile in networkx using read_shp

texts:
I am trying to use networkx for analysing networks which I represented in QGIS as a line shapefile and a point shapefile. If I import only the lines shapefile using networkx.read_shp('lines.shp') I can already do some of the analysis I want, but I cannot figure out how to do import also the point shapefile with the nodes. I need that because in some of my networks not all the edges go in a straight line from one node to the other and converting only the lines.shp networkx creates a node for each "bending" of the lines, and also because I would like to store in networkx both the coordinates and the attributes given to the points in the shapefile attribute table (exactly as when reading the lines.shp).

Both if I try networkx.read_shp('lines.shp','points.shp') or networkx.read_shp('points.shp') alone, I get the same error: "add_node() takes exactly 2 arguments (3 given)" referring to line 88 in https://networkx.github.io/documentation/stable/_modules/networkx/readwrite/nx_shp.html#read_shp ( net.add_node((g.GetPoint_2D(0)), attributes ).

Does anyone have suggestions on how to correctly read a point shapefile?

(In case that matters, I am using QGIS 2.18.17 with Python 2.7 and networkx 2.0)


texts: OSGeoLive Error: unable to locate package qgis-examples

texts:
I'm trying to upgrade my QGIS to the latest Testing (2.99) version via a complete reinstall.

I did this before by removing it first using

 sudo apt-get remove qgis* grass*

But now, I can get past this first step because I encounter an error saying something like:

E: Unable to locate package qgis-examples

, and there is a similar error about not being able to locate grassdata from apt-get. This prevents anything including a apt-get install -s. I've tried to revert back to a much earlier OSGeo Live VM (10.5) back in April/2017 when it was first released. Still, the same errors occur.

I don't understand how this could have happened. I've removed qgis and grass before without an issue (and had an older version of QGIS 2.99 successfully installed).

Has anyone encountered a similar issue when removing qgis, or know a fix for this error?


texts: Get catchment area of a lake polygon in R

texts:
I have several DEMs as geotiff files which I load into R using raster::raster. Additionally, I have lake outlines (located within the DEMs) in the form of shapefiles which I load into R using raster::shapefile.

What I would like to do is to identify the catchment of the lake within the DEM. Importantly, I would like to do this without making use of GIS tools (RGRASS, RSAGA, RQGIS), only R.

One approach I came up with involves the topmodel::subcatch function. This function allows me to get the catchment of any pixel in the raster. Therefore, I could use it to get the catchments of all the pixels which are within the lake polygon, or all the boundary pixels of the lake, alternatively.

However, since topmodel::subcatch only allows for a single pixel as input, I would have to perform it for all the lake boundary pixels and subsequently aggregate all the catchments. This seems inefficient. Is there a better way to do this (i.e. an R function similar to topmodel::subcatch which allows for multiple pixels, or a lake outline, as input)?

This approach would solve the problem, yet it uses R GIS tools, which I don't have access to.


texts:
I am trying to figure out how to add a search functionality to my openlayers map. I want the user to type in an address and the map to zoom into the location. I am looking at different geocoders and can't figure out a simple solution to this. Nominatim seems liek a possible solution. But I am trying to find something that I can develop on the front end. I have OSM data for my background and would prefer to use this. Any tip or ideas to point me in the right direction would be appreciated!

Thanks


texts: R: sf package points to multiple lines with st_cast

texts:
I want to create multiple lines out of given points as sf objects.

If I have a number of points as

library(sf)
pts <- st_multipoint(matrix(c(10, 10, 15, 20, 30, 30), nrow = 3, byrow = TRUE), dim = "XY")

and I am using st_cast to create lines of them

lines <- st_cast(pts, "MULTILINESTRING")

I will always get one sf object with multiple segments, but what I want to get is multiple lines (two in this example).


texts: Highlighted graphics resets on zoom in ArcGIS API for JavaScript?

texts:
I am highlighting graphic using graphic.setSymbol(pictureSymbol). It is working fine. But, when map is zoomed the symbol graphic is reset to default. If anyone has come over the issue, please help!
Below is the sample code which I'm using for zooming.

searchGraphic(e){
   var map = this.get('map');
   var layers = this.get('layers');
    var searchTerm = document.querySelector('.searchitem-node input').value;
   this.set('searchTerm', searchTerm);
   var selectedPictureSymbol = null;
   var pictureSymbol = null;
   var i = null;
   var j = null;
   var k = null;
   //if search button clicked
   if(searchTerm && e){
      e.preventDefault();
      searchTerm = searchTerm.toLowerCase();
   }
   for(i in layers){
      for(j in layers[i]){
         for(k in layers[i][j].graphics){
            if(layers[i][j].graphics[k].attributes.name){
               if(layers[i][j].graphics[k].attributes.name.toLowerCase() === searchTerm){
                  selectedPictureSymbol = this.highlightGraphics(layers[i][j].graphics[k], true);
                  this.set('graphic', layers[i][j].graphics[k]);
                  this.set('lineG', layers[i][j].graphics[k]);
               }else{
                  console.log('graphics ! found');
                  pictureSymbol = this.highlightGraphics(layers[i][j].graphics[k], false);
               }
            }
         }
      }
   }
   //if graphics not found then set graphics to blank for unselecting selected item
   if(searchTerm && e){
      if(!selectedPictureSymbol){
         console.log('set graphics blank');
         this.set('formClass', 'circuit-form');
         this.set('graphic', '');
         this.closeForm();
  }else{
      var graphic = this.get('graphic');
      if (graphic.geometry.type === 'point') {
         var maxZoom = map.getMaxZoom();
         map.centerAndZoom(graphic.geometry, maxZoom - 4);
      }else{
         var extent = GraphicsUtils.graphicsExtent([graphic]);
         console.log(extent);
         extent.xmax = extent.xmax + 50;
         extent.ymax = extent.ymax + 50;
         extent.xmin = extent.xmin - 50;
         extent.ymin = extent.ymin - 50;
         map.setExtent(extent);
      }
   }
}


texts: How do I make a picture marker symbol at a given latitude and longitude in ArcGIS for Javascript?

texts:
I have no idea how to use ArcGIS, so I spent about six hours trying to figure it out and got nothing.

The only thing I have working is showing the map.

arcgisUtils.createMap("39526eb31f1d42359ab82894999ae6a3","map").then(function(response){
      //update the app
      //dom.byId("title").innerHTML = response.itemInfo.item.title;
      //dom.byId("subtitle").innerHTML = response.itemInfo.item.snippet;

      var map = response.map;

      var graphicsLayer = new esri.layers.GraphicsLayer();
      map.addLayer(graphicsLayer);

      //test json
      json_object = {
             "displayFieldName" : "STATE_FIPS",
             "fieldAliases" : {
               "STATE_FIPS" : "STATE_FIPS"
             },
             "geometryType" : "esriGeometryPoint",
             "spatialReference" : {
               "wkid" : 102100
             },
             "fields" : [
               {
                 "name" : "STATE_FIPS",
                 "type" : "esriFieldTypeString",
                 "alias" : "STATE_FIPS",
                 "length" : 2
               }
             ],
             "features" : [
               {
                 "attributes" : {
            "STATE_FIPS" : "01"
                 },
                 "geometry" : {
            "x" : -113.5781,
            "y" : 37.0953
                 }
               }]};

      //end test json

      var point;
      var symbol = new PictureMarkerSymbol('img/Land/Lnd_icon_0.png', 200, 200);

      var graphic;

      for(var i=0; i<json_object.features.length; i++)
         {
          point = new esri.geometry.Point(json_object.features[i].geometry.x, json_object.features[i].geometry.y);
          graphic = new esri.Graphic(point, symbol);
          graphicsLayer.add(graphic);

         }

      //add the scalebar
      var scalebar = new Scalebar({
        map: map,
        scalebarUnit: "english"
      });

      //add the legend. Note that we use the utility method getLegendLayers to get
      //the layers to display in the legend from the createMap response.
      var legendLayers = arcgisUtils.getLegendLayers(response);
      var legendDijit = new Legend({
        map: map,
        layerInfos: legendLayers
      },"legend");
      legendDijit.startup();
          dojo.connect(map, 'onLoad', function(theMap) {
          //resize the map when the browser resizes
          dojo.connect(dijit.byId('map'), 'resize', map,map.resize);
        });

      dojo.addOnLoad(init);

    });


texts: Recovering output value in QGIS graphical modeler?

texts:
I use QGIS 3.0.3

In the graphical modeler, I use the algorithm "Basic statistics for fields".

Is there a way to recover the export values in the produced html file and using them in the modeler?


texts: Looking for a .NET Dijkstra component using Shapefile

texts:
I'm totally new to stackexchange and i'm in this beautiful informative place for one problem.
I'm a .net programmer and i need to write a program that reads a shapefile and performs routing operation. I'm looking for an open source solution.
Is there any idea?


texts: Looking for geocoding data for New York State

texts:
I have tried the following local data sets:

TIGER 2010
TIGER 2006
www.datasciencetoolkit.org
http://www.gisgraphy.com

They are all very inaccurate for data in New York State. I need the data to work at the street/house level similar to google maps.

I have searched all over the place and it does not seem this type of data exists.

Any recommendations?


texts: Need to give same color for Points and its label

texts:
I have 10000 points, using Graduate Color I have given different color for each points in 10000. I want to give a same color for each label according to each points colors using ArcGIS. is it possible, How to do it.

the above image is model for my work


texts: Storing earth observation derived attributes for polygons

texts:
I want to generate a number of vegetation related indices (NDVI, NDTI, etc., min/max/avg/std (pixel values)) from earth observation data for a large number of polygons during the entire growing season, roughly between March and November each year. The figures are about as follows: I have about one million polygons. Every day I will get new EO data (Sentinel-1/Sentinel-2) for about 20% of them. For each of these polygons I generate 10-20 indices based on the EO data. This gives me appr. 2-4 million records, every day. That makes appr. 500 - 1000 million during just one growing season (I´ll need to store at least 5 seasons).

The infrastructure within which I have to operate is predetermined and will have to be something based on either Oracle(Locator) or PostGIS. Personally I´d prefer PostGIS since OpenSource allows for much more flexibility.

My initial idea is to create a PostGIS database, which is partitioned based on year value. I thought about creating one attribute table where I create a new row for each date and each interpreted property (polygon geometries+id are stored in a separate table). It would look something like this:

Since I have to do different interpretations depending on the geographical zone where the polygon lies, I also thought about creating a separate table for each zone. This will however make querying more difficult.

My questions are hence:

Does Oracle (Locator) or PostGIS as base for all this make sense
at all or do I need to start asking for an account at an ESA
DIAS/Google Earth Engine/AWS in order to be able to use cloud
solutions?
If this indeed makes sense, wht is your opinion on my planned
table structure?


texts: Need help assigning values to field using OGR in python

texts:
I've got a similar question to this one How to programmatically populate a shapefile field with polygon areas in OGR?, but the solution doesn't seem to be working for me.  I am trying to import a shapefile from PostGIS into python using OGR, add a field, then populate the field with data from a csv based on an ID and the date.  When I print the fields for the feature during the loop, the field gets the correct data but it isn't writing to the layer o3_proj (last two print statements).  Am I using the SetFeature function incorrectly?  I'm new to PostGIS and ogr, so please speak slowly!

import os, ogr, gdal, csv
driver = ogr.GetDriverByName('PostgreSQL')
datasource = driver.Open("PG: host='localhost' port='54321' dbname='chicago' user='postgres'", update = 1)

o3_proj = datasource.GetLayerByName('o3_monitors(geom3528)')

date = []
maxozone = []
aqi = []
aqs_site_id = []

f = csv.reader(open('/ozone/ozone_data.csv','rU'))
x = list(f)

for iii in range(1,len(x)):
    date.append(x[iii][0])
    aqs_site_id.append(x[iii][1])
    maxozone.append(float(x[iii][2]))
    aqi.append(int(x[iii][3]))

daily_field = ogr.FieldDefn('temp_o3',ogr.OFTReal)
o3_proj.CreateField(daily_field)

for feat in o3_proj:
    for jjj in range(0,len(date)):
        if aqs_site_id[jjj] == feat.GetField(0) and date[jjj] == '03/06/2005':
            o3_proj.SetFeature(feat)
            feat.SetField2(6,maxozone[jjj])
            o3_proj.SetFeature(feat)
            print str(feat.GetField(0))+','+str(feat.GetField(1))+','+str(feat.GetField(2))+','+str(feat.GetField(3))+','+str(feat.GetField(4))+','+str(feat.GetField(5))+','+str(feat.GetField(6))
            print o3_proj.GetFeature(feat.GetFID()).GetField(6)

Result of the print statements:

17-031-0064,41.79079,-87.60165,2001,2012,None,0.036
None
17-031-0076,41.7514,-87.71349,2004,2012,None,0.038
None
17-031-4007,42.06028,-87.86323,2002,2012,None,0.035
None
17-031-4201,42.14,-87.79922,2001,2012,None,0.031
None


texts: Data structure for analysis of relationship between rainfall and vegetation indices?

texts:
I have several hundred raster files (NDVI) of the same area dating back to the 1980’s, as well as daily rainfall observations from a nearby weather-station over the same period.  The raster files are in the same folder and named by date, YYYY-MM-DD (I can re-name these if it will simplify things).  I want to look at the relationship between rainfall and change in vegetation index over time at various points in that landscape using R, but I don’t know how to process the data efficiently.

Ideally the raster files could all be called from a single object, like a rasterstack, in which was also recorded the year and Julian day number (ie 1 to 365/6) for each image.  I could then extract information from rasters from specific time periods or from the whole series, and / or relate change in vegetation index to specific rainfall events (which data I already have in Julian day and year in a data frame).

I can bring the raster files into R, named as they are in the original folder and stored in a list object, from which I can then make a rasterstack, but this doesn’t attribute them a date.  Can I append a RAT to each raster when it is processed into the list, and automatically add year and Julian day attributes?

Can anyone suggest appropriate structures / functions to perform these tasks?


texts: Adding constraints for least cost path analysis

texts:
I want to calculate least cost path between two points A and B using a DEM from which I got my slope, cost distance and back link raster. However, the resulted least cost path crosses through two lakes and I don't want that. I want them to be constraints, so the line should find another way to the destination point.

I created two shapefiles for the lakes, but how do I make the least cost tool take them into consideration?


texts: severe thunderstorms warnings dataset for europe

texts: I'm looking at the possibility of using georef data on historic thunderstorms alerts. Does anyone know of a good (free) source for this information?


texts: Scale factor for EVI in Google Earth Engine

texts:
I'm new to GEE. I'm  working on EVI data analysis. So this is my code for MODIS evi. I want to know how to multiply all the images there inside ImageCollection by scale factor 0.0001 to get real EVI values.
Because my values are -2000, 10000

How do I get this work done?

var clipper = function(image){
  return image.clip(sureste);
};

Map.addLayer(sureste);
print(sureste);

var modis = ee.ImageCollection('MODIS/006/MOD13Q1')
  .map(clipper);
var months = ee.List.sequence(1,12);

// Group by month, and then reduce within groups by mean();
// the result is an ImageCollection with one image for each
// month.

var byMonth = ee.ImageCollection.fromImages(
      months.map(function(m) {
        return modis.filter(ee.Filter.calendarRange(m, m, 'month'))
                    .filterBounds(sureste)
                    .select('EVI').mean().multiply()
                    .set('month', m);
}));

print("Monthly images",byMonth);


texts: Orientating Contour Labels Uphill in QGIS?

texts:
I'm quite new to QGIS and I'm trying to style topographic contours. I have managed to style the lines the way I'd like but I'm having trouble with the line labels - I can't style them so the labels are orientated uphill.

This question is similar to How to place elevation numbers on contours with uphill orientation and position? except I want to be able to export the style as an SLD as well. I hope the image below explains what I'm trying to achieve:


texts: Transforming ArcGIS date field to integer?

texts:
I have a featureclass with a 'date' class field named 'ACQ_DATE' formatted as:

1-8-2013
24-7-2013

etc. As you can see the number of digits is not harmonized. I want to display the dates on a graduated colour scale in the symbology (i.e. earliest date: red, latest date: green), but I understand that graduated colours are not applicable to a date field. Therefore, I want to transform the date to an integer which is supported by graduated colours. The integer should be in the format:

20130801
20130724

I've tried the following python code block, but it doesn't work.

date=str(!ACQ_DATE!)
year= date.split("-")[2]
month= date.split("-")[1]
day= date.split("-")[0]
dateInt= year * 10000 + month * 100 + day

This code gives a syntax error, but I'm not sure what's wrong with it.

How do I convert a date field to an integer field?

A good solution was provided by a colleague. Write in the calculator:

datetime.datetime.strptime(str(!ACQ_DATE!), "%d-%m-%Y").strftime("%Y%m%d")


texts: Converting geodatabase public lands to a shapefile

texts:
Can a whole geodatabase be converted to a single shapefile? This is my code, not sure what I am doing wrong...

# Import system modules
import arcpy

# Set environment settings
arcpy.env.workspace = "D:/UW/Winter_2016/501/pythonProject/publicLandsData/pLDWorkspace"

nfc = "publicLandsDataToGDB"
fp = "D:/UW/Winter_2016/501/pythonProject/publicLandsData/pLDWorkspace/publicLandsDataConversion.gdb"

arcpy.CreateFeatureclass_management(fp, nfc, "Polygon","","","", "WGS 1984")
arcpy.FeatureClassToShapefile_conversion("D:/UW/Winter_2016/501/pythonProject/publicLandsData/waParksGDB/WA_PublicLandsInventory_2014_r1.3_UWSEFS_PublicInfo.gdb",
                                         "D:/UW/Winter_2016/501/pythonProject/publicLandsData/newShapefiles")

the error is saying Execute Error which I believe is an arcpy function, that the input features: Dataset: "variouspathname.gdb" does not exist or is not supported..


texts: Calculating point proximity in ArcMap 10.1

texts:
I have a geocoded point file which (theoretically) contains one point for every zip code in the U.S. I want to create a smaller list of zip codes from this file. The list should include only those zip codes that are 100 miles from one another. In other words, if I were to look at this on a map, I would only see one point for every 100 miles, and any points in between would not be included. I'm having trouble getting this part of my project started in ArcMap 10.1. How would I go about creating this output?

I thought that I could use the Near or Point Distance tools, but I always run into trouble when trying to use the same file as both my inputs. Also, I thought that I could drape a 100x100 mile grid on top of the point file, and then grab the closest one to the center of each panel, but I'm not sure how to set that up either.

Any ideas?


texts: Extract certain area of contours/satellite images

texts:
I have 2 satellite images (jpeg2000 format) and 2 xyz-files from the same area. I am trying to extract certain region from the area including contour lines and the satellite image in order to import these into Google Sketchup for terrain modeling.

I have been able to create contour lines from xyz-files. Anyway, I would like to first merge contours and then satellite pictures and finally extract only certain part of the data (both contours and satellite pictures). Is there some kind of instruction available that would help in a problem like this?


texts: Adding dynamic picture elements to mapbook

texts: I have a feature class full of traffic signal data including up to 3 pictures of each traffic signal.  I want to create a series of data driven pages that include 3 picture/data frames that dynamically update the pictures based on HTML link in the attribute table.


texts:
I need to do routing for a vehicle. A map with detailed turn-by-turn driving directions should be generated. For example, it should be like car race game in which the map with vehicle tracking is provided and driving instructions are shown by arrow keys. Is it possible with QGIS?

EDIT:

QGIS gps tracking plugin is useful to track vehicle. Please, can anyone clearly explain how I can provide driving direction to the driver (for example: turn left, turn right for given start point and end point).


texts:
I’m experimenting with some rule based styling which is working really well to a point but I’m hoping it can do more. I work in school transport I’ve plotted all pupils in a certain school on the map. With categorised styling I’ve then made pupils on different transport routes different colours with the route number as text in the middle of the symbol.  This categorised styling has also created ‘sub layers’ for each route after clicking classify which is great.

What I’m hoping to add is a bit more complex I think. I can include information in my data to show if a pupil is a wheelchair user. If they are I want to show them on the map using a wheelchair symbol. So if they are a wheelchair user they will show up as a wheelchair symbol. If they are not then it will revert to my previously mentioned categorised styling.

I’d quite like to add a few more rules like this potentially as well so there may be a couple of rules which will determine the symbol and if none are met then categorized styling is reverted to.


texts: Writing a stand-alone PyQGIS script in Python 3.6 instead of 3.7?

texts:
I am currently using Python version 3.6. Previously used QGIS 2.99 as it was throwing errors, shifted to QGIS version 3.4. With the help of this "How can I discover the Python version in QGIS?" , I found that QGIS 3.4 has Python version 3.7.0.
By this I understood QGIS 3 will use Python 3.

Can Python 3.6 be used in QGIS version 3.4, will it support a standalone script?


texts:
On my Mac, running OS X 10.10.5 and I installed Grass with homebrew:

brew tap osgeo/osgeo4mac
brew install grass-70

Unfortunately, I cannot choose the pg driver in the db.login menu.

Can anyone tell me how to install it or how to brew Grass with the pg driver activated?


texts: Raster to Point in ArcMap: export raster data to a point shapefile with geometry of a given point shapefile

texts:
I want to export a calculated raster to a point shapefile with predefined geometry (i.e. positioning and extent of the points). Points not covered by the raster should return empty values (e.g. NoData). It works for the positioning of the points, however it does not work for the extension of the points.

I do have:

point shapefile MyTemplatePoints.shp as a template
raster file mytempr calculated from template point shapefile above
raster file mycalcr calculated from various input data in an arcpy-script
point shapefile MyCalculatedPoints.shp exported from mycalcr with Raster to Point (Conversion)

The arcpy-script to calculate mycalcr contains the following command line:

arcpy.env.snapRaster = mytempr

By this, the centroids of the cells of mycalcr snap to those of mytempr. The two rasters have different extensions which is not a problem. Export of mycalcr to MyCalculatedPoints.shp is performed with Raster to Point (Conversion) with input parameters as shown to the left in the image. The points of MyCalculatedPoints.shp fit to the positioning of the points of MyTemplatePoints.shp as visible to the right in the image. However, the extension of MyCalculatedPoints.shp does not match the extension of MyTemplatePoints.shp as I would expect by Environment Settings > Processing Extent.

It seems to me that this processing is unsuitable to generate a point shapefile with larger extension than the raster itself (and to fill these points with empty values, e.g. NoData).

What am I missing?


texts: Deleting selected features from Vector/OGR in GDAL/Python?

texts:
I run the following code to first select specific attributes using a query and then I need to delete them

from osgeo import gdal, ogr
import os

##shapefile = "D:/New.shp"
##driver = ogr.GetDriverByName("ESRI Shapefile")
##dataSource = driver.Open(shapefile, 0)

shapefile = ogr.Open( "D:/New.shp")

layer = shapefile.GetLayer()
layer.SetAttributeFilter("Area < 5000")

for i in layer:
    print i.GetField("Area")
    layer.DeleteFeature(i)

I get the following error :

layer.DeleteFeature(i)   File
"C:\Python27\ArcGIS10.2\lib\site-packages\osgeo\ogr.py", line 1499, in
DeleteFeature
return _ogr.Layer_DeleteFeature(self, *args) TypeError: not an integer

I changed the last line of my code to ogr.Layer_DeleteFeature(i) to which i get the error

ogr.layer_DeleteFeature(i)
AttributeError: 'module' object has no attribute 'layer_DeleteFeature'