texts: Georeferencing imported DWG layer using QGIS

texts:
I would like to georeference a DWG layer using QGIS 3.
I have imported the DWG layer saved by Microsoft Office Visio Professional.
After import I have a nice lookup like below.

The problem is, that I would like to plot this layer with an existing GeoTiff raster layer in QGIS.

I found some hints:

Importing a georeferenced DWG/ DXF into QGIS
https://www.sigterritoires.fr/index.php/en/how-to-integrate-dwg-data-in-qgis/
Georeferencing vector layer with control points using QGIS?
http://osgeo-org.1560.x6.nabble.com/how-to-georeferenced-an-autocad-dxf-file-td4530764.html
https://www.youtube.com/watch?v=KkPdDR87kXo

and most of them indicate firstly the conversion to DXF or even SHP and next georeferencing.
I am not happy with that solution, because this DWG layer is good-looking and keeps all shapes and colourations the same as they appear in Visio.
Is there any way to georeference them similarly to the raster georeference in QGIS, where I can pick up some points from the maps and assign them to my DWG layer?
I am interested in the vector georeferencing in this case. I checked the Vector bender plugin, but it doesn't work this way.


texts:
I converted shapefile to GeoJSON, the shapefile is 10mb and the result of GeoJSON was around 30. In order to make this conversion, I used the following code:

import geopandas as gpd

SHP_FILE = './00ent.shp'
JSON_FILE = './mexico.geojson'

geodf = gpd.read_file(SHP_FILE)
geodf = geodf.to_crs(epsg='4326')
geodf.to_file(JSON_FILE, driver="GeoJSON")

How can I shrink this file with GeoPandas or another Python library? I used an online tool, and I can reduce it to 2MB, but I need to do this programatically.

Any solution?


texts: Difference between 'point' and 'WKSpoint' in ArcGIS

texts: I have to add 3-D points to make a polygon in ArcGIS using VB.NET. I have defined an IpointCollection and want to add points to this. I have options to add 'point' or 'pointcollection' or 'WKSpointz'. Which one should I use? Please explain this as I am new to this field.


texts: Finding regions or zones in raster DEM of similar slope or aspect values in ArcGIS for Desktop?

texts:
I have a detailed 5m DEM for my study area. I want to find regions or zones where the slope or aspect (or preferably both) varies little. In other words, a nice evenly sloped hillside with a consistent aspect should show up high value, whereas complex, rugged terrain should show up as low. I want to identify these areas of smooth terrain, and preferably rank them by some kind of area. I guess that means I have to specify some kind of threshold or tolerance of values which I consider to be similar - I don't expect the slope or aspect to be exactly the same at any point, particularly not with floating point values.

What kind of filter do I need to pass over my data to identify these rough areas? I have spatial analyst and 3D analyst.


texts: How do I edit a Direct Connection Workspace using arcpy.da.Editor?

texts:
How to start editor to edit version data for ArcGIS 10.1 direct connection. I tried as given below code but it is throwing an error "Operation only allowed by the owner of the version [SDE.DEFAULT]". We are able to edit features manually.

workspace = "Database Connections/Connection.sde"
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()


texts: Going from degrees to meters in projection conversion using GDAL?

texts:
I have raster data in a Lambert-Conformal conic projection created from a NetCDF file by running sub-processes in a Python app that runs the following (Linux) command line argument:
gdal_translate -ot "Float32" -of "GTiff" -a_srs ${MY_WKT} nc_file.nc tif_file.tif

where:
]$ echo ${M_WKT}
]$ PROJCS["unnamed",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",13.34249973297119],
PARAMETER["standard_parallel_2",47.31750106811523],
PARAMETER["latitude_of_origin",30.33000183105469],
PARAMETER["central_meridian",-86.16000366210938],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]

and a snippet of the "gdalinfo" dump is:
Attributes:
grid_mapping_name:              lambert_conformal_conic
earth_radius:                   6370000.0
standard_parallel:              (13.3425, 47.3175)
longitude_of_central_meridian:  -86.16
latitude_of_projection_origin:  30.330002

pressure_interp#coordinates=XLAT Time level XTIME XLONG
  pressure_interp#FieldType=104
  pressure_interp#grid_mapping=LambertConformal
  pressure_interp#missing_value=9.969209968386869e+36
  pressure_interp#units=hPa
  pressure_interp#vert_units=m
  pressure_interp#_FillValue=9.96921e+36
  south_north#axis=Y
  south_north#units=m
  south_north#_FillValue=nan
  west_east#axis=X
  west_east#units=m
  west_east#_FillValue=nan
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -937500.000,  937500.000) ( 97d20'18.12"W, 38d41'55.14"N)
Lower Left  ( -937500.000, -937500.000) ( 95d28'43.54"W, 21d 9'29.33"N)
Upper Right (  937500.000,  937500.000) ( 74d58'53.91"W, 38d41'55.14"N)
Lower Right (  937500.000, -937500.000) ( 76d50'28.49"W, 21d 9'29.33"N)
Center      (   0.0000000,   0.0000000) ( 86d 9'36.01"W, 30d19'48.01"N)
Band 1 Block=150x13 Type=Float32, ColorInterp=Gray
  NoData Value=9.96920996838686905e+36
  Unit Type: hPa

I want to both be able to apply a color scheme for the pixels based on the floating-point data and convert the projection to web mercator (EPSG:3857), so the goal is to use "gdalwarp" and then "gdaldem color-relief" with my color file.  I seem to be missing a step as I try the following:
gdalwarp -ot "UInt16" -t_srs "EPSG:3857" tif_file.tif new_tif.tif

and receive the following error:
ERROR 1: The transformation is already "north up" or a transformation between pixel/line
and georeferenced coordinates cannot be computed for /path/to/tif_file.tif

I was suggested to add "-to SRC_METHOD=NO_GEOTRANSFORM" but this is not what I am trying to accomplish (and doing so just to see yielded the same error).  When I use "gdal_translate" instead the same error comes up because it creates a virtual warped file.  Adding "-a_ullr" seemed to work and the color shading code is fine, but it wants meters and not degrees. Thus putting the degree extent of the data as the corners as arguments for "a_ullr" puts the image basically at the Equator/PM intersection since -97 longitude is considered 97 meters west of the PM.  Is there an intermediate step I need where this projection can be converted into meters and then using "gdal_translate" or "gdalwarp"?


texts: Fractal Net Evolution Approach (FNEA) Citation

texts: I've been looking for some original source to cite for FNEA (Fractal Net Evolution Approach), instead of just citing eCognition.  Does anyone know of a paper on the origin of FNEA?  My searches just turn up papers using FNEA.  Perhaps the FNEA name was coined later?


texts: How to move all elements of PageLayout to left (right, right, up, down..) 5 cm (C#)?

texts: Are there anyone know how to move all element (text, line,) of a PageLayout (in GraphicsContainer) to left (or right, up...) 5 cm using C# ?


texts: How to read attributes in pyQGIS and QGIS 1.8 without setSelectedFeature

texts:
I add a vector layer to a map with

QgsMapLayerRegistry.instance().addMapLayer(self.vlayer)

later the attributes are read with

def select_all(layer):
layer.select([])
layer.setSelectedFeatures([obj.id() for obj in layer])

and

layer = self.vlayer
select_all(layer)
dict={}

for elem in layer.selectedFeatures():
    attrs = elem.attributeMap()
    for (k,v) in attrs.iteritems():
        dict.update({k:v})

The problem is that it selects all the feature in the map, too (which highlights them yellow). I tried to use

layer = self.vlayer
layer.select()
dict={}

for elem in layer:
   attrs = elem.attributeMap()
    for (k,v) in attrs.iteritems():
        dict.update({k:v})

but i get only an empty dictionary. Is it possible to get the values without selecting everything?


texts: Consuming Network Analyst Service not Returning Features in ArcGIS Server?

texts:
I am attempting to consume a rest service that I created here. Its just a simple network of streets in the Utah County, Utah area. Basically, it just shows the user the closest EMS service to the location that they choose. The model works in ArcMap 10.6, I simply click a point on the map and the service adds a model to my map with the EMS and closest route layer (Pictured below).

When I try to consume it as a service it doesn't quite work as I would expect though. I am just following the workflow that I found, but the object that is returned from the "Closest_Facility__2_" parameter does not seem to have any feature data in it. I have included a codepen here that shows all my code, unfortunately it cannot be run on codepen.io because of https/http errors, but if exported as a .zip, it can be run locally and it will work.


texts: Landsat 8 - NDVI Global Servises - Clouds

texts: I'm looking for Landsat 8 NDVI Global Services, like ArcGIS Online, but including scenes with the high presence of clouds (50-95%). Fresh scenes with clouds needed because they can contain the areas of interest without clouds.


texts:
What is this error? It clearly says we have a FeatureCollection then it says Cannot flatten an object that is not a collection.

FeatureCollection (Error)
Cannot flatten an object that is not a collection.

Reproducible code:

///////////////////////////////////////////////
///
///                           functions definitions start
///
/// Function to mask clouds using the Sentinel-2 QA band.
function maskS2clouds(image) {
    var qa = image.select('QA60');

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

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

    // Return the masked and scaled data, without the QA bands.
    return image.updateMask(mask).divide(10000)
                          .select("B.*")
                          .copyProperties(image, ["system:time_start"]);
}

////////////////////////////////////////////////
///
/// add Day of Year to an image
///

var addDate_to_image = function(image){
  var doy = image.date().getRelative('day', 'year');
  var doyBand = ee.Image.constant(doy).uint16().rename('doy');
  doyBand = doyBand.updateMask(image.select('B8').mask());

  return image.addBands(doyBand);
};

////////////////////////
///
/// add Day of Year to an imageCollection
///

var addDate_to_collection = function(collec){
  var C = collec.map(addDate_to_image);
  return C;
};

////////////////////////////////////////////////
///
///
/// add NDVI to an image

var addNDVI_to_image = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

////////////////////////////////////////////////
///
///
/// add NDVI to an imageCollection

var add_NDVI_collection = function(image_IC){
  var NDVI_IC = image_IC.map(addNDVI_to_image);
  return NDVI_IC;
};

var extract_sentinel_IC = function(a_feature){
    // var start_date =
    // var end_date =
    // var cloud_percentage =
    // var geom = ee.Feature(feature_col).geometry();
    var geom = a_feature.geometry();
    var imageC = ee.ImageCollection('COPERNICUS/S2')
                .filterDate('2012-01-01', '2019-12-31')
                .filterBounds(geom)
                //.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', "less_than", 10)
                .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 10));

    imageC = imageC.map(maskS2clouds);
    imageC = imageC.select(['B8', 'B4', 'B3', 'B2']);

    imageC = addDate_to_collection(imageC);
    imageC = add_NDVI_collection(imageC);
    //imageC = imageC.map(add_NDVI_collection)
  return imageC;
};

///
///                         functions definitions end
///
///////////////////////////////////////////////

var city = ee.Feature(ee.Geometry.Rectangle(-122.42, 37.78, -122.4, 37.8),
                      {label: 'City'});

var city_IC = extract_sentinel_IC(city);
print (city_IC)
var city_IC_flatten = city_IC.flatten()
print (city_IC_flatten)


texts: Adding latitude/longitude coordinates of a map item in QGIS3 Print composer?

texts:
I can already create lat/lon or UTM grid in Print Composer but I was also asked about possibility to have a mix of both:

1) standard map grid using UTM - I can do this :-)

and

2) latitude/longitude (EPSG:4326) coordinates of map item corners.

The reason is to have UTM grid to easily see the distance on printed map, but also the lat/lon coordinates for a GPS.

I would probably be able to get close to this solution with Map Corners Coordinates Plugin

https://ctu-geoforall-lab.github.io/qgis-map-coords-plugin/

and copy the coordinates from resulting text file in a text object in Print Composer.

But the plugin is still QGIS2 only and the map canvas corners might be quite different from the map item extent in Print composer.

The item properties in composer also displays some coordinates so I might be able to copy it from there but would probably need to switch to EPSG:4326 for map canvas which would distort my maps.

any idea how to do it?

PS: if there is no direct solution, does Print Composer support some sort of custom scripts or plugins we could make someone create?


texts: Logging incoming request HTTP headers in GeoServer

texts:
I have a GeoServer 2.0.2 and I send POST requests using mobile devices. I would like to see how the incoming POST request HTTP headers are formatted, especially what kind of character encoding and response type are provided by the client.

How can I setup GeoServer and/or Jetty to provide me with these informations?


texts:
I'm having lots of trouble with my jpg showing up on my qgis template. I have saved the image in the same folder as my .qgs file, so it should come with associated files. As you can see, the thumbnail is loaded, but only the portion of the thumbnail in the far upper left of the image is transferred to my canvas.

Georeferencing the image does not work either - the image loads in the georeferencer and I'm able to associate 4+ points to my coordinate system on the canvas, but when it converts it to my canvas, nothing shows up. I assume this is because the only part of the image that shows is the tiny portion in the upper left of the image.

This problem is driving me crazy!!! Please help???


texts: I've clicked around on a few of the related question. I am not using qgis or anything "fancier" than ArcMap 10.2.2 for desktop and Microstation V8i. Is there a way to successfully convert shapefiles to a dgn file? I've been able to convert and open a shapefile in a dgn, but the coordinate system is always off. Let me know if more information would be helpful!


texts: How to navigate the REST endpoints and retrieve Field information programatically

texts:
I'm adding a Dynamic Map Service to an ArcGIS Server JavaScript API map, and I need to get hold of each component layer and its fields, so I can run a query.

See the sample ArcGIS Server services directory at http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/HomelandSecurity/operations/MapServer

If you open this URL in a browser you can see the list of available layers (0,1,2) and each layer's fields.

How can I retrieve the list of fields programatically, within a JavaScript function? The code snippet below shows how to retrieve the layer endpoints, but I can't see how to get to the fields.

(One option is to add the layers as Feature Layers, but I'd prefer to avoid this if possible).

var url = "http://sampleserver3.arcgisonline.com/ArcGIS/rest/services/HomelandSecurity/operations/MapServer/";
var dynLayer = new esri.layers.ArcGISDynamicMapServiceLayer(url);
map.addLayer(dynLayer);

dojo.connect(dynLayer, "onLoad", function() {
   var infos = dynLayer.layerInfos;
   for (var i = 0; i <= infos.length - 1; i++) {
       var layerId = infos[i].id;
       var restEndPoint = url + layerId;
       //restEndPoint is the layer's URL - how can I retrieve its fields?
   }
});

Thanks,
Steve (cross-posted to the ArcGIS Server forum)


texts: Duplicating point within shapefile using ArcPy?

texts:
How do you duplicate a point within a single feature class using Python? I have a partial sample script below, but in a nutshell, I want to duplicate a point, copy it into the same shapefile, and populate the fields with user-provided strings. I suspect this is elementary, but I'm new at this and am frustrated with loops and cursors.

# Script tool will take a point file and allow the user to duplicate a single point
# while updating necessary field information. The purpose is not to delete or overwrite
# the original point object, but rather to add a duplicate point with the new information.
# This allows historical context for a point that spatially never changes, but who's
# fields are dynamic. No information is hard-coded.

import arcpy
# set up input variables for target folder and target dataset.
targetFC=arcpy.GetParameterAsText(0)  # this is the point file
targetKeyField=arcpyGetparameterAsText(1)  # this is the key field in the point file
newKeyID=arcpy.GetParameterAsText(2)  # must be duplicate of a string in targetKeyField
newField2 = arcpy.GetParameterAsText(3)
newField3 = arcpy.GetParamaterAsText(4)

# Add new fields to table
in_table = targetFC
field_name2 = "Field_2"
field_name3 = "Field_3"
field_type = "TEXT"
arcpy.AddField_management(in_table, field_name2, field_type)
arcpy.AddField_management(in_table, field_name3, field_type)

# Loop through targetFC to find a match between newKeyID and a string in targetKeyField
with arcpy.da.SearchCursor(targetFC, [targetKeyField, "SHAPE@"]) as cursor:

     # insertCursor based on the match; duplicate the point in the table
          # Populate the new cursor with the new parameters (newField2, newField3, etc)
     # Error if the newParcelID does not match a string in the targetKeyField
     # User is brought back to enter a new newKeyID string.


texts: MODIS image reprojection/ change coordinate system in ArcGIS Desktop?

texts:
I have downloaded Utah state shapefile in NAD_1983_UTM_Zone_12N projection and also downloaded the MODIS image of the area (MOD09Q1) in .hdf format which has sinusoidal projection.

the .hdf format was converted to .img using ERDAS imagine 2011.

Upon adding these two layers in ArcGIS 9.3, the boundary shapefile does not fall on the image rather lies far away. I have also tried projecting the image to NAD_1983_UTM_Zone_12N...its still not working.

Any suggestion on what to do?


texts: Merging all feature classes with same name nested in multiple geodatabases using ArcPy?

texts:
I am using arcpy and I want to extract a specific feature class called "building" nested in a number of geodatabase within my workspace.

Here is the code:

import arcpy
import os

workspace = "C:/Wiley/P1/gis"
search = "Building"
outdir = "C:/Wiley/P1/gis/HK80.gdb"
fc = []

walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Polygon")

for dirpath, dirnames, filenames in walk:
          for filename in filenames:
             if search == filename:
                     fc.append(os.path.join(dirpath, filename))
if fc:
     output = os.path.join(outdir, os.path.basename(search) + "_merge")
     arcpy.Merge_management(fc, output)

This code was successful.

But there are many other featureclass in the geodatabase apart from "Building" So I expanded the code a bit and try to loop the merge command through all these feature class:

    import arcpy
import os

#--------------------------
#example list workspace
arcpy.env.workspace = "C:/Wiley/P1/gis/HKU_Job69.gdb"

#Generate exmaple list dataset
datasets = arcpy.ListDatasets(feature_type='All')
datasets = [''] + datasets if datasets is not None else []
fc_list = []

#--------------------------

#where all the gdb is nested
workspace = "C:/Wiley/P1/gis"

#define output
outdir = "C:/Wiley/P1/gis/HK80.gdb"

#List of Feature classes
fc = []

#define walk for looping
walk = arcpy.da.Walk(workspace, datatype="FeatureClass", type="Any")

for ds in datasets:
     for fc_example in arcpy.ListFeatureClasses(feature_dataset=ds):
        fc_list.append(fc_example)

for f in fc_list: #Loop through the list
   for dirpath, dirnames, filenames in walk:
        for filename in filenames:
           if f == filename:
                     fc.append(os.path.join(dirpath, filename))
   if fc:
        output = os.path.join(outdir, os.path.basename(search) + "_merge")
        arcpy.Merge_management(fc, output)

This is where I hit a wall. Is there anything wrong with this code? I am new to Python so for a trained eye you guys may spot out the error in the logic immediately. The f variable didn't seem to loop at all.


texts: Unable to process "subbasin parameters from raster"

texts:
I was wondering why do i have a different interface for subbasin parameters from raster window. .

i image is displayed below.

instead of this. .


texts:
I have a very good experience with Garmin eTrex Legend HCx -
here and here they write that HCx has a very sensitive chip, that works in forest, deep valleys and in cities - and I can only confirm - with HCx I was able to get my position in plane!

How sensitive is eTrex 30, related to HCx? Have you been able to catch satellites in room / plane / in a narrow street? eTrex 30 doesn't say to be that sensitive at garmin website.

(Could HCx then be more precise/usable in these conditions despite eTrex 30 supporting GLONASS?)


texts:
I'm definitely a 'geonewbie' and I'm trying to test if a point is contained inside a region.
I have taken some certified data here:

https://www.istat.it/it/archivio/124086

getting WGS84 UTM32N for 2016. Then I have imported the Reg_2016_WGS84.shp file with shp2pgsql-gui. The import was successful:

Importing with configuration: reg_2016_wgs84_new, public, geom, [myHome]/Limiti_2016_WGS84/Reg2016_WGS84/Reg_2016_WGS84.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=4326
Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Shapefile import completed.

I'm interested on region 6, so check for validity:

select ST_IsValid(a.geom) from
public.reg_2016_wgs84_new a
where a.cod_reg = 6

Query says: true

Then trying to check a point that for sure is outside the region:

select a.*, St_contains(a.geom, ST_SetSRID(ST_MakePoint(45.462360, 9.185939),4326))
from
public.reg_2016_wgs84_new a
where a.cod_reg = 6

Query says: false

But also this point that for sure is inside the region got false:

SELECT a.*, St_contains(a.geom, ST_SetSRID(ST_MakePoint(46.062340, 13.237270), 4326))
FROM
public.reg_2016_wgs84_new a
where a.cod_reg = 6

So ST_Contains seems to return always false even if the MULTIPOLYGON is valid. May be I'm doing a big mistake somewhere...

My PostGis version:

POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.5.0-CAPI-1.9.0 r4084" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.0.1, released 2015/09/15 GDAL_DATA not found" LIBXML="2.9.4" LIBJSON="0.12" TOPOLOGY RASTER


texts: Leaflet put coordinates of every point in array

texts:
I get points using GeoJSON.AJAX and I want to store coordinates in array to pass it later to function that builds routes between points. Now I use this code

<script>
    function our_layers(map){

    var coords = []; //array to store coordinates
    var points = new L.GeoJSON.AJAX("{% url 'incidences' %}",{

        onEachFeature: function (feature, layer) {
            layer.bindPopup(a lot of code here);

            x = feature.geometry.coordinates
            coords.push(x);//here I add coordinates to array
        }
    });
    points.addTo(map);

    //there is button "getcoords" I created to test some things
    document.getElementById("getcoords").addEventListener("click", logCoords);

    function logCoords() {
         console.log("Test " + coords);
     }

</script>
{% leaflet_map "gis" callback="window.our_layers" %}

That way it works perfect but I need to get coords array without pressing button, obviously. But if try putting concole.log(coords) after points.addTo(map) it returns empty array. I guess it's because it need some time to calculate (because if I press button immediately after refreshing page I get empty array too). Is there any way to print array right after pushing coordinates done and without setting timeout?
I also tried

        function drawPoints(p, callback) {
                    p.addTo(map);
                    callback()
                }
        function logCoords() {
                    console.log("Test " + coords);
                }

        drawPoints(points, logCoords);

But it still giving me empty array. So, is there any way to get my array without setting timeout?


texts: Buffering with physical barrier using ArcGIS Desktop?

texts:
I don't think it is possible, but I'm looking for a way to stop my buffer at a certain physical barrier.

For example, I have point near a river and I'm creating a buffer around my point that is larger than the river-point distance. So, I'm buffering on the other side of the river.

I can erase my buffer polygon with my river polygon and then erase what's on the other side, but I have ~50 sites, with a couple thousands points total so that's a lot of cleaning.

Any ideas?


texts: problem with raster calculator results when subtracting two rasters with different resolutions

texts:
I have two elevation datasets which I am trying to subtract to calculate thickness (using ArcGIS raster calculator). The first datset is a tif which is oriented north-south and has a cell size of 5 ft x 5 ft. The second data set is a surfer grid which has a cell size of  10 ft x 10 ft but is rotated by approximately 30 degrees to the east. When I subtract the two rasters, I select bilinear re-sampling in the environments section of raster calculator. However, when I spot check my resulting output raster of thickness, I often get a thickness that is vastly underestimated (e.g. thickness should be 10 ft but the raster output equals a thickness of 2). Strangely enough, other cells that I spot check will have a thickness that is close to what I expect. I also re-ran the same calculation but used nearest neighbor re-sampling; I got the exact same results.

What is causing this discrepancy? Is it the rotated grid? How can I fix this?


texts: Creating Match Table using ArcPy gives ERROR 000090: Error creating table, workspace may be read-only?

texts:
I am creating an attachment table linking attachments from my local PC, saved in "C:\Temp\Match" to a Featureclass housed in a SDE Database on our GIS Server.

I can run the GenerateAttachmentMatchTable Tool from within ArcGIS and run it multiple times with no issue.
The problem comes in when I try and run the tool from within a Python Script.
I keep getting an Error:

  ExecuteError: ERROR 000090: Error creating table, workspace may be read-only.

The folders are not read-only, I have checked and even changed the location of the output.
The script has been modified and copied from the Geoprocessing Results Window and modified to be acceptable to a stand-alone script.
Please see the code below:

import arcpy
InFeat = "Database Connections/xxx/_Enviro.DBO.MonPnt"
AttFolder = "C:\\temp\\Match1\\JPG"
OutTable = "C:\\temp\\Match1\\MatchTable"
LinkF = "Name"

arcpy.GenerateAttachmentMatchTable_management(in_dataset= InFeat, in_folder=AttFolder, out_match_table=OutTable, in_key_field=LinkF, in_file_filter="", in_use_relative_paths="RELATIVE")

The folder in question is not read-only and in fact has just been created to test this. There are no ArcMap documents linked to this folder, it only contains from images and a few text files I used to run some checks elsewhere in the script.

Running windows 7, 64-bit and ArcMap 10.4.1


texts: Creating large amount of random points in binary raster?

texts:
I want to create a point vector dataset of 10000 points (or larger) within a binary raster, where the points should be constrained to areas where the raster value is 1.

I tried the following steps.

Polygonize raster
QGIS: Vector -> Research Tools -> Random Points

This works fine up to 2000 points but anything above just causes QGIS to crash.

Is there a way to create a vector dataset with a large number of point features constrained by a binary raster (or the polygonized version of it)?

The following tools are at my disposal, ranked from most to least favourable: QGIS, Python, R, ArcGIS

This is what I'm going for, only with 10x the point features.


texts:
Environment: Win10 and GDAL 2.1.0 in the OSGEO4W shell, 64 bit. Edits made below.

I have temporary access to a whole-of-country dataset in ESRI's Arc/Info binary grid format, which according to this - http://www.gdal.org/frmt_various.html - should be supported by gdal. Its accessible through webDAV.

The whole dataset is about 200GB, but I only need a small portion of it, so i would like to use something like

gdal_translate -projwin_src EPSG:4283 -projwin xmin ymin xmax ymax src_dataset=http://user:pass@path_to_file/w1001001.adf dst_dataset=local_path

to get the portion I need. The command fails out with error 4 'failed to create temporary file', because I don't have write access to the remote location.

I've tried incorporating "--config CPL_TMPDIR local_folder" into the command, and I've tried setting it separately (set CPL_TMPDIR=local_path), but that doesn't seem to help much. The next problem is error 4: "local_folder\w1001001.adf is not recognised as a supported file format', which contradicts the documentation.

SOLVED: per comment below, one must use vsicurl, e.g.

gdal_translate --config CPL_TMPDIR D:\local_folder -projwin_srs EPSG:4283 -projwin ulx uly lrx lry /vsicurl/http://user:pass@webdav-url/w001001.adf D:\local_folder\output.tif

works on both tif and adf source files.


texts: Keeping registry window open while executing R script in QGIS?

texts:
I'm trying to execute an R script in QGIS that already worked in another PC. However now it's not working and the window is automatically closed when the processes of the script finishes.

Is there any piece of code to paste at the end of an R script to keep this window in order to see what is been executed?

This is the screen shot that I get if I'm quick enough to catch it. How could I keep this window opened?


texts: How to use GDAL /vsizip/ to access GML file(s) from Sentinel-2 zip package

texts:
Which is the correct way to write the path through the *.Zip S2 imagery package to rich the *.GML vector file(s) ?
Does someone knows how to correctly write the path to *.GML file for using it with /vsizip/ ?
I would like to target the *.GML vector file(s) from *.SAFE/GRATICULE/QI/ subfolder. I did a try but it doesn't work.

Here is also a link to one of my previous possible related questions.

Thanks in advance for any tip!


texts: Why does QGIS try to jam some of labels to end of line geometry and not others?

texts:
My east-west (B) line is labeled as I would expect but the label for the north-south line (A) has the label stuck at the end.  This is the result of the label engine and I think (??) I have tried every combination of settings to try to figure this out.
Does anyone know why QGIS does this and how to remedy it?


texts: SDA4PP not found? How to do Kriging in QGIS?

texts:
I need to make kriging analysis in QGIS (v1.7.4). I found out, that it should be a part of SDA4PP plugin, but I just couldn't find it in any repository (even third party). It looks, like it just disappeared from the web.
Is it possible that I did something wrong?


texts: How to make OpenLayers example work in jsFiddle

texts:
I'm working on a jsFiddle with a OpenLayers example in it. http://dev.openlayers.org/releases/OpenLayers-2.11/examples/snapping.html

At the moment it's not working because it's not getting a response for the http request to get data. How do I fix that?

The jsFiddle is here: http://jsfiddle.net/TcuxA/6/ Go to the line "// create three vector layers" in the script. There are three requests for data. If you type the URLs in your browser you get the JSON, but my Firebug gives 3 errors when I run the jsFiddle.

I tried fixing with jsFiddle echo, but that didn't work. I don't know how to change the script to load the data otherwise.

Why can't I get the JSON by these URLs? What are good solutions?


texts:
This screenshot is an example of what I have to do :

This picture was taken on a different and older version of what I have now, since all of the present images were NOT markers bound to the map

After my images became elements of the map, I had to rewrite this script :

var angle = 0;
var width = 60;
var height = 40;

$(document).keydown(function(e){
    if(e.which == 65){ //"A" key
        angle-=5;
        left();
    };
    if(e.which == 69){ //"E" key
        angle+=5;
        right();
    };
    if(e.which == 83){ //"S" key
        height+=5;
        smaller();
    };
    if(e.which == 90){ //"Z" key
        height-=5;
        bigger();
    };
    if(e.which == 81){ //"Q" key
        width-=5;
        thinner();
    };
    if(e.which == 68){ //"D" key
        width+=5;
        wider();
    };
});

function left(){document.getElementById('rotate').style.transform = "rotate("+angle+"deg)";};
function right(){document.getElementById('rotate').style.transform = "rotate("+angle+"deg)";};
function smaller(){document.getElementById('rotate').style.height =  height + 'px'};
function bigger(){document.getElementById('rotate').style.height = height + 'px'};
function wider(){document.getElementById('rotate').style.width = width + 'px'};
function thinner(){document.getElementById('rotate').style.width = width + 'px'};

and it was working the way I wanted to : having only one image that can rotate, can be resized at a time.

Where as now, I'm using this script (with leaflet.rotatedMarker for the rotation):

var opened = 0;
var angle = 0;
var markerNumber = 0;
var divider = 0;
var markerList = [];
var width = 60;
var height = 40;
var marker = L.marker(coordinates, {icon: L.icon({iconUrl: imagePath, iconSize: [60,40]}),draggable: true}).bindPopup('Use A, E, Z, Q, S or D to adjust me !').addTo(mymap);
divider += 1;
markerList += markerNumber;
$(document).keydown(function(e){
    if(e.which == 65){angle-=5/divider;setAngle();}; //"A" key
    if(e.which == 69){angle+=5/divider;setAngle();}; //"E" key
    if(e.which == 90){height-=5/divider;smaller();}; //"Z" key
    if(e.which == 81){ width-=5/divider;thinner();}; //"Q" key
    if(e.which == 83){height+=5/divider;taller();}; //"S" key
    if(e.which == 68){width+=5/divider;wider();}; //"D" key
});
marker.on('popupopen', function(e){
    opened = 1;
    if(opened == 1){
        $('.leaflet-marker-icon').addClass('adjustable');
        marker.setRotationOrigin('center center');
    };
});
marker.on('popupclose', function(e){
    opened -= 1;
    markerNumber +=1 ;
    if(opened == 0){
        $('.leaflet-marker-icon').removeClass('adjustable');
    };
});

function setAngle(){marker.setRotationAngle(angle)};
function smaller(){$('.adjustable')[markerNumber].style.height = height + 'px'};
function taller(){$('.adjustable')[markerNumber].style.height = height + 'px'};
function wider(){$('.adjustable')[markerNumber].style.width = width + 'px'};
function thinner(){$('.adjustable')[markerNumber].style.width = width + 'px'};
};

I identified two or three problems:

The setRotationAngle impact ALL markers on the map [FIXED]
The getElementsById('rotate') only makes the FIRST marker resizable [FIXED]
after each popupclose :

the width and height vlaue wil automatically be increased by 5 [FIXED]
I can still rotate my markers, which isn't what I want [FIXED]

What I need:

Using popupopen : having only the marker bound to the opened popup resizable and spinnable. DONE
Using popupclose : removing the resizable and spinnable functions of the marker bound to the closed popup. DONE

All I'm trying to do now is to remove the rotationOrigin option when I close the popup

I know this sound extremely familiar to what I've asked before, but now, I'll need to apply this to multiple markers


texts: How to create raster attribute table to a raster data using gdal python?

texts: I classified a remote sensing data and saved it as a raster using Gdal python. I want to create a raster attribute table for this classified result.. my classified raster has five classes with 1,2,3,4,5 representing classes A,B,C,D,E respectively... Is it possible to assign color to each class.. since it is a single band raster image..


texts: Downloading old version of qgis to Windows

texts: I have bought the book "Learning QGIS 2.0" by Anita Graser and wish to work through the tutorials therein in order to improve my knowledge of QGIS.  I currently use QGIS 1.8 and wish to download version 2.0, as I think I would get confused if I were to download the latest (2.6?) version.  What's the easiest way to get hold of version 2.0?


texts: Converting from ECEF to geodetic coordinates

texts:
I am trying to convert an ECEF coordinate to a geodetic coordinate using this procedure,  but I don't fully understand the process. It is prefaced by stating the geodetic parameters {a,b,e,e'} are assumed to be known, but does not state what they reference to.

I assume that a and b are the equatorial and polar semi-axes, and that e is just Euler's number. Am I correct in this? And what is e' supposed to be?


texts: Using VBScript to label shapefile in ArcMap?

texts:
If you're reading this for the first time, it might be best to focus on section beneath ANOTHER UPDATE.

I'm a beginner with VBScript and I have a question about using an VBScript expression for a label. I want to label street address numbers and the fractional address number, if it exists, in my map. This data is stored in an address points shapefile.  I go Properties > Labels > Expression... in Esri software to enter this expression:

[HSE_NO]&" "& [FRAC_NO]

With this expression the labels will be created next to the address points as follows:

 412 1/2 412

See how both labels are on the same line? I would like to have the labels look like this:

412
412 1/2

I want the first label to be the address value without a fraction, and the following lines would be address values with fractions. I welcome your suggestions.

UPDATE

I think what I need is more complex than I realized. As I mentioned, I want first label to be address street value (HSE_NO) without a fraction (where FRAC_NO is NULL). I'm working on getting this to occur with what's below. It's not correct! Any ideas?

Function FindLabel ( [HSE_NO], [FRAC_NO] )
if (IsNull[FRAC_NO]) then
    FindLabel = [HSE_NO]
    End If
End Function

ANOTHER UPDATE

I have a point shapefile of addresses. Some address points contain the address number and street name. Some points contain address number, address fractional number, and street name. Points containing fractional values will overlap points containing not containing fractional values. I want to label these overlapping points like this:

1159
1159 1/2

My current VBScript labels addresses in this manner.

1159 1/2 1159

Does anyone know how to get the desired result? Here's my current expression:

Function FindLabel ( [HSE_NO], [FRAC_NO] )
  if (IsNull([FRAC_NO])) then
    FindLabel = [HSE_NO]
  else
    FindLabel = [HSE_NO] & " " & [FRAC_NO]
  End If
End Function


texts: I have a large raster layer with wind speed values assigned to each cell. I also have point data that displays locations where wind farms have been constructed. I would like to add the raster values to the wind farm location data. The goal is to figure out when each wind farm was built over each wind speed value. Any ideas?


texts:
I am using data driven pages and arcpy.Graph() object to automate production of synchronised map/long-section series similar to this:

I'd like to keep vertical scale for long-section part. In order to do so I have to insert a couple of dummy points at the end of profile, because I cannot access these graph's properties via arcpy:

Does anyone know how to change min/max values for axis using script?


texts: How do I serve data in a sqlite database on a shared web host, like GoDaddy, onto a GIS web map?

texts:
I am new to the web mapping arena.  I currently am managing a website and would like to add a mapping portion to it, mostly for some simple viewing and querying of a small number of layers.
For simplicity, I would like to use SQLite as the data store.  I will use Quantum GIS or some such to load the spatial data into the SQLite database.  Then ftp that onto the web site's directory structure.
I know that I will use OpenLayers to actually develop the map that is integrated into the web page.
I am thinking that FeatureServer would be appropriate to access the data from the SQLite database and convert it to a format that OpenLayers is able to consume.  Is this correct and will it work?  If so, how do I install this on a godaddy type server?  If not, what would be appropriate software to use?
I would like to use Open Source software to accomplish this task.

Addtl Info - 2011/02/02 11:30AM PT
To answer amercader
 - Yes, I can install software.  I also have php and python installed, and can run CGI scrips.

To all
 - It seems that Django is similar to Featureserver, I guess.  I seem to be stuck at where or how to go about installing these programs.  I only know how to access the web site via ftp, I don't know how to access via command line.  Part of my problem is not knowing the correct questions to ask to get the correct answers.
- Perhaps the most direct question to ask is how do I install Featureserver on a shared hosting account like godaddy?

To underdark
 - a text file would be an easy solution, indeed.  I would like to make this scalable for the longer term, as well as having to reduce conversion steps to get my data posted.  In the future, I would potentially like to accept user input that will modify the spatial data being stored, thus making a text file a less desirable option.

Thanks all for input so far, it is appreciated.


texts: ArcGIS Pro recieving error 000405: No records within table when records do in fact exist?

texts:
I selected a subset (blue) of points from the feature class (red) which all have populated fields in the attribute table, and can be manually edited (by clicking each individually and typing a value) but cannot be batch edited by calculating field. When counted using tool, zero records are recognized and geoprocessing error 000405: no records within a table is displayed.

Is anyone else having this problem in ArcGIS Pro?


texts:
I've downloaded grib files from the new Carra reanalysis from ECMWF's CDS. The data covers a really big region in the Arctic.
What I need is to extract a subregion of the whole region to lower filesize and make it more manageble to handle the files. It does not seem to be possible to get a subregion directly from cds so I must use a tool and do it myself.
I've tried to use CDO's sellonlatbox, but it does not support the grid type of the files. I've also tried loading it into python with cfgrib, but xarray can not select region on a curvilinear grid.
I haven't found anything useful in Eccodes either. Is there a simple way to do this?


texts: How to delete polygons less than a certain area?

texts:
I'm trying to delete all polygons with an area smaller than a certain give value. Is this possible in qgis?

The goal here is to delete all the islands off a map.


texts: LinearRegression with Google Earth Engine

texts:
I ve been doing simple linear regressions with GEE and I have found a possible error that I would like to share and in case that I am wrong, please let me know:

Is it possible that when LinearRegressions function is used, GEE always takes equidistant distances in the independent variable(x)?

I ve tried to do regressions with NDVI (y) and time (x) and GEE always takes the same distance between observations and it is not true because between images(dates) there is no the same distance, sometimes are 1 day but sometimes are 3, 10 or whatever.

Does anyone know how to introduce the real distance between independent variable?

For example: image1 (from 1st of January), image2 (from 2nd of January), image3 (from 10 of February) would be: 1-2-10 instead of 1-2-3 (what GEE does).

The code is as follows:

// Simple regression of year versus SAVI.

// Define the start date and position to get images covering TDPS,
// Bolivia, from 2000-2010.
var start = '2000-01-01';
var end = '2014-12-31';

var region = table4

// Filter to Collection
// time of year to avoid seasonal affects, and for each image create the bands
// we will regress on:
// 1. A 1, so the resulting array has a column of ones to capture the offset.
// 2. Fractional year past 2000-01-01.
// 3. SAVI.

var col = ee.ImageCollection("MODIS/MOD13Q1")
  .filterDate(start, end)
  .filter(ee.Filter.calendarRange(1,1,'month'))
  .filterBounds(region)
  .map(function(image) {
    var date = ee.Date(image.get('system:time_start'));
    var yearOffset = date.difference(ee.Date(start), 'year');
    var savi = image.expression(
      '(1 + L) * float(nir - red)/ (nir + red + L)',
      {
        'nir': image.select('sur_refl_b02'),
        'red': image.select('sur_refl_b01'),
        'L': 0.5
      });
    return ee.Image(1).addBands(yearOffset).addBands(savi).toDouble();
  });

print(col)
// Convert to an array. Give the axes names for more readable code.
var array = col.toArray();
var imageAxis = 0;
var bandAxis = 1;

// Slice off the year and savi, and solve for the coefficients.
var x = array.arraySlice(bandAxis, 0, 2);
var y = array.arraySlice(bandAxis, 2);
var fit = x.matrixSolve(y).clip(region);

// Get the coefficient for the year, effectively the slope of the long-term
// SAVI trend.
var slope = fit.arrayGet([1, 0]).clip(region);

Map.addLayer(slope);

// Export the image, specifying scale and region.
Export.image.toDrive({
  image: slope,
  description: 'SAVI_slope_jan_MODIS_00_14',
  scale: 250,
  maxPixels: 1e9,
  region: geometry
});


texts: What's the most efficient way to load feature attributes into a numpy array?

texts:
I am currently working on an add-on that relies on the user choosing fields of tables (that can be from shapefiles, csv or dbf) that are open in QGIS and then making computations with such data.

Right now I am looping through all the features and collecting the data using the standard cookbook techniques to loop through features and to read attributes, but it is taking too long.  Some tables have about 25,000,000 features, and the actual data loading takes much longer than the computing. Is there any shortcut?  Can I load whole vectors at once or some other similar thing?

It occurs to me that I can just get the file information and go directly to the database in disk and load it from there, but I REALLY want to avoid it.

It seems that the question was not very clear, so here goes what I have now.

all_matrix=[]
for f in matrix.getFeatures():
    a=feat.attributes()[origin]
    b=feat.attributes()[destination]
    c=feat.attributes()[flow]
    all_matrix.append([a,b,c])
 matrix=np.array(all_matrix)


texts: Changing particular value of raster dataset to zero using ArcPy?

texts:
I try to change a particular value of  raster dataset to zero using python script.
When I deal with one raster, it works. Python scripts as follow:

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
# Set environment settings
env.workspace = r'D:\MATLAB\TEST'
# Set local variables
a="states1166.tif"
inRaster = Raster(a)
outCon = Con(inRaster == 1.000000020040877e+20, 0, inRaster)
# Save the outputs
outCon.save("D:\MATLAB\TEST\prif11166.tif")

But when I deal with multiple raster, it failed. Python scripts as follow:

# Import system modules
  import arcpy
  from arcpy import env
  from arcpy.sa import *
  import sys
  import os
  arcpy.CheckOutExtension("Spatial")

  env.workspace = r'D:\MATLAB\PriForest'

  outdir = os.path.join(os.path.dirname(env.workspace),"Replace")
  if not os.path.exists(outdir) : os.mkdir(outdir)
  inRasters = arcpy.ListRasters("*")

  print "Begin Replace"

  for inRaster in inRasters:

      outCon = Con(inRaster == 1.000000020040877e+20, 0, inRaster)

      filename = 'PriForest' + os.path.basename(inRaster)
      outfile = os.path.join(outdir,filename)

      outCon.save(outfile)

   print('finished')[enter link description here][1]

Does anyone know the reason?


texts:
I have a problem with QGIS as here Previous question.  I think QGIS Temporal is corrupted and the only way to correct is to uninstall and reinstall QGIS (QGIS 3.16.1 Ubuntu 20.04). However if I try this via Ubuntu Software I get the error

Unable to remove "QGIS Desktop: The following packages have unmet
dependencies

but no packages are listed.  Is there a way to remove QGIS or reinstall/repair via Terminal please?
It would be nice if people were intent on answering questions rather than nitpicking wording


texts: How do I create a negative to positive value raster with neutral 0?

texts:
I have a raster layer with a band with values from -137 to 9229. I want to create a color scheme showing how far from 0 the positive numbers are in one color, and how far from zero the negative numbers are in another color, with a neutral 0 (that, unmistakenly, doesn't lean towards the positive or negative colors, nor diverges from the color scale).

Moreover I'd like a "-55"-area to be just as far from the neutral color (0) as a "+55"-area, but in different colors. So if neutral would be #333, then -5 might be #338, and +5 would be #833, even though the array of values, as in this case, is heavily biased towards positive numbers.

Hope that explains my problem. I'd greatly appreciate any help I get! This is for my Bachelor thesis project on the creation of the metro and highway system in DC during the 50s to 80s.

EDIT: Oh, gosh. I forgot to tell you I'm using QGIS 1.7.4 on an OSX Lion. And have very limited access to Cartographica and ArcGIS full suite.

Best regards,
Henrik


texts: Generate or create OGC WMS using Python mapnik

texts:
I want to create an OGC WMS server from mapnik PostGIS using Python.
I have successfully post table on PostGIS using:
shp2pgsql -s 4326 /Users/apple/Downloads/111SHAPEFILE/Point.shp points | psql -d gis -U apple -h localhost

then I successfully connect mapnik with PostGIS database using:
from mapnik import PostGIS
postgis = PostGIS(dbname="gis", user = "apple", host = "localhost", port = '' ,password = "", table= 'points')

I want to create OGC WMS server from mapnik, how can I do this?


texts:
I'm fairly new to ArcGIS and I'm trying to string along several queries. individually they work in same query they don't.
background

I want to query all permit types  rc, re, re/ap, and db; less than 30 days old, and only in certain counties.

I'll be making 3 layers with ranging dates <30, 30-60, 60-90, >90

whats working
Layer 1:

"PERMITTYPE" ='RC' OR "PERMITTYPE" = 'RE' OR "PERMITTYPE" = 'RE/AP' OR "PERMITTYPE" = 'DB' AND "DATEFILED" >= CURRENT_DATE - 30

** This returns about 60 results in the below counties

Layer 2:

"COUNTY" IN ( 'BORDEN (TX)', 'GARZA (TX)', 'LYNN (TX)', 'KENT (TX)', 'SCURRY (TX)' )

Layer 3 (looking for different permit type but works)

"PERMITTYPE" = 'OD' AND DATEFILED >= CURRENT_DATE - 30 AND "COUNTY" IN ( 'BORDEN (TX)', 'GARZA (TX)', 'LYNN (TX)', 'KENT (TX)', 'SCURRY (TX)' )

New layer (not working )

"PERMITTYPE" ='RC' OR "PERMITTYPE" = 'RE' OR "PERMITTYPE" = 'RE/AP' OR "PERMITTYPE" = 'DB' AND "DATEFILED" >= CURRENT_DATE - 30 AND "COUNTY" IN ( 'BORDEN (TX)', 'GARZA (TX)', 'LYNN (TX)', 'KENT (TX)', 'SCURRY (TX)' )


texts: Creating polygons from line segments using PostgreSQL and PostGIS

texts:
Is there a PostGIS function that creates polygons from multiple line segments?

I want to convert the pictured thirteen lines to four polygons. All lines belong to the table "boundary" (gid, geom). I've categorized them by the gid column for a better illustration.

Is ST_MakePolygon what I'm looking for? Can anyone help me with the SQL syntax that QGIS requires for loading the polygons?

Edit: changed "closed line geometries" to "multiple line segments"


texts: Zonal statistic that fits area portion of polygons using ArcMap

texts:
As you can see below, I have a polygon feature with boundary of seven regions. I need to extract the raster value intersect to it that fit portion of every region area touches to the pixel. For example, region 4 touches some % area of pixel f and the rest % area of pixel g.
So the final results value I want for region 4 is equal to (% area of pixel f x raster value of pixel f) + (% area of pixel g x raster value of pixel g).
I'm using ArcMap.


texts: Debugging ArcMap Add-In .Net 4.5 in Visual Studio 2013: Unable to start ArcMap.exe

texts:
I am developping an ArcMap 10.3 Add-In targeting .NET framework 4.5 on Visual Studio 2013. I used ESRI's wizard to build a simple button add-in and the projects debug properties are configured to start the external program "C:\Program Files (x86)\ArcGIS\Desktop10.3\bin\ArcMap.exe".

The add-in ran properly when I started debuggig the Add-In for the first time but I couldn't get the breakpoints to hit. I followed ESRI recommandation "How to debug" on ESRI online help and modified the ArcMap.exe.config file to enable supported runtime 4.0.30319

<startup>
     <supportedRuntime version="v4.0.30319"/>
      <!-- <supportedRuntime version="v2.0.50727"/> -->
</startup>

The next time I started debuggin the Add-In I get an error message from Visual Studio : "Error while trying to run project: Unable to start program 'C:\Program Files (x86)\Desktop10.3\bin\ArcMap.exe'."

This error message comes right after ArcMap's splash screen. The following screenshot is not from the actual Add-in because the splash screen appears and disappears way too fast and I can't get screenshot.

I tried different things but without success :

build the project using different platforms :"Any CPU" and "x86"
changed "ArcCatalog.exe.config" file to use the same supported runtime
restarted computer
used different supported runtime version: version="v4.0.30319", version="v4.0", version="v4.5"

I dont know what to do next to be able to debug the Add-in. I don't want to develop without this feature!


texts:
I want to add mouse on-hover effect when I am about to select a particular property Lot (section of land). E.g. highlight and change its background color. Is this even possible and where to modifiy the code exported from qgis2web package?

I found a style property variable inside [layer_name]_style.js under 'style foler' but not sure if this is supposed to be the one to change.


texts: Why is the use of memory mapped buffers in shapefile data stores not recommended on Windows?

texts:
When creating a shapefile data store, you have the option to enable 'memory mapped buffers'. Both GeoServer and GeoTools recommend to disable this on Windows.

GeoServer docs:

  Enables the use of memory mapped I/O, improving caching of the file in memory. Turn off on Windows servers.

GeoTools docs:

  Optional: memory map the files (not recommended for large files under windows, defaults to false)

What is the reason for this recommendation? The only seemingly related resource I found is a discussion from 2004, but I don't know if it still applies or if it is actually related to the mentioned documentation.

  I've discovered why the shapefile data store test fails on windows.
  Basically, the datastore is closing a reading channel on a file in
  order to overwrite it with a new file that was build while using
  the writer. Now, this triggers an exception because you cannot open
  a file for writing that still has a memory mapper part, and the
  memory mapped part is still there by contract (sigh!).

  [...]


texts: How is StreetDirectory.com web map application developed?

texts:
I wonder can you tell me how is developed map application like on this site? http://www.streetdirectory.com

What API is used?


texts: How to show up label more than one properties of a shapefile

texts: How to show a label that includes more than one properties of a shapefile on the layers, let say I'd like to show up to properties of a shapefile in ArcInfo / ArcGIS 10.0.


texts: Splitting polyline at multiple points using Mapinfo Pro?

texts:
I have a table of points and used the mapbasic app "Connect the Dots" to create a polyline.   The problem is it creates a single polyline and I want to split it up into separate lines so I can delete a few line segments.

Is there a way to take a single polyline and split it up using multiple points as the breaks?  (Mapcad did not have what I needed.)

Backstory: I am using this to make lines of GPS tracking data and sometimes we get a bad GPS read.  That bad read skews the line away then back.   Deleting the point can be problematic since we'd need to run the tool then delete the problems and the line created then run the line again.  I'd prefer to split the line up and select and delete the segments as needed.


texts: Datum transformation with velocity corrections

texts:
I have a coordinates with velocity values of some points in ITRF2014 at 2010.0 epoch and I want to convert these coordinates to ITRF96 1998.0 epoch.
I apply the velocity corrections to my coordinates before apply transformation to ITRF96.
X(t) = X(t0) + Vx * (t-t0) % to X coordinate , x(t0) i ITRF2014

I used the formula as above, but I'm not sure if it's correct.
I used MJD values and then divided the (t-t0) by 365.25:
t = 1998 ---  MJD = 50814
t0 = 2010 --- MJD = 55197

(T-T0) = (50814-55197)/365.25

Is it correct to do the velocity correction this way? Is my use of time correct?


texts: Calculating azimuth from two points, both having latitude/longtiude?

texts:
How to calculate the Azimuth angle between two points on the surface of the earth if i ONLY have the lat/long of them?

Here is an example image: the azimuth angle is X which i want to calculate AND the angle to the magnetic pole.


texts: Odd Color Shift after Merging Rasters in QGIS

texts:
I have a large number of raster orthophotos that I'm trying to merge into a single image using QGIS. The merge completes but the color is significantly off afterwards, as shown here:

This seems like it should be a very simple operation, so I must be missing some key piece of information.

I've placed two of the original orthophotos here in case someone wants to try this and help me figure it out:

https://www.dropbox.com/s/mnssxgmb2jfhrj3/Aerials.zip?dl=0


texts: Choose between two points with same name but different attributes in ArcGIS

texts: I have a dataset of thousands of points, made by the intersection of flowlines with waterbody polygons. Each polygon has two (or sometimes more, if the waterbody has multiple tributaries) points made by these intersections. I would like to select the point with the lowest elevation, for each polygon. If there an automated way to do this?


texts: Changing ArcSDE version of *.mxd in ArcPy?

texts:
Is there a way to switch the version a .mxd data source is pointing at without using arcpy.mapping.findAndReplaceWorkspacePaths?

I have a script that automates map book production that works very nicely. The source .mxd used by this script points as SDE.Default. Occasionally editors want to see how their edits will look in the final product before they rec and post changes to SDE.Default. To do this I want to call the existing script but have it's .mxd point to the same version that the editor is currently working on. I can't use the editors .mxd because it will not have the correct layout elements and data driven pages.

The only way I can think of to do this is to pragmatically create a new .sde connection file pointing at the edit version and use arcpy.mapping.findAndReplaceWorkspacePaths with this new file to change the source.

But it would be easier and cleaner to just switch version on the map layers. Is this possible? Or is there another approach that I'm missing.


texts: Converting osm.pbf to shp-format using ogr2ogr to extract major highways

texts:
I want to extract all major roads (motorway, trunk, primary, secondary) of the world from OSM data. So far I downloaded the continent.osm.pbf files. Now I need to convert them into .shp-format.
Searching for a way to do so I found the ogr2ogr-tool and so far, it extracts all lines to a sqlite database.
ogr2ogr -f SQLite continent.sqlite continent.osm.pbf lines
Now I also wanted to edit the osmconf.ini file, so that it only includes highways with the tag: motorway, trunk, primary, secondary. I deleted other keys such as waterway,aerialway from the file, but it does not exclude them from the sqlite-database.
# keys to report as OGR fields
attributes=name,highway

Further, I just want highways with the tags "motorway, trunk, primary, secondary", but the only line where it seems to be possible to select them is the following line, which I do not understand fully
z_order_sql="SELECT (CASE [highway] WHEN 'minor' THEN 3 WHEN 'road' THEN 3 WHEN 'unclassified' THEN 3 WHEN 'residential' THEN 3 WHEN 'tertiary_link' THEN 4 WHEN 'tertiary' THEN 4 WHEN 'secondary_link' THEN 6 WHEN 'secondary' THEN 6 WHEN 'primary_link' THEN 7 WHEN 'primary' THEN 7 WHEN 'trunk_link' THEN 8 WHEN 'trunk' THEN 8 WHEN 'motorway_link' THEN 9 WHEN 'motorway' THEN 9 ELSE 0 END) + (CASE WHEN [bridge] IN ('yes', 'true', '1') THEN 10 ELSE 0 END) + (CASE WHEN [tunnel] IN ('yes', 'true', '1') THEN -10 ELSE 0 END) + (CASE WHEN [railway] IS NOT NULL THEN 5 ELSE 0 END) + (CASE WHEN [layer] IS NOT NULL THEN 10 * CAST([layer] AS INTEGER) ELSE 0 END)"

I found the wiki.openstreetmap blog entry, describing how to edit the osmconf.ini, but it also did not clarify things for me.
So my question is: Can I edit the osmconf.ini file in such a way that it only extracts major roads from the OSM-file or do I have to do this in a separate step in my sqlite-database and export it just then to QGIS?


texts: Using processing.runalg in python script

texts:
I'm trying to run the following python script:

from qgis.core import *
import sys
import os
import qgis.utils
from qgis.analysis import *

app = QgsApplication([],True,"")

# supply path to where is your qgis installed
QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/MacOS/", True)

# load providers
QgsApplication.initQgis()

sys.path.append("/Applications/QGIS.app/Contents/Resources/python/plugins")
import processing

layer3 = QgsVectorLayer("../output/precinctsZipcodes.shp", "intersect", "ogr")
if layer3.isValid():
    print "Loaded layer 3"
else:
    print "Layer 3 failed to load!"

processing.runalg("qgis:exportaddgeometrycolumns", layer3, 0, "../output/precinctsZipcodesArea.shp")

QgsApplication.exitQgis()

However, the line processing.runalg("qgis:exportaddgeometrycolumns", layer3, 0, "../output/precinctsZipcodesArea.shp") isn't saving any output. The strange thing is that this syntax works when I run it directly in the python console. Does anybody know what might be happening?


texts: Geoserver GET request to the REST API without authentication

texts:
I am using Geoserver to serve a number of WMS layers.  I am using jQuery's $.getJson('/geoserver/rest//featuretype.json') to programmatically populate a table with the available layers and their abstracts (Geoserver REST API).  I am also iterating over a Geoserver reflect request to generate thumbnails.  At this point, this requires authentication.

Is it possible to access this information through the REST interface without authentication?

Should I be making these request through PHP, Python, or cURL, instead of directly calling the URL in the JS?

I understand the need to authenticate some of the REST URLs as they provide the login and password for the postgis backend in plain text.


texts: How to create my own remote sensing index?

texts: Using Landsat images how do I create my own indices to extract particular landuse? Are there any tool/ software for creating indices?


texts: Create a new connection in QGIS 3.4.3 DB Manager

texts: I am trying to establish a new database connection from QGIS 3.4.3 to Postgres/PostGIS. I have two connections established but I can't find any options to establish another connection. I cannot recall how the other connections were established and I can't find any information on how to do this.


texts: Problems moving POSTGIS database from postgresql 8.3 to 9.3 using Navicat

texts:
I am using Navicat 11 to move a POSTGIS database from 8.3 to 9.3 on Mac OS X 10.9. I have installed PostgreSQL 9.3 and PostGIS with homebrew (brew install postgresql postgis) and everything looks ok. I then create the db locally with template "postgis_template".

Then with Navicat I connect to the remote 8.3 PGSQL, choose the remote db and choose "Data Transfer", choose the remote db, schema public and the local db, schema public. When the transfer starts I get this error:

[Msg] [Dtf] Drop table: "geometry_columns"
[Err] [Dtf] Transfer Data [DROP TABLE IF EXISTS "geometry_columns";]:
ERROR:  "geometry_columns" is not a table
HINT:  Use DROP VIEW to remove a view.
[Err] [Dtf] Terminated

In PGSQL 8.3 there is a table "geometry_columns" while in 9.3 there is a VIEW with the same name. If I try to delete the VIEW in my new db I get an error that it is need by POSTGIS.

Also after several tries I was able to delete the view but now I get this error and wonder if everything is transfered and if I have to do something:

[Msg] [Dtf] End transaction on target server
[Err] [Dtf] Transfer Data [ALTER TABLE "kulturminner" ADD CONSTRAINT "enforce_dims_punkt_geom" CHECK ((ndims(punkt_geom) = 2)) NOT DEFERRABLE INITIALLY IMMEDIATE;]:
ERROR:  function ndims(geometry) does not exist
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.
[Err] [Dtf] Terminated

Any idea how to fix this?


texts: Not working formula in Query Builder (Filter)

texts:
I have been trying to filter a big dataset based on multiple values. The idea is to extend the formula as put in the image below. But even the first part doesn't work.
When I do: "kadastraleGemeenteWaarde" IN 'Ridderkerk' I get all results. If I do "kadastraleGemeenteWaarde" = 'Ridderkerk' I get no results.
A tip in the right direction would be appreciated a lot.


texts:
I'm trying to use a for loop in a list to assign different buffer widths to the list depending on the length of the line being buffered.  I have my list of lengths, and tried to use a for loop to iterate the list with an output variable to hold the new buffer width values.  When I say print buffer_width, the result is "300 Feet" so the for loop is not going though every item in the list, just printing the last item. Finally, I'd like to place my buffer_width variable into the buffer_analysis() function.

length = [1990, 2016, 3051, 3279, 1855, 3432, 4323, 3016, 4218]

for le in length:
    if le < 3000:
        buffer_width = "75 Feet"
    elif le > 4000:
        buffer_width = "300 Feet"
    else:
        buffer_width = "175 Feet"

print buffer_width


texts:
I followed the instructions here. My input layers are 7 .tif images, each representing a spectral band of a Landsat image. I'm attempting to extract these bands using GDAL:

raster_dataset = gdal.Open(raster_data_path, gdal.GA_ReadOnly)

bands_data = []
for b in range(1, raster_dataset.RasterCount+1):
    band = raster_dataset.GetRasterBand(b)
    bands_data.append(band.ReadAsArray())

Issue is, raster_dataset.RasterCount == 1 and I'm getting arrays full of zeros.


texts: Unchecking "Scale symbols when reference layer set" by default in ArcMap?

texts:
In Layer Properties > Display, the "Scale symbols when a reference layer is set" is checked by default, and 90% of the time I don't want/need the symbols scaled.

Is there any way to change the default setting so this box is unchecked when a new layer is added?

I am using ArcGIS Desktop 10.3.1.


texts: Android Google Maps average polyline paths

texts: I'm making an android app that maps my path through plotting polylines using the google map API. If I have several of these paths, it just looks like a bunch of paths plotted on top of each other. Is there a way, through the Google Maps API to "average" these paths out? I was starting to head down a complicated path of computing the headings (geography util) and using those headings to get nearby polylines with similar headings and then averaging out the lats and lons to draw a polyline that was the average of all of the similar polylines. Before I do that I was wondering if anyone knew of a better way and/or if I should be using polygons instead of polylines? Here's a visual of what I am trying to do (sorry I used paint):


texts: How to re-implement a Dojo application using Leaflet?

texts: I apologize in advance if this is not the correct forum to ask my question. I have an web map application developed using the Dojo framework. However I would like to create the same web map application using leaflet. I see that Dojo uses widgets, does leaflet API have something similar to the Dojo widgets that I can use? Is there an example that I can find where someone is using Dojo and leaflet together?


texts:
I am trying to develop a java application that displays : Google maps, Openstreetmap and Bing Map Data.

So which is the best SDK to use in this contet Geotools or uDig ?


texts: Add GML Vector Layer to OpenLayers from String (Variable)

texts:
I am trying to add a vector layer to a map, wich is working fine if included like this:

var gmlfeature = new OpenLayers.Layer.Vector("GML", {
                      protocol: new OpenLayers.Protocol.HTTP({
                        url: "gmltest.xml",
                        format: new OpenLayers.Format.GML()
                      }),
                      strategies: [new OpenLayers.Strategy.Fixed()]
})

map.addLayer( gmlfeature );

I wanted to be able to upload .gml Files to the webpage, though. After parsing the GML content and reprojecting the coordinates, I would like to include the features from a string (wich has the same content as 'gmltest.xml). I tried to follow the method described in this post with no luck.
Here is what I've come up with so far: jsFiddle, the Features seem to load but they are not projected as expected.
A sample GML can be downloaded here.


texts: How to find the closest Edge from a Line

texts:
Currently I'm working on a basic 3D model viewer in MATLAB, which can import .OBJ files (modeled in e.g. Blender). The purpose of this viewer is to test different subdivision schemes (like Catmull-Clark, or Doo-Sabin). A subdivision scheme is basically a method to convert a mesh to a smoother (refined) mesh.

In order to do this properly, I would like to implement some functionality to select Vertices and Edges in my viewer.

Vertex selection isn't very hard -- the way this works is that when I click somewhere in the viewer, MATLAB provides me with two points, the begin- and end points of the line perpendicular to the screen (by clicking on a 2D screen that depicts a 3D scene, you actually click on a line, not a single point). Using this method it is relatively simple to find the closest Vertex.

However, I'm not quite sure how to do the same for an Edge. Imagine that the imported geometry is a simple cube, so 8 Vertices and 12 Edges. If I know click on an arbitrary point in the viewer (i.e. on some line perpendicular to the screen), how can I then determine which of the 12 Edges is the closest one?

[First Edit] Just to be clear, I'm not really looking for a piece of (MATLAB) code, but rather for some theoretical explanation (or just a reference to an article/book/website). Thanks!

[Second Edit] By the way, is this the correct place for such a question? Or might the Mathematics or StackOverflow boards be better suited?


texts: Clipping a buffer from network analyst shapefile to a base shapefile

texts:
I have exported my polygon shapefile which I got after I ran my Network Analyst into a shapefile. It has 6 buffer zones (delimited by time in increments of 10 minutes). When I check the attribute table it has my FromBreak, ToBreak and all the increments, etc.

Now I want to clip this buffer shapefile I have exported to my base shapefile (which is basically a municipality shapefile). Everything that is not included in this municipality shapefile, it's water.

Obviously I do not want my time buffer to run over water, however, when I clip the buffer to my municipality shapefile, the result shapefile loses it's buffer increments. Actually if I check the attribute table, it has the same attribute table my base shapefile had (the municipaility shapefile). How do I make the resulting shapefile after clipping to keep my buffer shapefile attributes instead?


texts: Grid tracing not wanted when using Trace tool for editing in ArcGIS for Desktop?

texts:
So I am using the trace tool to complete polygons. The problem I am running into is that the trace tool will run along my grid lines. This is easily avoided if the polygon is easily distinguishable from the grid line but if my polygon crosses a grid line at say 1-5 meters and I am collecting/creating polygons at 1:5,000 then theres no way for me to know that I just created the polygon along the grid line instead of the polygon. Does anyone know of a way to prevent the trace tool from tracing on my grid layer/layers. I have converted the grids to line work (to use production difference polygon) and then thought maybe converting to graphics would solve the problem but it did not. I also made the grids layer unselectable. Any suggestions?

The lime green line is the grid and the dark yellow is the polygon.

This is one that I probably would catch but...


texts: I have attached my stream network where the green dots shows the location with discharge gage now I want to create a raster all along the stream network to have a discharge data all along the river. Which Interpolation will be better. I have tried with IDW method but it's not working well. Please provide me some suggestion.


texts: Geoserver CSS line with different colors

texts: I want to visualize a line with Geoserver with multiple colors with CSS. Each 3 px, another color has to show. So I have red, green and yellow and I want a line with this 3 colors, each 3 px color should change. Is this possible with CSS?


texts: Seeking public transport data for Berlin in csv or shp?

texts:
I am looking for a shapefile or an Excel document (csv) with all the public transport stops (bus, underground, tram and train). I have searched on OpenStreetMap but I have problems downloading the data.

Do you have any suggestions for me or know where I could find such data?

When I go to https://overpass-turbo.eu/ and enter for example:

node
  [operator=BVG]
  ({{bbox}});
out;

node
  [operator="S-Bahn Berlin GmbH"]
  ({{bbox}});
out;

I get some results, but let us say 70% of the stations are missing. And the weirdest part is that I'm getting only amenity=vending_machine, not the "stops".

To sum up, my goal is to get a csv file, that contains:

Name of the station
What kind (U-Bahn [metro], S-Bahn, Bus or Tram)
Lat/Lng

I found this wiki article about the Berlin Transportation, but it is not complete.


texts: Changing dataSource URI of WMS RasterLayer

texts:
I'm trying to change the data source of a QgsRasterLayer to point to a different URI.

I can change the dataProvider's datasource with wmsLayer.dataProvider().setDataSourceUri(uri) and the change is correctly shown.

But as it's just changing the dataProvider the change is not saved and will default back on reloading the project. wmslayer.source() also still returns the old source.

There is a QgsRasterLayer::setDataSource() function, but it doesn't appear to have a python binding.

How would one go about changing it correctly?


texts: Calculating Area for data in WGS 84 using ArcGIS for Desktop?

texts:
I am using ArcGIS 10.1.1 for Desktop on my laptop.

When I try to calculate the area in WGS 1984 the area is disabled.

If I work in one projection, say WGS 1984, and calculate area in another Projection, is it ok to overcome this problem or it will effect the area calculation?

Can I calculate area using other projection and working in one projection?


texts: Land Classification with Spectral Python (SPy)

texts:
Does anyone have experience working with the Spectral Python module for land classification?  The following k-means clustering example takes an iterative approach to classification through pixel clustering.  Any thoughts, examples and/or visuals from your own work with SPy would be appreciated.

Source: SPy

>>> (m, c) = kmeans(img, 20, 50)
Using single-pass cluster algorithm to initialize clusters.
Clustering image...done
Starting iterations.
Iteration 1...done
Iteration 2...done
        18743 pixels reassigned.
Iteration 3...done
        3387 pixels reassigned.
Iteration 4...done
        2512 pixels reassigned.
Iteration 5...done
        1960 pixels reassigned.
---// snip //---
Iteration 28...done
        149 pixels reassigned.
Iteration 29...done
        135 pixels reassigned.
Iteration 30...done
        107 pixels reassigned.
^CIteration 31... 15.9%KeyboardInterrupt: Returning clusters from previous iteration
>>>


texts:
I need to build an application with the following requirement:

shows a map .png picture to user.
then, users select several points around the map (I can get the x/y pixels from each point)
based on these points I want to convert x/y pixels to latitude/longitude values.

The question is: how can I do this? what do I need to consider?

I think that I can relate a point (x,y) with a predefined lat/lon and then calculate the other points based on the first one. But I'm not sure how can I achieve this.

This picture shows more or less what I need.

green point: point with predefined lat/lon
red points: points I need to determinate


texts: Python to get ids from shapefile don't return all records

texts:
I have some Python code to get the ids and georeference from a shapefile:

canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()
nos = QgsFeature()
for i in allLayers:
    if i.name()=='path_nodes':
        path_nodes = i

nodedata  = path_nodes.dataProvider()

while nodedata.nextFeature(nos):
    i = i + 1
    print str(no.gid())

It does not return all of the ids. Can anyone help?


texts: Disable button on User Form upon opening MXD?

texts:
I have a toolbar with two buttons that allow users to display an edit attribute window of a selected feature.  I would like the buttons to be disabled until the user selects one feature.

Current, I somewhat have this feature working. The button is enable whened opening the document and then it is disabled after the first time opening the window (I have it in the OnUpdate Event). How can I change this so that the button is disabled when the user opens the MXD and then is enable when only one feature is selected?

I think my problem is just that i have it in the OnUpdate event?? Using Visual Basic 2008 Express, ArcGIS 10

Protected Overrides Sub OnUpdate()

Dim pFeatLyr As IFeatureLayer
pFeatLyr = GetLayerByTOC("AJ_RecommendedImprovementsLines")

'Enable Button
If GetFeatSelCount(pFeatLyr) = 1 Then
    Me.Enabled = True
Else
    Me.Enabled = False
End If
End Sub


texts: GEE: Creating a buffer stratum to improve accuracy

texts:
I want to create a buffer (post-)stratum in addition to my other four strata (forest, nonForest, degradation, deforestation), in order to increase the accuracy of the auto-classification. My wish is to re-classifiy pixels into this buffer stratum, if they are both 1) classified in the 'forest' or 'nonForest' strata and 2) in a 2-pixel (60 meters) radius from a disturbance event (classified as 'degradation' or 'deforestation'). Has anyone an idea how I could do this? Any tips on creating a buffer stratum in general would be very helpful already.
Here the code:
// The results are masked in areas of non change, so first unmask
var unmasked = image.unmask()

// Get each stratum and assign correct label
var forest = unmasked.select('forest_2004').eq(1)
  .and(unmasked.select('dist_1').eq(0))

var nonForest = unmasked.select('forest_2004').neq(1)
  .and(unmasked.select('dist_1').eq(0)).remap([1],[2])

  var magThreshold = 20

var degradation = unmasked.select('post_1').eq(1)
  .or(unmasked.select('mag_1').lte(magThreshold))
  .remap([1],[3])

var deforestation = unmasked.select('post_1').gt(1)
  .or(unmasked.select('mag_1').gt(magThreshold))
  .remap([1],[4])

// Combine the strata into a single band image
var strata = ee.Image.cat([forest,nonForest,degradation,deforestation])
  .selfMask()
  .reduce(ee.Reducer.firstNonNull())

Map.addLayer(strata, {min: 1, max: 4, palette: ['#013220','#e5d8bd','#08e8de','#e80812']}, 'strata')

var attributes = {
  forest: 1,
  nonForest: 2,
  degradation: 3,
  deforestation: 4,
  years: '2004-2020'
}

var geo = image.geometry()
Map.addLayer(geo,{},'geo')

Export.image.toAsset({
  image: strata.toByte().setMulti(attributes),
  scale: 30,
  region: image.geometry(),
  maxPixels: 1e13,
  description: '1501_strata'
})

Following the script in the GEE:
https://code.earthengine.google.com/eaea56da0d7de873b8e4fe11d61ad20a?accept_repo=users%2Fbullocke%2Fcoded


texts:
I am writing a web service which in return gives json format data can anyone help me to convert this data into geojson format.

{
 "STATE_ABBR": "AK",
"file_name": "addresses",
"AREA_WATER": "2.45383480336E11",
"AREA_LAND": "1.47795321158E12",
"loc_address": "Kenai Spur Hwy, Kenai, AK 99611, USA",
"STATE_NAME": "Alaska",
"PERSONS": "710231",
"address_type": "church",
"MALE": "369628",
"location_type": "GEOMETRIC_CENTER",
"loc_lng": "-151.1820715",
"loc_name": "Kenai Fellowship",
"FEMALE": "340603",
"loc_id": "10672"
}


texts: Requires Datatype for inserting features

texts: I have imported postgis data in geoserver. In geoserver the spatial data column is displayed as MultiLineString, Does it affect when inserting features into geoserver or does it have to be explicitly defined as geometry?


texts: Cluster symbol color based on categorised variables

texts:
I'm making a categorized point map with QGIS. I'm using the categorized symbol and qualitative color palette. I need to group the points into clusters. Unfortunately I can't create clusters based on categories in a single layer. (the color assistant accept only numeric values).

However, on this post and this screenshot it seems possible...

https://north-road.com/qgis-cluster-renderer-crowdfunding/


texts: Land masking in ArcGIS Raster Calculator?

texts: I'm trying to mask land in a raster and just wondering if anyone knows how to do this. I tried using the raster calculator by producing a new raster based on a threshold value but my output raster is classified to 0 and 1 - 0 being my mask. I still want to retain the original values in my pixel.


texts: How to edit line shapefile with many segments to single segment?

texts:
As title, I have one line shapefile with many segments:

How to use QGIS to quickly edit this line shapefile and output a new file looks like the original one, but with only one segment (its attribute table only has one row)?


texts: Working with LiDAR data using other than Esri software?

texts: I will be working with LiDAR data soon and I am wondering what the alternatives to using Esri software to do this are, including open source solutions?


texts: Importing JSON or GeoJSON in Google Earth?

texts:
I was wondering if it was possible to import layers such as GeoJSON layers on the web-based version of Google Earth?

I didn't find any resources online but for the pro version and only for shapefiles.


texts: Make the name of BOOKMARKS unchangeable in QGIS

texts:
I am working in a large area and I need to use Bookmarks in QGIS.
When I click on them to zoom in in the interested area, accidentally I can change its name.

Is there any option to make Names as ReadOnly or to not have the opportunity to change them?


texts: Using gdal.Warp() on netCDF file. Simple case fails

texts:
I'm using the python API to GDAL to try to do a simple geographic warp of a netCDF4 satellite file to a small sub region.  My python program fails to produce any output but does not throw any exceptions or errors.  Any ideas?
The command line operation of gdalwarp works fine and produces meaningful netCDF output (for a single band or variable).  But the python version does not. THe output 'ds' is "None' and no netCDF output file is created. Here is a snapshot:
    from osgeo import gdal, gdalconst

kwargs = {'format': 'netCDF', 'copyMetadata': True,
          'outputBounds': '[-29 -2 -25 2]',
          'dstSRS':'+proj=longlat +datum=WGS84 +no_defs'
         }

src_filename = './subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc'
nc_file = 'NETCDF:' + src_filename
src = gdal.Open(nc_file, gdalconst.GA_ReadOnly)

subDatasets = src.GetSubDatasets()
#print('subsdatatsets', subDatasets)

ds = gdal.Warp('./output_lonlat.nc', src, **kwargs)
print("ds is ", ds)

del ds

Here is the command line operation that does work:
% gdalwarp -t_srs "+proj=longlat +datum=WGS84 +no_defs" -te -29 -2 -25 2 NETCDF:"./subset_20180729030000-JPL-L2P_GHRSST-SSTskin-MODIS_A-N-v02.0-fv01.0.nc":sea_surface_temperature -overwrite -of netCDF out_lonlat.nc


texts:
I have a legend created in ArcMap 10.1 which appears to be generating an abnormal white gap to the right of the legend content.

How can I remove this?


texts: Can ModelBuilder be used to export different versions of the same mxd with various constraints?

texts:
I am currently trying to use modelbuilder to go through each row in an attribute table, select the streets that intersect the polygon (which is the row in the attribute table), zoom to the selected features, apply a definition-query to show only the currently selected polygon and the roads, then export the current layout as a pdf.

For example, say you have a 4x4 group of polygons labeled like so

A B C D
E F G H
I J K L
M N O P

What I am trying to do is:

select the row that correlates to F
select the streets that intersect it,
zoom to the polygon with the streets selected,
apply a definition query so that I can only see F and the streets that intersect it, and then
export the map layout as a pdf and then
go on to do the exact same thing for G, then H, and so on.

Is this feasible to do with modelbuilder/python, or would it just be easier to use other methods like data driven pages and manually apply the definition-query?

I have included the portion of the code that I have managed to get working in model builder below.

import arcpy

arcpy.ImportToolbox("Model Functions")

cotgis11_COTGIS_Block_Watch__2_ = "cotgis11.COTGIS.Block_Watch"
Streets_59_999_to_10_000 = "Streets 59,999 to 10,000"
I_cotgis11_COTGIS_Block_Watc = "I_cotgis11.COTGIS_Block_Watc"

arcpy.IterateRowSelection_mb(cotgis11_COTGIS_Block_Watch__2_, "Block_Watc #", "false")

arcpy.SelectLayerByLocation_management(Streets_59_999_to_10_000, "INTERSECT", I_cotgis11_COTGIS_Block_Watc, "50 Feet", "NEW_SELECTION")

arcpy.CalculateValue_management("zoomToSelectedFeatures()", "import arcpy\\ndef zoomToSelectedFeatures():\\n        mxd = arcpy.mapping.MapDocument('CURRENT')\\n        df = arcpy.mapping.ListDataFrames(mxd,\"Layers\")[0]\\n        df.zoomToSelectedFeatures()\\n        arcpy.RefreshActiveView() ", "Feature Layer")


texts: Detailed Bathymetric/Bathymetry contour data (0-200M)

texts:
I am looking for free Bathymetric contour data in shapefile format for depths between 0-200 meters for North America. I prefer contours at 10 meter increments, but will take whatever I can get.

I have found Bathymetric layers for 0 and 200 meters from http://www.naturalearthdata.com/downloads/10m-physical-vectors/. However, I need more detail within that range.

My questions are:

Where can I find freely available Bathymetric data between 0-200 meters for North America?
How did you find the data and what did you use to search for it?


texts: Adding elevation data to a network dataset

texts:
I am trying to add elevation data to my road network to more accurately calculate walking distance to and from schools. I have figured out how to do the network with out elevation, but for the life of I can't figure out how to add in the elevation (hills etc).

I have currently have:

Road Network Feature Class
Road Infrastructure Feature Class
Contour Lines 1-5m Feature Class
TIN and DEM models produced from the
contour lines


texts: Get Feature that was snapped to in New Snapping Environment

texts:
I asked this over at the ESRI forums but lack of activity over there made me post here as well.

Does anyone know if there is a way to get the feature that something was just snapped to in the Snapping Environment in ArcGIS 10? I have code that will do this in the classic snapping mode but I would like to take advantage of the new snapping environment.

In classic snapping you could loop through the IFeatureSnapAgent.FeatureCache and perfrom a HitTest to see which feature was just snapped to.

I see you can get the location of the snapped point via ISnappingEnvironment.PointSnapper.Snap(), but this does not expose any information about what was snapped to. I see you can add/clear/update a geometry cache on the ISnappingEnvironment.PointSnapper but you can't get any of the items in cache. I'm not too sure that would be helped either way since its just geometries and not features. I guess in theory you could perform a spatial query against each layer in the map with the location of the ISnappingResult but I don't see this as something that is feasible performance wise.

I am missing something or has anyone else dealt with this limitation of the new Snapping Environment?

Matt


texts: Reading label positions from attribute fields with Python in QGIS 3?

texts:
I'm trying to label a shapefile with coordinates from attribute table with pyqgis.

I've been messing around with labeling shapefiles with python, I can do most things that I need, but the main setting still eludes me. How do I set data defined coordinates with python from attribute table (label_x, label_y).

Below is just a sample code for triggering labels but it does not have anything in it to trigger labels from coordinates.

myLayer = QgsProject.instance().mapLayersByName('test_poligon')[0]
settings  = QgsPalLayerSettings()
settings.fieldName = "Fid"
settings.placement = 3
settings.enabled = True
settings = QgsVectorLayerSimpleLabeling(settings)
myLayer.setLabelsEnabled(True)
myLayer.setLabeling(settings)
myLayer.triggerRepaint()


texts: Histogram charts rendering away from data

texts:
I am having issues with producing Histograms in my file. Initially I thought they were not showing up, but it turns out they are, just extremely far away from my actual data. For reference:

Here are the attributes I have set up. Is there anything in here that could be causing this issue? Is it possible that there is something else going on that is causing this to happen?


texts: How can I generate polygons like in a Voronoi diagram/Thiessen Polygons and get the nodes for each of the polygons?

texts:
I want to be able to generate a set of polygons where I can extract all nodes for every polygon. For example:

Polygon 1 - (0, 0),(0,2),(2,0)
Polygon 2 - (0, 2),(2,2),(2,0)
Polygon 3 - (0, 2),(5,5),(8,5),(8,0),(2,0)
And so on...

I'm not interested in available tools since this is just a part of a project. I want to be able to insert random points to generate this random dataset with polygons and its nodes coordinates.
Where do I start? Is there an algorithm I can implement in a programming language?
BTW: the data is supposed to be used in a PostgreSQL database as geometries.


texts:
Can the grid labels be written along the ends of the grid lines within the data box without Cartography tools and with E / N at the end?

I've been playing with offsets and other options. The only result I got were labels within the box but not along the lines. The chart is not north up.

Editing the labels is not what I'm looking for as there is quite a lot of these charts.


texts: GeoServer: serve series of GeoTIFF files in YYYY/MM/DD/ structure

texts:
I have a growing archive of satellite data in GeoTIFF format that is stored on disk in a YYYY/MM/DD/ directory format. Each day, the most recent data is added to this archive.
I've spent the last few weeks familiarizing myself with GeoServer and know how to add a single file to a store, and visualize the data through a WMS connection. But know I am ready for the next step and point GeoServer to my data archive.
How could I serve such an archive with GeoServer and set up a WMS service for it?


texts:
Is there a way to enable map tips for a DEM raster layer? I want to hover over a raster and get the values.

I usually do this in ArcMap and I am trying to learn if it is possible in QGIS.


texts:
I am writing in the hopes that there is somebody here that might know/use a GIS server solution, that has WPS functionality which would allow me to use data from a WFS 2.0.0 GIS service, which serves GML 3.2 data.

An example of service I am trying to use is this INSPIRE download service:

http://inspire.biodiversity.ro/WFS/RO_ENV_PS/wfs?service=wfs&version=2.0.0&request=GetFeature&TypeName=ps:ProtectedSite&featureid=RONPA0056

From what I read, WPS services can't process data from a WFS directly, and first needs to read the geometries in memory, to be then able to geoprocess.

I have tried several tools from the Geoserver 2.7, using the demo WPS interface to create processing requests however there seems to be no support for WFS 2.0 and GML 3.2 geometries. I have tried using the process vec:CollectGeometries but it only supports WFS 1.1 or WFS 1.0 and not the WFS 2.0 standard.

Does anyone have a way of working with the service above?


texts:
I'm using QGIS 3.16.1.  I'm creating a heatmap using the symbology.  At first it appears like a normal heatmap, but it is still rendering.  When it finishes the whole heatmap turns red.  Here's a before and after screenshot.
What is happening here?


texts:
If tableField = "Process ID" and layerField = "PID" then this works:

IRelationshipClass relationshipClass =
    memoryRelationshipFactory.Open(
    "Join",
    table as IObjectClass,
    tableField,
    geoFeatureLayer.DisplayFeatureClass as IObjectClass,
    layerField,
    "forward",
    "backward",
    esriRelCardinality.esriRelCardinalityOneToOne );

If tableField = "Process ID, Year Made" and layerField = "PID, YM" then it fails.

(Table is SQL Server table.)

I have tried brackets and quotes to no effect. Is it possible to pass a compound key?


texts: Adjust raster value if cell contains a point in ArcMap?

texts:
My raster consists of 100x100m cells that all have the value 0, and a shape with point data. I'd like to update all the cells that contain a point to 1, but can't seem to find a tool that does this. It is possible to do this using a tool (or multiple tools) from the toolbox?

What I mostly find when trying to look this up is 'point to raster', which I don't think is what I need, and I believe 'con' does not work with shapefiles.

I'd like to do it in ArcMap (10.4). I'm trying to avoid using a python script, as I'm having lots of trouble loading the arcpy module on my computer.


texts: How would I group values in a field based on range of number using python, and then update the attribute table?

texts:
An example of what I'm trying to accomplish is described below, however they don't describe how to use the data update cursor to update the attribute column.

data - Group values based on range of number in python - Stack Overflow
https://stackoverflow.com/questions/20437394/group-values-based-on-range-of-number-in-python

fc = r'C:\Users\Test.gdb\Test1a'
print "fc set"
class_field = "Numbers"
name_field = "GroupedID2"
lst = sorted(arcpy.da.SearchCursor(fc, [class_field,name_field]))
class Delta:
    def __init__(self, delta):
        self.last = None
        self.delta = delta
        self.key = 1
    def __call__(self, value):
        if self.last is not None and abs(self.last - value[0]) > self.delta:
            # Compare with the last value (`self.last`)
            # If difference is larger than 2.5, advance to next project
            self.key += 1
        self.last = value[0]  # Remeber the last value.
        return self.key

import itertools
for key, grp in itertools.groupby(lst, key=Delta(2.5)):
    for tup in grp:
        print(tup + ('project{}'.format(key),))
with arcpy.da.UpdateCursor(fc, "GroupedID2") as cursor:
    for row in cursor:
        cursor.updateRow(row)
        print (row[0])

The trouble that I have is that no matter where I place the update cursor, I can't get it to record each row correctly.  It will make the entire attribute the last project number or the first project number instead of rows 1-3 project1, 4-6 project2, and 7-8 project3.

This is what the code prints.  The feature class remains project 1 but it isn't updated to include the correct project number from the 3rd column.

2014-07-16 07:26:59.135000
fc set
(1279015.0, u'project1', 'project1')
(1279017.5, u'project1', 'project1')
(1279020.0, u'project1', 'project1')
(1279022.5, u'project1', 'project1')
(1279027.5, u'project1', 'project2')
(1279030.0, u'project1', 'project2')
(1279032.5, u'project1', 'project2')
(1279037.5, u'project1', 'project3')
(1279040.0, u'project1', 'project3')
(1279042.5, u'project1', 'project3')

when I alter the code to read this:

import itertools
for key, grp in itertools.groupby(lst, key=Delta(2.5)):
    for tup in grp:
        print(tup + ('project{}'.format(key),))
        test = ('project{}'.format(key))
        with arcpy.da.UpdateCursor(fc, "GroupedID2") as cursor:
            for row in cursor:
                test = row[0]
                cursor.updateRow(row)
                print (row[0])
print(datetime.now()-startTime)

It still keeps it as project 1.


texts: Batch clipping all files in folder and then save using ArcPy

texts:
I'm trying to use a script to batch clip all files in a folder, then saved with a _clipped.shp at the end of the original, I need a script that is something like this:
arcpy.env.workspace = "E:\\folder\\test_input"

# Set local variables
in_features = arcpy.env.workspace
clip_features = "E:\\folder\\mask\\boundary.shp"
out_feature_class = "_clipped.shp"

# Execute Clip
for f in arcpy.env.workspace
    arcpy.Clip_analysis(f, clip_features, f(name) + out_feature_class)

Where the f(name) 'points' to the orignal shapefile name and then 'appends' _clipped.shp at the end and then saves and performs this process iteratively for each file.
Also should I have a different output folder for the clipped features?
The reason I'm doing this is because the actual batch clip feature of the clip tool is not working for some reason.