texts: Why would FME imported data into SQL Server not being found by QGIS?

texts:
Why is FME imported data into SQL Server not being found by QGIS?

I've imported some MapInfo format data into SQL Server 2014 using FME.

In SQL Management Studio (SSMS), I can do a query and the spatial results tab shows the various polygons as expected

But if I switch to QGIS, create a connection to the database, and then open the list of tables in the dbo schema, the data set isn't showing

I have also loaded data using ogr2ogr and those data sets are showing, and can be displayed on a map in both QGIS and SSMS

So why aren't FME loaded data sets showing?

Additional Info

I deleted the old MSSQL connection, and when I created a new one, I noticed the the "Only look in the geometry_columns metadata table" and
"Also list tables with no geometry" options. Unchecked first and checked the second, and now all the tables are showing, both the ones loaded using ogr2ogr and using FME.

Now when I try and add a FME layer, I get the error ...

"Layer is not valid: The layer dbname='GisData' host=(local) user='localdev' password='password' estimatedmetadata=true srid=27700 type=LineString table="dbo"."FmeFootpaths" (GEOM) sql= is not a valid layer and can not be added to the map"

The thing that caught my eye was the srid=27700 ... the ogr2ogr has srid 32768  (Should point out I don't understand srid's?)

Also I checked this geometry_columns, and it only has the tables loaded using ogr2ogr


texts: Finding parcel polygon features with frontage on multiple streets using ArcGIS Desktop?

texts:
The Department of Public Works would like to know how many properties in the city have a sidewalk at both the front and back.  To get at this, I'm looking for a method to find parcels with frontage on more than one road.  I'm working in ArcInfo 9.3.1 (file geodatabase).  I'm looking for a solution that doesn't require scripting or downloading other tools. Street centerlines and parcel polygons are available; there is no sidewalk layer.

Using some of the ideas from this forum, I have developed and tested a process.  Some of the results are puzzling (corner lots may or may not be included in the final results-I think DPW doesn't want corner lots included) and I'm concerned about scaling it up to the full 42000 parcels (I tested 1223 parcels).
Are there steps to add to this method to improve the results?  Is there another way to accomplish the same thing that could be more automated through model builder?

Run Feature to Line tool on parcel polygons feature class (two lines are created for boundaries shared by polygons and lines retain parcel id)
Build a geodatabase topology with this newly created parcel line feature class using the rule “must not have overlaps” then validate topology (interior lines will be in error)
In ArcMap, add topology and feature class and open error inspector then search for errors > select all errors > right click > select features > delete (to delete interior lines)
Run Frequency on remaining features in parcel line feature class with parcel id as frequency field – select records in table where Frequency > 1 and export to new table
Join Frequency>1 table to parcel polygons based on parcel id, keeping only matching records – export join results


texts: Reading value and coordinates of every pixel post NDVI calculation in ArcGIS Desktop?

texts:
I have a pre forest fire image and a post forest fire image. I would like to make a comparison between the coordinates and the values of both images pixel by pixel.

How would I access the value and coordinate of each pixel, post NDVI calculation in ArcGIS?


texts: Sentinel-1 SAR image is displayed upside down

texts:
I am looking at Sentinel-1 SAR image downloaded from ESA Copernicus - dhus location around south of Tokyo, Japan. File format is GeoTIFF.

It looks good in QGIS.

However, if I use SAGA,

or SNAP, or OTB-Monteverdi,

the image looks upside down. Maybe "upside-down" is not an appropriate word, as the image lost its ENE-WSW orientation. The image seems simply drawn in XY space swapping lower-left and upper-left as origin.

In addition to QGIS, ArcGIS shows it correctly (according to my friend).

Sentinel-1 images are huge (several GBs) so I cannot test many images, but I reproduced this with randomly selected image in Australia.

Perhaps there are some geotags missed by these software?

Is there a way to modify the image (geotags, to be more specific) to properly read and display it in packages other than QGIS and ArcGIS software?


texts: Openlayers2 how to display WFS with projection other than EPSG:4326

texts:
Here is some OpenLayers code:

var projection = "EPSG:3009";

map = new OpenLayers.Map('map', {

//    projection: "EPSG:900913",
//    displayProjection: projection,
//    displayProjection: "EPSG:4326",
});

var osm = new OpenLayers.Layer.OSM();

var protocol = new OpenLayers.Protocol.WFS({
    version: "1.1.0",
    srsName: projection,
    url: "http://" + window.location.hostname + "/geoserver/test/wfs",
    outputFormat: "application/json",
    readFormat: new OpenLayers.Format.GeoJSON(),
    featurePrefix: "test",
    featureNS :  "test",
    featureType: "eksjo_trad",
    maxFeatures: 1000,
    geometryName: 'geom',
});

wfs = new OpenLayers.Layer.Vector("WFS", {
    projection: new OpenLayers.Projection( projection ),
    strategies: [
        new OpenLayers.Strategy.BBOX(),
    ],
    protocol: protocol,
});

map.addLayers([ osm, wfs ]);
map.zoomToMaxExtent();

My problem is that OpenLayers draws the points in France instead of Sweden. I am including the proj4.js library, although this appears to make no difference.

If I redefine the 'projection' variable to be "EPSG:4326" the points appear in the right place. However, a user needs to be able to edit the co-ordinates in EPSG:3009 - I get the right co-ords when 'projection' = "EPSG:3009".

Here is the JSON with the EPSG:3009 co-ords for 1 feature returned by the WEF:

{"type":"Feature","id":"eksjo_trad.1","geometry":{"type":"Point","coordinates":[146224.14140031,6392371.57291367]},"geometry_name":"geom","properties":{"trad_id":1,"layer":"VG_TRAD_BARRTRAD__DIGITALISERAD 11000","subclasses":"AcDbEntity:AcDbPoint","extendedentity":null,"linetype":null,"entityhandle":"146B","text":null}},

Can anyone help me?


texts: Automatic numbering of points during digitization in QGIS?

texts: Does somebody know how to automatically number the points during digitization in Qgis. To explain, when I create the new point there is the window where I must to type the point number. Is there any way to automatic number without manually typing of point number?


texts:
I am new in GIS.

I am trying to use the curves of level (downloaded from this site in SHP format) to use them into Tableau Desktop.

The problem is that if I loaded the lines, they are not correctly setted, as shown in this picture from the same file loaded with QGIS.

How could I fix it?

In the archive I downloaded from the site linked above there are also other 3 files, with extension dbf, prj and shx that maybe can be useful to center the dataset.


texts: NHGIS census files get lat/lon from polygon

texts:
I've never worked with shapefiles before and I'm having a lot of trouble with something that I think should be pretty simple.

So I have the census tract files from NHGIS for 1980 and the structure looks like:

NHGISST     NHGISCTY    GISJOIN     GISJOIN2    SHAPE_AREA  SHAPE_LEN   geometry
0   560     0210    G56002100016    56002100016     2.569850e+09    241261.569958   POLYGON ((-732462.5590209048 499607.0109478103...
1   560     0210    G56002100013    56002100013     9.690180e+06    14271.148305    POLYGON ((-730612.5504814894 448835.8542251377...
2   560     0210    G56002100014    56002100014     1.686311e+07    19833.251952    POLYGON ((-728941.3490845333 443959.3627771207...
3   560     0210    G56002100008    56002100008     4.748756e+06    10873.293718    POLYGON ((-733509.072211165 444652.6344217169,...
4   560     0210    G56002100009    56002100009     3.887032e+06    10440.687964    POLYGON ((-734860.1103011537 445048.4184498419...

I am trying to use the geometry polygon structure to extract a lat/lon centroid but nothing seems to work.  Whenever I use something like MultiPoint or Proj (I'm using python by the way) I get coordinates that are completely off.

I'm completely stumped.

def to_lat_long(t):
    inProj = Proj(init='epsg:3857')
    outProj = Proj(init='epsg:4326')
    return transform(inProj,outProj,t[0],t[1])
break_count, break_count_max = 0, 1
for k,v in shape_file_store_1.items():
    if type(v[1990]) == str:
        x = v[1990][1:]
        for s in shp_1990:
            if s == 0:
                first = next(iter(shp_1990))
            else:
                first = s

            shp_file = first['properties']
            items = list(shp_file.items())
            s_file = items[3][1]

            coords = first['geometry']['coordinates'][0]

            points = MultiPoint(coords)
            points_c = points.centroid
            points_r = points.representative_point()

            print(f'Centroid: {points_c}')

            centroid_latlong = to_lat_long(points_c.coords[0])
            df_townhalls.at[df_townhalls['geo'] == v['geo'], '1990_lat'] = centroid_latlong[0]
            df_townhalls.at[df_townhalls['geo'] == v['geo'], '1990_long'] = centroid_latlong[1]

            print(v['geo'])
            #print(f'Representative Point: {points_r}')
            #print(df_townhalls.loc[df_townhalls['geo'] == v['geo']])

            break
            if s_file == x:
                v[1990] = first
        #break
        break_count += 1
        if break_count > break_count_max: break
    #break


texts: Different performance of counting features in ogr.Layer returned from ExecuteSQL

texts:
I am working on a PyQt5 desktop application GDBee which allows executing SQL queries against Esri file geodatabases using GDAL. My development environment is Windows 10 64bit, Python 3.5, and Python bindings for GDAL 2.1 / 2.2. I use the built-in OpenFileGDB driver (read-only access).

This is the minimal code you need to reproduce the problem. The data are coming from the PostGIS workshop materials by Boundless (shapefiles have been loaded into a file geodatabase as is).

import time
import ogr
ds = ogr.Open(r"NYC.gdb")

res = ds.ExecuteSQL('select * from streets', dialect='SQLITE')

start_time = time.time()
print(('Number of feats using len', len(res)))
end_time = time.time()
print(('Time with len', end_time - start_time))

start_time = time.time()
print(('Number of feats using GetFeatureCount', res.GetFeatureCount()))
end_time = time.time()
print(('Time with GetFeatureCount', end_time - start_time))

res = ds.ExecuteSQL('''SELECT
    Polys.Name AS NeighborhoodName
    , Count(*) AS CrimeCount
FROM
  Homicides AS Points
JOIN
  Neighborhoods AS Polys
ON
    ST_Contains(Polys.Shape, Points.Shape) = 1
AND
    Points.YEAR in (2008, 2009, 2010)
GROUP BY
    Polys.Name
ORDER BY
    CrimeCount DESC
LIMIT 10
''', dialect='SQLITE')

start_time = time.time()
print(('Number of feats using len', len(res)))
end_time = time.time()
print(('Time with len', end_time - start_time))

start_time = time.time()
print(('Number of feats using GetFeatureCount', res.GetFeatureCount()))
end_time = time.time()
print(('Time with GetFeatureCount', end_time - start_time))

The output:

('Number of feats using len', 19091)
('Time with len', 0.022000551223754883)
('Number of feats using GetFeatureCount', 19091)
('Time with GetFeatureCount', 0.021004199981689453)
('Number of feats using len', 10)
('Time with len', 5.007046937942505)
('Number of feats using GetFeatureCount', 10)
('Time with GetFeatureCount', 5.008031845092773)

I have noticed that the performance of getting the number of features returned by a complex SQL query takes a lot longer (5 secs) even though the number of features is tiny (just 10 features). It can have something to do with the spatial functions, but I am not completely sure. The performance of the len and the GetFeatureCount is practically identical.

Does anyone know why getting the length of features in the Layer object returned is different depending on the SQL query executed and if there is any way to make it faster?


texts: Converting HDF to TIFF using batch processing with gdal_translate

texts:
I have a huge MODIS dataset and I need to translate them to tiff from hdf.
My input path is: C:\Users\HP\Desktop\deneme\ham_veri
and hopefully, the output path is: C:\Users\HP\Desktop\deneme\yeni
This is what I'm trying to run for batch processing:
FORFILES /m MOD10A2.A2016321.h20v05.006.2016340133152.hdf C:\Users\HP\Desktop\deneme\ham_veri "cmd /c gdal_translate -sds -of HDF4_SDS:hdf:@file:01 C:\Users\HP\Desktop\deneme\yeni\MOD10A2.A2016321.h20v05.006.2016340133152.tif"

and this is what I got from command prompt:
ERROR 6: Too many command options '%C:\Users\HP\Desktop\deneme\yeni%\@MOD10A2.A2016321.h20v05.006.2016340133152.tiff'

Where am I doing wrong?


texts: Monthly average of multiple Tropomi CH4 files using ArcGIS Pro

texts:
I have 30 images for CH4 emissions for one month but in float32 type. Some images have no data at all and some images have few pixels with data such as this.
When clicking the black pixels it shows NaN values.
I used cell statistics and raster calculator but the output is not what I want. I tried three ways:
1- I downloaded image without mentioning anything of how the pixels which have no data should be shown in ArcGIS Pro, when opened in ArcGIS Pro those pixels show nan values, I then tried averaging the 30 images but the result of that is no image and value as shown below.
2- Then I downloaded images with -9999 and Nan values given to pixel with no values and the result of averaging for NaN is this  and for -9999 it is different, it shows the image plus the values but in negative and addition of all the values rather than their mean.
Is there a possible solution or an error caused by me?


texts: Categorize difference between max and min value for features with same attribute

texts:
I have a points layer. Some points have the same value in FIELD_1. I want to give different colors (value in FIELD_3) for different sizes of ranges of values in field FIELD_2 among points with same FIELD_1.

I'll try to be clearer with an example:

If points have same FIELD_1 and FIELD_2 ranging from 3 (min value) to 10 (max value) return (in FIELD_3) 10-3 = 7

If points have same FIELD_1 and FIELD_2 ranging from 2 (min value) to 24 (max value) return (in FIELD_3) 24-2 = 22

etc...


texts:
I have two polygon layers, Land-use and Mesh.

I had this mesh layer in specific geometric way. Anyways, I want the mesh layer to have the landuse value of the maximum area it contains. Lets say that a polygon in the mesh contains three land use types (40% urban, 30% roads, 30% green areas). I want it to assign the value of urban because it has the biggest area.

I tried Spatial Join tool, but it did not give me what I was looking for. I used the (mode) as a Field Map of Join Features which should Use the value with the highest frequency but it did not work.


texts: Defining indices for moving/sliding window rasterbrick in r

texts:
I have a rasterbrick of bioclimatic variables for a period, e.g.
 rasterBrick
class      : RasterBrick
dimensions : 3875, 8017, 31065875, 10  (nrow, ncol, ncell, nlayers)
resolution : 0.04490319, 0.04490319  (x, y)
extent     : -180.0393, 179.9495, -90.04088, 83.95898  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : /tmp/RtmpdAQr8j/raster/r_tmp_2020-10-29_083822_35311_30585.grd
names      : biovars2010x, biovars2011x, biovars2012x, biovars2013x, biovars2014x, biovars2015x, biovars2016x, biovars2017x, biovars2018x, biovars2019x
min values :     221.4524,     221.9133,     220.9857,     221.2995,     220.7355,     221.3381,     219.8078,     221.7039,     220.2983,     220.9756
max values :     306.6502,     306.7648,     306.5869,     307.1276,     307.3617,     307.3223,     307.1781,     307.0450,     307.1749,     307.3326

I need to do a sliding window where I calculate the mean over a period, e.g., 3 years, and then for the next iteration, start the calculation from the second year in the previous period.
E.g.: 2010-2013; 2011-2014; 2012-2015 and so forth.
So I set some z-values for the brick as follows:
s = setZ(rasterBrick, 1:nlayers(rasterBrick))

 s
class      : RasterBrick
dimensions : 3875, 8017, 31065875, 10  (nrow, ncol, ncell, nlayers)
resolution : 0.04490319, 0.04490319  (x, y)
extent     : -180.0393, 179.9495, -90.04088, 83.95898  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
source     : /tmp/RtmpdAQr8j/raster/r_tmp_2020-10-29_083822_35311_30585.grd
names      : biovars2010x, biovars2011x, biovars2012x, biovars2013x, biovars2014x, biovars2015x, biovars2016x, biovars2017x, biovars2018x, biovars2019x
min values :     221.4524,     221.9133,     220.9857,     221.2995,     220.7355,     221.3381,     219.8078,     221.7039,     220.2983,     220.9756
max values :     306.6502,     306.7648,     306.5869,     307.1276,     307.3617,     307.3223,     307.1781,     307.0450,     307.1749,     307.3326
time       : 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

And then I made a vector of the moving window indices as follows:
v <- 0:7
v<- c(sapply(v, `+`, 1:3))
 v
 [1]  1  2  3  2  3  4  3  4  5  4  5  6  5  6  7  6  7  8  7  8  9  8  9 10

And then did the following:
biovar.mean = zApply(s, by=v,fun=mean)

But I get the following error:
Error in setZ(b, zval, name) : length(z) == nlayers(x) is not TRUE
In addition: Warning message:
In ind[] <- indices :
  number of items to replace is not a multiple of replacement length

...but that obviously doesn't work. Any ideas how I can better specify the sliding window?


texts: pyqgis change or update attribute data with real-time data or Timer update，Memory has been growing

texts:
I update my shapefile attribute data by QTimer ,and labling it,
keep the lable visiable,change attribute by 1 minute,
but the Memory has been growing

1

    layers =self.iface.legendInterface().layers()
    selectedLayerIndex=self.dlg.labelfield_2.currentIndex()
    sl=layers[selectedLayerIndex]
    CurrentSecond=QTime.currentTime().second()
    self.ldp = sl.dataProvider()
    features = self.ldp.getFeatures()
    fdindex=self.ldp.fieldNameIndex(self.dlg.labelfield.currentText())

    fdindex=self.ldp.fieldNameIndex(self.dlg.labelfield.currentText())
    for feature in features:
        self.UPDATEMAP[feature.id()] = { fdindex: str(CurrentSecond) }
    self.ldp.changeAttributeValues( self.UPDATEMAP )
    self.iface.mapCanvas().refresh()

2

    layers =self.iface.legendInterface().layers()
    selectedLayerIndex=self.dlg.labelfield_2.currentIndex()
    sl=layers[selectedLayerIndex]
    CurrentSecond=QTime.currentTime().minute()
    self.ldp = sl.dataProvider()
    features = self.ldp.getFeatures()
    fdindex=self.ldp.fieldNameIndex(self.dlg.labelfield.currentText())
    sl.startEditing()
    for feature in features:
        feature[fdindex]=str(CurrentSecond)
        sl.updateFeature(feature)
    sl.commitChanges()
    self.iface.actionSaveEdits()
    if iface.mapCanvas().isCachingEnabled():
        sl.setCacheImage(None)
    self.iface.mapCanvas().refresh()

boththis two methods  have the metmory growth,so what can I do?


texts: Which topics are essential for QGIS beginners?

texts:
What topics do you think are essential in an introduction for QGIS?

The target audience is those with zero or rudimentary knowledge in mapping and geospatial software.


texts:
I have a data set with 1000 different latitudes-longitudes. I wish to extract average annual temperature and annual precipitation for each of these coordinates. These data can easily be obtained from WorldClim and processed using DIVA-GIS.

Is there anyway to do this on R?

I want my final output to be a dataframe with the annual temperature and precipitation for each coordinate. I am new at GIS in R, so I seek a basic code chunk along with the required libraries for this output.


texts: Performing cs2cs datum transformation using GDAL?

texts:
I'm using GDAL 3 to perform datum transformations. I'm trying to transform coordinates from a local datum to WGS84 using the following command:

cs2cs +proj=latlong +ellps=clrk80 +towgs84=-166,-15,204 +to +proj=latlong +datum=WGS84

The datum has parameters: deltax=-166, deltay=-15, deltaz=204 and the ellipsoid is Clark 1880.

When I put the coordinates "5, -15" in then I receive the transformed values of:

5d0'0.984"W 15d0'2.12"N 0.00

This point comes from test data, and the actual transformed coordinate should be - 4° 59' 51.88272", 15° 0' 0.92437", - 71° 52' 38.87483". Is there something wrong with my GDAL installation, or did I input the command incorrectly?


texts: Is there a way to fill a column in the attribute table with the same value, similar to a paste .  Is there something in the calculator that can be used to fill selected records/attributes.


texts:
I would like a user to select from a form a particular thing to measure, and for the form to close and show the Draw toolbar as the toolbar currently focused. My current code is (when a button is clicked):

function measureFeature(layer_name) {
    /* Activates the layer for editing and closes the form */

    closeForm(); // function defined elsewhere to close form
    Application.Map.Layers(layer_name).Editable = true;
    // Focus 'Draw' toolbar here
}

All I am missing is something along the lines of:

Application.Toolbars.Item('Draw').Focus();

or

Application.focusedToolbar = 'Draw';

Does this even exist, maybe under a different method or object?

Also asked here: https://geonet.esri.com/thread/55469, but with no answer in 2 years


texts: Using wild cards in GRASS to remove files

texts:
I have a bunch of files on a GRASS location, all named temp_something, that I want to delete. I tried to use the g.remove tool using the multiple option and typing temp* in, but it doesn't work out, and a character * not allowed error is reported.

Am I doing something wrong or is there just no way GRASS works with wildcards?


texts: WFS not in position when using a Google Base Layer

texts:
I use two base layers and user can switch from OSM MapQuest to Google from a drop down list. However my WFS layer is not in position (shift of about 5 degrees to the East) when the base layers is a Google Base Layer. Position is fine when the OSM layer is used.

I have read a few posts about related problems which were caused by wrong coordinate systems for the map (not using spherical Mercator for the map projection) but I tried to set it properly.

    geographic = new OpenLayers.Projection("EPSG:4326");
    mercator = new OpenLayers.Projection("EPSG:3857");

    OLmap = new OpenLayers.Map('map', {

        projection: mercator,
        //maxExtent: extent,

        controls: [
            new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.MousePosition(), //{div:document.getElementById("coordinates")}
            new OpenLayers.Control.ScaleLine(),
            new OpenLayers.Control.Zoom()
        ],
        zoomDuration: 2,
        numZoomLevels: 10,
        displayProjection: geographic
    });

    arrayOSM = ["http://otile1.mqcdn.com/tiles/1.0.0/map/${z}/${x}/${y}.jpg",
                "http://otile2.mqcdn.com/tiles/1.0.0/map/${z}/${x}/${y}.jpg",
                "http://otile3.mqcdn.com/tiles/1.0.0/map/${z}/${x}/${y}.jpg",
                "http://otile4.mqcdn.com/tiles/1.0.0/map/${z}/${x}/${y}.jpg"];

    baseOSM = new OpenLayers.Layer.OSM("Political map", arrayOSM,OpenLayers.Layer.OSM,
        {'maxExtent': new OpenLayers.Bounds(-20037508.34,-20037508.34,20037508.34,20037508.34)});
    OLmap.addLayer(baseOSM);

    var gphy = new OpenLayers.Layer.Google(
        "Google Physical",
        {type: google.maps.MapTypeId.TERRAIN}
        // used to be {type: G_PHYSICAL_MAP}
    );

    OLmap.addLayer(gphy);

    Feature_centroid_wfs = new OpenLayers.Layer.Vector("site",
    {
        projection: geographic // layers' original projection
        ,styleMap: StyleObject
        ,rendererOptions: {zIndexing: true}
        ,strategies: [new OpenLayers.Strategy.Fixed()]
        ,protocol: new OpenLayers.Protocol.WFS({
            version: "1.1.0"
            ,url: "http://myserver/geoserver/wfs"
            ,featurePrefix: "ramsar_sdi" //geoserver worspace name
            ,featureType: "features_centroid_published" //geoserver Layer Name
            ,featureNS: "http://myserver/geoserver/ramsar_sdi" // Workspace URI
            ,geometryName: "geom" // field in Feature Type details with type "Geometry"
            ,srsName : "EPSG:4326"
            ,outputFormat: "application/json"
            ,readFormat: new OpenLayers.Format.GeoJSON()
            //,maxExtent: extent  // necessary if BBOX()
        })
    });

    OLmap.addLayer(Feature_centroid_wfs);


texts: Pgrouting: Visualizing analyzegraph ring geometries

texts:
I have a sample network that I analyzed using pgr_analzyeGraph. The outputs shows that I have 3 ring geometries. How do i get which nodes is part of the ring geometries?

The official document only shares how to see dead ends and where potential gap problems exist. Need to pinpoint where the ring geometries exist at in the tables.


texts:
I am trying to flip (mirror) a raster upside-down. Using following command incoporating numpy function flipud()

gdal_calc.py -A input.tif --calc="flipud(A)" --output output.flipud.tif

resulted with output.flipud.tif remains the same orientation as input. However, with

gdal_calc.py -A input.tif --calc="fliplr(A)" --output output.fliplr.tif

results with output.fliplr.tif which is flipped left-right from the input.tif

Does anyone know why flipud() does not work in this case?


texts:
I'm working with Sentinel 1 radar imagery. For some images in the collection in parentheses next to the angle band the is a percentage (see image)

When I do some further steps in my analysis these images return a masked pixel and I therefore need to select a different image. What does this percentage refer too? And is there a way to filter images where this occurs so I can remove them from the ImageCollection I couldn't find an obvious image property?


texts: In QGIS I'm attempting to create a buffer zone around a number of connected lines, then use the buffer layer to clip another layer. The lines features were provided as a shape file, and buffering works fine. However when I try to clip with the buffered layer I get an error stating 'invalid geometry'. Upon closer inspection, there is a line not connected to the others (see image), which appears to be the source of the error. How do I connect this line to the other? I've experimented with edit features (extend, reshape etc.) but can't seem to find the function I require to link them.


texts: Using Magellan Meridian and QGIS GPS Tools

texts:
I am using QGIS 2.0.1 on Ubuntu 12.04.  I have a Magellan Meridian GPS and I trying to use the GPS Tools plugin to download waypoints to a layer in QGIS.  I have successfully used the Gebabbel interface with GPSBabel to extract waypoints from the GPS to a .gpx file.

The only way I have successfully used the GPS Tools to download waypoints and add to a layer is by using the following command in the GPS device editor.  This is the same command that Gebabbel issues.

gpsbabel -w -i magellan,baud=9600 -o gpx -f /dev/ttyS0 -F /home/robert/GPSpoints/waypoints.gpx

I also use “waypoints” for the Layer name and Output file text boxes on the main GPS Tools dialog window.

This method seems to only work once.  If I remove the layer from QGIS, delete the .gpx file, and repeat the download request in the GPS Tools, I get an error message or QGIS becomes unresponsive,

I know the Meridian is an old unit, but is anyone successfully using it with QGIS GPS Tools.


texts:
----------x-----------------x----------------x---------------
--------------x-------------------x--------------------------

This is a Two-way road, each line stands for one way, and the x means a split point in the way. That's to say, we will split each way to multiple segments before we save them to postgis. Each segment have a name same as the road name.

For example, the name of the above road is St1, so the 7 segments will all have the same name.

Now for a given name like St1, we want to combine all the matched segments as one (if one way) or two(if two way) linestrings.

After googling, it seems that the pg function ST_Union (https://postgis.net/docs/ST_Union.html) looks like what we want, after the test, it worked.

However there are some different roads(in different position) may own the same name.

As shown, suppose the red road in the top right have the same road name as "West Chester Pike", but they are two different roads in deed.

So if we may extract all the roads by a given name "West Chester Pike", ST_Union will merge all of them, so we wonder if it is possible to distinguish the roads for example by distance?


texts: Summarizing column in QGIS field calculator?

texts:
I have data with a column containing a number value, and I want to divide the value of each feature by the total sum of all the values in the column (to get a percentage).

So I would like to know if there is some simple syntax in the QGIS field calculator to get the sum of a column, preferably without resorting to for loops and what not.

I have a shapefile where feature A has the value X in column B. How can I summerize the values for all features in column B in order to divide value X to get a percentage for each feature?


texts: Understanding KML To Layer performance in ArcGIS Desktop?

texts:
I've got a very simple code using an ArcPy tool to convert KMLs to Layers in a directory:

import arcpy, os
from arcpy import env

arcpy.env.workspace = r"C:\KML\Test"
outLocation = r"C:\KML\Test\Convert"

for kml in arcpy.ListFiles('*.KML'):
    print "Converting {0}".format(os.path.join(arcpy.env.workspace, kml))
    arcpy.KMLToLayer_conversion(kml, outLocation)

Each KML is no less than 20,000 KB and probably has 2000-3000 records per. On average, it is taking about 1.4 - 2 minutes for each KML to convert using the arcpy.KMLToLayer_conversion().

Any thoughts on if this is a memory issue? Not sure why it appears to be running so slow.

Also of note, import arcpy causes the script to delay for about 5 seconds when it starts.

This code is run outside of ArcMap through IDLE.


texts:
I'm looking for a shapefile for the coastline of the UK.

I'm guessing the Ordnance Survey 'geo-portal' is the best place to get this, but whenever I go to the Ordnance Survey site I get the shakes. (Seriously, trying to find any data on that site is really traumatic. I don't have a spare three days.)

Does anyone just have a simple URL for such a shapefile?


texts: How to create an accurate 5 km buffer in ArcGIS 10.1?

texts: I need 5 km buffer of  a basin boundary of Brazil using ArcGIS 10.1. The basin boundary is in GCS. I have tried it giving buffer distance 5 km and 0.0454545 degree decimal (assuming 1degree=110km). I have got two different results for two distances. Which one is more accurate?


texts: Sites that demo freely available Basemaps?

texts:
I've seen such sites in the past.
I remember seeing one put together by Wm Lerer at a conference a couple years ago.

What I'm looking for is a map that demonstrates all (well, at least a lot) of the freely available basemaps (tilesets, web services, etc) that I could include in a mapping application (using OpenLayers, Leaflet, or whatever).


texts: Calculate road distances between points in QGIS 3.4?

texts:
I would like to calculate road distances between a series of points. I don't have any road network layer, I need to rely on OpenStreetMap.

What is the best plugin available right now? And how should I organize the layer with points in order to calculate the road distances between them? I would like the outcome (distances) to be expressed in meters (I am working at the urban scale).


texts: Importing a shapefile or geopackage from QGIS 3.16 to MSSQL Server 2019

texts:
I can only find PostGIS, SpatiaLite and Oracle in the database management to import layers to the respective databases.
However, I can connect to MSSQL Server in QGIS and draw and save geometries in it as well. Manual copying of the attribute table also works.
Why can't I save complete layers or shapefiles in MSSQL with QGIS?


texts: New features to shapefile in QGIS end up with previous features attribute values - why?

texts:
When I digitize a new feature to a polygon shapefile in QGIS the attribute data of the new polygon is replaced by the data entered into the last feature created.  Has anyone else ran into this issue?  Is there a way to fix it?

For example, I may generate a polygon and give it the NAME value of "Gizmo", and dateMade value of "2016-06-12" but when I hit "okay" to confirm it, the complete feature will have a NAME and dateMade value of the previously made feature.  To correct this I have to generate the feature, save it, and then go and edit the feature adding the correct information.


texts: Multiband GeoTIFF loading slow in browser

texts:
Here's the problem: http://www.doctorharapos.com/webMap/geotiff.html

My GeoTIFF files are light, just 4MB
My process: starting from an ASCII Esri file I transform it to GeoTIFF using gdalwarp, because I also do a change of projection
I stack items in groups of 4 GeoTIFF files with gdalbuildvrt
I add LZW compression with gdal_translate and overviews with gdaladdo

An example of the answer of gdalinfo:

Driver: GTiff/GeoTIFF
Files: mP2_apr_p50.tiff
Size is 2978, 1167
Coordinate System is:
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"]]
Origin = (-8.830684460645021,42.286367649318294)
Pixel Size = (0.000057568360165,-0.000057568360165)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -8.8306845,  42.2863676) (  8d49'50.46"W, 42d17'10.92"N)
Lower Left  (  -8.8306845,  42.2191854) (  8d49'50.46"W, 42d13' 9.07"N)
Upper Right (  -8.6592459,  42.2863676) (  8d39'33.29"W, 42d17'10.92"N)
Lower Right (  -8.6592459,  42.2191854) (  8d39'33.29"W, 42d13' 9.07"N)
Center      (  -8.7449652,  42.2527765) (  8d44'41.87"W, 42d15'10.00"N)
Band 1 Block=2978x1 Type=Int32, ColorInterp=Gray
  NoData Value=0
  Overviews: 1489x584, 745x292, 373x146, 187x73
Band 2 Block=2978x1 Type=Int32, ColorInterp=Undefined
  NoData Value=0
  Overviews: 1489x584, 745x292, 373x146, 187x73
Band 3 Block=2978x1 Type=Int32, ColorInterp=Undefined
  NoData Value=0
  Overviews: 1489x584, 745x292, 373x146, 187x73
Band 4 Block=2978x1 Type=Int32, ColorInterp=Undefined
  NoData Value=0
  Overviews: 1489x584, 745x292, 373x146, 187x73

Any clue of what is happening or what can I do to speed up the loading?


texts: Generation of documentation for QGIS plugin

texts:
I've written rather much docstrings etc. in my QGIS plugin, mostly in order to keep track of the functions for myself, but it would be nice to have some more documentation, like a sphinx page.

I have never generated any documentation previously but I know that pb_tool runs sphinx when I deploy the plugin. But the index.html in help\build\html folder is rather empty.

Any suggestions of which settings to update etc.?


texts: gdal_polygonize.py - Silly Size requested

texts:
While running this:

./gdal_polygonize.py GBSM_JKR_ALL_WGS_GEOREFERENCED.tif -f "ESRI Shapefile" to_shapefile/

When the process memory size reaches 32GB it dies with:

ERROR 1: CPLRealloc(-6442451024): Silly size requested.
Segmentation fault (core dumped)

I am trying to convert this file to a shapefile to load it into a database. The .tif file is 101gb. The machine has 3TB of memory. I am (very!) new to GIS systems and have been googling each error.

Is it possible to

run this to exceed a 32gb memory allocation?
chop up the .tif file so that it can be converted piecemeal?


texts: Using negative radius buffer for ST_Buffer?


texts:
I'm trying to put together a command on tippecanoe that will batch-create tilesets that upload as individual layers. All of the layers are keyed to the same spatial boundary, which is in a geojson along with the unique fields for the individual layers. Right now I'm uploading them individually which is taking a long time and I have quite a few layers to create.

Is there a particular command input that I'm missing for this?


texts: Image Collection minimum and maximum Extraction in Google Earth Engine?

texts:
I've wrote this code but it makes error about min and max image?

var modis = ee.ImageCollection("MODIS/006/MOD13Q1")
.filterBounds(geometry)
.filterDate("2000-01-01","2001-01-01")
.select("NDVI");

print(modis);

var mod13 = modis.map(function(img){
  return img.multiply(0.0001)
  .copyProperties(img,['system:time_start','system:time_end']);
});

print(mod13);

var vci = mod13.map(function(img){
 var id = img.id();
 var min =  img.reduceRegion(ee.Reducer.min(),geometry,250);
 var max = img.reduceRegion(ee.Reducer.max(),geometry,250);
 return img.expression(
   "(NDVI-min)/(max-min)",{
     "NDVI" : img,
     "max" : max,
     "min" : min
   }).copyProperties(img,['system:time_start','system:time_end']);
});


texts: Need Projection Help to Buffer Lines in QGIS

texts:
I need to create buffers (1km, 3km, and 5km) around the rivers of this shapefile from NOAA, something I remember being relatively simple in my GIS class.

However no matter what I do the buffers I create cover the entire United States.  I've tried a number of different CRSs (for the layer I did Right-Click>Save as and set the CRS I wanted to try, I also set the project CRS not just the layer), but the buffer always covers most, if not all, of the continent.  Additionally if I use the distance tool to measure this shapefile from coast to coast, it comes out in meters.  I could be wrong, but I'm fairly sure that the country is larger then a couple dozen meters across.


texts: How to use Heatmap plugin and its output?

texts:
I have a layer of points which I would like to depicted as a heatmap concentration of points for visualization.

In QGIS 1.8 I tried the built-in heat map plugin, leaving the default values as is, and got this grey box.

What's the next step after getting this layer produced by the heatmap plugin? I'm hoping to produce something similar to this output from here. http://www.sethoscope.net/heatmap/

P.S.
I actually tried this and the output was beautiful but couldn't get the image as a layer into qgis.


texts: Add cumulative sum column to attribute table

texts: I have an attribute table in QGIS in which I would like to calculate the cumulative sum. The table is sorted by N (number of fatalities) and I would like to calculate the cumulative sum on F (frequency of accident), i.e. F_cumulative should be a sum of the current row and the rows above. How should I do this? Can I use a formula in the field calculator?


texts: How to join MySQL database to QGIS polygons

texts: I Have an *.sql file that has no spatial information. For every record there is an id. I also have a Qgis project with polygons without any id.  How can i connect the polygons with the database? I know that i have to add the same id to the shapes to make the connection. I don't know how to make the connection. Is there any tutorial that can help me?


texts: Using Select Expression and Field Calculator in ONE STEP with QGIS

texts:
I have an Open Street Map Shapefile and i want to use Field Calculator to add information to a field called a800900 according to the field  TYPE of street, by different area (field called ID_2).

The field TYPE of street has the following values

motorway
trunk
primary
secondary
tertiary
residential

And the field ID_2 is a numeric code for 85 municipalities:

8527
8300
8302
...

The information added to a800900 it's based on the TYPE of street and will be different per each ID_2. So, I've been selecting each ID_2 and then performing calculations with Field calculator in two steps.

STEP 1:"Select By Expression"

ID_2 IS '8283'

STEP 2: "Field Calculator"

CASE

WHEN TYPE IS 'motorway' THEN 2405.396*0.01847

WHEN TYPE IS 'trunk' THEN 2305.396*0.0154

WHEN TYPE IS 'primary' THEN 1807.226*0.013

END

The following figure show this process:

Then, i do similar calculations for another ID_2, varying the values in Field Calculator:

STEP 1:"Select By Expression"

ID_2 IS '8302'

STEP 2: "Field Calculator"

CASE

WHEN TYPE IS 'motorway' THEN 2405.396*0.75

WHEN TYPE IS 'trunk' THEN 2305.396*0.34

WHEN TYPE IS 'primary' THEN 1807.226*0.47

END

At the end I'll get values in a800900 different for each TYPE of street and ID_2. As example, if ID_2 is 8300 and TYPE of street is motorway, then a800900 is 2450, if ID_2 is 8302 and TYPE of street is motorway, then a800900 is 2300, and so on.

I have to make the estimation with Field Calculator for each ID_2, being 85 areas, and thats's why i'm looking a faster procedure. I would like to know if there is a ONE STEP procedure with Field Calculator, like:

CASE

WHEN ID_2 IS '8283'AND WHEN TYPE IS 'motorway' THEN 2405.396*0.01847

WHEN ID_2 IS '8283'AND WHEN TYPE IS 'trunk' THEN 2305.396*0.0154

WHEN ID_2 IS '8283'AND WHEN TYPE IS 'primary' THEN 1807.226*0.013

WHEN ID_2 IS '8302'AND WHEN TYPE IS 'motorway' THEN 2405.396*0.75

WHEN ID_2 IS '8302'AND WHEN TYPE IS 'trunk' THEN 2305.396*0.34

WHEN ID_2 IS '8302'AND WHEN TYPE IS 'primary' THEN 1807.226*0.47

END

Is there a way to do this in one step?


texts: Displaying CSV file with centre point data but no XY data?

texts: How do I convert this csv data to points in arcgis? Please see screenshot below.


texts: How to modify mapnik's buffer-size by kosmtik's localconfig.json/.js?

texts: I want to override mapnik's buffer-size in kosmtik. How to override it in localconfig.json/.js?


texts:
I am running ArcGIS 10.2.2 Version 10.2.2.3552, using ArcMap.

I have searched and haven't found anything on this or if it is even possible.  I'm trying to add something to my script that allows the user to define what pages to export themselves, so if they export the maps each week and only need 3 or so pages, they can input this themselves.  I tried throwing something together but anything i throw in there isn't working.  i think this is as close as i get:

Can someone offer me a hint or a nudge if i am in the right direction?

import arcpy
mxd = arcpy.mapping.MapDocument(r"E:\User Folders\Recreation.mxd")
ddp = mxd.dataDrivenPages
pageNumbers = raw_input("Enter Page/s:")
for pageNum in pageNumbers:
    ddp.currentPageID = pageNum
    pageName = mxd.dataDrivenPages.pageRow.Name
    a = ddp.pageRow.getValue('COUNTY_NAM')
    b = ddp.pageRow.getValue('RECREATION')
    arcpy.mapping.ExportToPDF(mxd, r"E:\User Folders\PDFs\\" + a + " " + "(" + b + ")     "+".pdf", image_quality="FASTEST", resolution = "100")
del mxd

Its getting hung up on "The attribute currentPageID is not supported on this instance of DataDrivenPages - even though i have used it successfully in the past.


texts: How to install .tab nativeX driver for QGIS

texts:
How to Install support of extended .tab format in QGIS? I have already installed MapInfo 17.0.3 on my notebook, it supports .tab nativeX format, but QGIS doesn't. But I saw in this web-site
https://support.precisely.com/product-downloads/item/mapinfo-efal-sdk-download/
information about EFAL driver for GDAL 3.0.X and want to try it. Please, explain me what steps I should do to make my QGIS support this .tab extended format


texts: Does OpenGeo Suite give their full capability for free?

texts: Can you help me to decide if I should install Geoserver, PostgreSQL and OpenLayers by myself or should I use OpenGeo Suite as it involves all these? I would need full capability which all these technologies offer when installed outside OpenGeo Suite and I would need to be able to host my own data.


texts: I need to convert a .tif file into a .map format. How can I do it? Is it possible using GDAL?


texts: How to SELECT one day's scene (Landsat 8 scenes)

texts:
I have a huge dataset of raster in PostGIS. Each scene is divided into tiles (tile-size 500x500). I was trying to SELECT rast FROM landsat8 WHERE ST_Intersects(rast, polygon). It runs very slow, because it returns longitudinal scenes of raster. Where I just want one day's scene.
 But the problem is if I query with fixed date, adjacent scene with different date won't be selected. That's not true.
So I think about some solutions:

Because Landsat 8 re-visiting every 16 days, so I will insert every 16 days into different table. So I can query from raster that cover all the earth without duplicate.
I'm not sure about ST_Within(), ST_DWithin() then SELECT rast that doesn't overlap other >50%

What would you do in this situation?


texts: GRASS GIS 7 v.net.iso in the processing toolbox fails to iterate over nodes

texts:
I have successfully created an iso map for a small area of a road network and three points with forward and backward arc costs. However when trying to iterate over the three nodes to produce a separate iso for each, after the first completes I receive the error;

"local variable 'filename' referenced before assignment See log for more details"

And the process terminates. While I could select each point in this example, I am testing for a much bigger exercise. Could anyone point me in the direction of a solution?

Specs - QGIS 2.18.4 64bit, Windows 7 64bit. 2 .shp files stored on a network drive.


texts: Converting multiple .txt files in a folder to shapefile

texts:
I want to convert my .txt files to shapefile. I want the system to be able to move to the next .txt file in the folder, other than me inputting one after the other. This is whatI have done but I have to input the files one after the other which takes a lot of time

 import arcpy

 newfcName ="532724.shp"
 outpath = r"C:/FLOOD/IKEJA AND KOSOFE"

# Declaration
arcpy.env.overwriteOutput=True
arcpy.env.workspace= outpath

# Create new Shapefile and add FIELDS
newfc = arcpy.CreateFeatureclass_management(outpath, newfcName, "Point")
arcpy.AddField_management(newfc, "X", "STRING", field_length = 50)
arcpy.AddField_management(newfc, "Y", "STRING", field_length = 50)
arcpy.AddField_management(newfc, "Z", "STRING", field_length = 50)
arcpy.AddField_management(newfc, "I", "FLOAT", field_length = 50)

# Reference Cursors
cursor=arcpy.da.InsertCursor(newfc, [ "SHAPE@XY", "X", "Y", "Z", "I"])

# Read File
a = open("C:/FLOOD/IKEJA AND KOSOFE/DEM/532724.txt","r")
inputF = a.readlines()

for line in inputF:
xCoordinate, yCoordinate, zValue, iValue = line.split(" ")
xy = (float(xCoordinate), float(yCoordinate))
newRow = (xy, str(xCoordinate), str(yCoordinate), str(zValue), float(iValue))
cursor.insertRow(newRow)

a.close()
infc = r"C:/FLOOD/IKEJA AND KOSOFE/532724.shp"
sr = arcpy.SpatialReference("Minna UTM Zone 31N")
arcpy.DefineProjection_management(infc, sr)


texts: Saving all symbols to style using ArcPy?

texts:
It is possible to save custom symbology to a style file within the Format Symbol pane, but is there an ArcPy command to save symbology to a style file?

I have rich symbology which needs to be transferred to another project and would like to loop over all feature layers (each feature layer has a unique value symbology) which are in a Map, and save each symbol to style.

Below is example symbology of just one feature layer.


texts: While adding a polygon feature in QGIS 2.2.0-Valmiera, the representation symbology of that particular feature(which is in-progress) appears to be "filled" until the completion(rightclick)  of digitization, eventhough that layer symbology is given as 'Outline: Simple line' as shown in attached screenshot. How to fix this?


texts: CartoDB distance measurement/best route along street network

texts: I would like to find out how I can perform a distance measurement between 2 known points along a street network, similar to the way google maps gives directions and distances. Is this possible in CartoDB?


texts: GIS mapping software for visualizing and publishing maps

texts:
Currently I am writing my Bachelor thesis and for testing purposes I need software where I can visualize and publish maps and put vector data (point, line, polygon) on them. Can anybody recommend such software?

I have already founded GIS Cloud. I need something like this.


texts: Accessing the Quality Band for MCD12Q1 with MODISTools GetBands

texts: I am using MODISTools (https://github.com/seantuck12/MODISTools) to process MCD12Q1 data, which, according to the LP DAAC website (https://lpdaac.usgs.gov/dataset_discovery/modis/modis_products_table/mcd12q1), has 16 bands. I believe one of them is a quality band (Land Cover QC). However, when I use the MODISSubsets function to download MCD12Q1 and the GetBands function to view the bands, only 15 bands are returned (without the Land Cover QC band). Without this band, I can't use the MODISSummaries function. Does anyone know how to work around this issue?


texts: I imported polylines and annotations from DGN file. Polylines are pipelines. Annotation layer contains pipe diameters. There is no common ID's. The goal is to assign diameters from annotation to nearest pipe. Anybody knows how to do that?


texts: Divide polygon evenly based on lines and intersection points

texts:
I have this purple polygon illustrated in the figure below. I would like to split the polygon so that each line segment (green) split the polygon evenly with the surrounding line segments. A line segment is split at intersection points.

I have tried to illustrate by painting a red line where I think the polygon should be divided.

Is this a normal and simple process to do or does this require a lot of working?

The terms Voronoi and Delaunay comes to my mind, but I'm not sure if that is what I want here.

Edit:
I'll try to describe how I think of evenly in this case. I want each green line segment to take its share of the purple polygon. It is equal to buffering the green line segments, but instead of creating one big polygon (just like the purple one) each buffered line segment continue to grow until it meet another one - just like inflating a lot of ballons in the same room, they can grow as long as there are vacant space.


texts: Creating route and adding that route as layer to map using ArcPy?

texts:
I wrote a little script for ArcGis, and when I make it run by GUI it happens that it seems to create routes but the routes do not appear in the TOC.

Consequently, when the script tries to execute operation on those routes, the dataset cannot be found.

How can I solve this issue?


texts: Finding zipcodes within line using GeoPandas

texts:
I have a shapefile of Houston zip codes and gdb folder of lines (within Houston).
The lines are in a gdb folder containing over 125 files. The types of files included in this folder are GDBINDEXES, GDBTABLE, GDBTABLX, ATX, SPX, and etc. When I add this whole folder as a layer on ArcMap I get a series of different lines and ArcMap labels the layer as edges. Each line starts at one coordinate and ends at another coordinate.

I want to figure out which zip codes each line crosses using GeoPandas. I have over 500,000 lines. I have no previous experience using GeoPandas.


texts: Why doesn't r.kappa work for me in QGIS 3.2.16?

texts:
I have 2 rasters reclassified from 1 to 6 from 2 interpolated IDWs. I apply the following steps in r.kappa:

I select the raster of the prediction in classification result
The other base raster that contains the reference classes
I do not set the extent of the GRASS GIS 7 region (optional)
Size of cells in the GRASS GIS region 7 = 10
Error matrix and kappa = save as file
Execute
I verify in records where the txt file is saved (confirm that it was saved where indicated)

The problem is that I do not generate the txt in the assigned folder.


texts: Reading tile package file(.tpk) in ArcObject C#?

texts:
I wish to read and show some tile package file using ArcObject C# 10.x.

However, I failed to find how to read an tile package file on My Application using ArcObject C#.

Can you tell me how to do this?


texts: Geoserver GeoWebCache caches tiles of incomplete layers

texts:
I have data in a PostGIS database that is being pulled into Geoserver using a PostGIS (JNDI) store. The dataset I am using is a continuous polygon grid, but the dataset itself isn't all that large. When I try to have Geoserver's GeoWebCache cache this layer, at a certain zoom level, the layer appears to only be partially complete at the time the GWC "takes a snapshot" of the data. If I zoom in or zoom out one level, the layer tiles appear correct.

The layers are stored as 4326 in PostGIS, but are served as 900913. I have tried converting the table geometry to 900913 so there wouldn't need to be any re-projection, but this has no effect on the cache tiles. I have also tried caching using 4326, but those tiles have the same incomplete layers.

The dataset has a geometry index and the table is clustered on the index.

Is there some kind of optimization I'm not thinking of that would allow GeoWebCache to cache the complete layer instead of the incomplete one?

EDIT: I tried changing the styling to the default "polygon" style, and the layer appears correctly when cached. So it doesn't appear to be an issue with the geometry, it is instead an issue with the styling. I tried making the styling as simple as possible by not needing Geoserver to do any searching for values between X and Y, but instead replacing bins with integers. I also applied an index to the column being used for styling. So all Geoserver really needs to do is style bin_val = 1 blue, bin_val = 2 red, and so forth, but even that isn't returning a complete cached layer.

For example, this is what the cache looks like at zoom level 5:

And this is zoom level 4, notice the Alaskan data is only the outermost grid polygons, while the more inshore ones are not visible:

Geoserver version 2.4.5


texts: Create default ArcToolbox parameter with current folder location of ArcToolbox .tbx

texts:
I have an ArcToolbox tool I created based on a script I wrote that needs an input parameter that equals the folder location of the ArcToolbox .tbx file from which the tool runs.  I am distributing the toolbox, so users can run it from a variety of locations.  How do I create a folder parameter that autofills with a default of the current .tbx location?

I have tried writing a bit of python under the Validation tab of the Tool properties:

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

    stringFilter = self.params[1].filter

    list = []

    pathname = os.path.realpath(".")

    list.append(pathname)

    stringFilter.list = list
    if not self.params[1].altered:
        self.params[1].value = stringFilter.list[0]

return

But when I run the tool, the only text that autofills into the parameter reads, "RemoteDatabase"


texts: How to show label over markers on leaflet?

texts: I need to show label over markers on leaflet to make it readable. How can I implement it? See the picture below label is overlay-ed by marker


texts:
I'm using QGIS LYON 2.12.1, and I want to turn a plain gray scale raster file (that I got from Google search) into a DEM.

The aim is to apply the DEM (terrain models) function to the newly created DEM.

How would I convert the .tiff or .jpeg into a DEM with elevation data?


texts: Correct method to save and restore leaftlet drawing into mysql database

texts:
I have the following code for saving my coordinates of what I draw. So what I do on drawing created I have the following function.

 map.on(L.Draw.Event.CREATED, function (e) {

        var type = e.layerType,
                layer = e.layer;
        var points = JSON.stringify(e.layer.toGeoJSON());

        console.log("coord : "+e.layer.toGeoJSON().geometry.coordinates);
        var latlngs = [];
          latlngs = LatLngsToCoords(e.layer._latlngs, 1);

        if (type === 'marker') {
            layer.bindPopup('A popup!');
        }
        if (type === 'polygon'){

        console.log("coord : "+e.layer.toGeoJSON().geometry.coordinates);

        }

        drawnItems.addLayer(layer);
    });

This output e.layer.toGeoJSON().geometry.coordinates gives me example so I will then transfer it into 175.28162956237796,-37.77007019461929,175.28506278991702,-37.77271618103962,175.27849674224854,-37.77458188377217,175.27781009674072,-37.77376776473304,175.28162956237796,-37.77007019461929

into POLYGON(( 175.28162956237796 -37.77007019461929,175.28506278991702 -37.77271618103962,175.27849674224854 -37.77458188377217,175.27781009674072 -37.77376776473304,175.28162956237796 -37.77007019461929 ))

I know I could do this and save into the MySQL database.

For the reverse when I need to restore the polygon I will need to extract the POLYGON and now reverse the latitude and longitude back.

I find there is a lot of steps. Is there any direct method ?


texts: Can't create time-series raster WMS with Geoserver due to time regex expression issue

texts:
I am attempting to create a raster time-series WMS server using GeoServer and have been working off this tutorial provided here:

http://docs.geoserver.org/stable/en/user/tutorials/imagemosaic_timeseries/imagemosaic_time-elevationseries.html

Even though I created the timeregex.properties file as indicated then loaded the store and enabled the time dimension, the WMS getCapabilities shows only that the YYYY-MM-DD portion of the filenames were ingested and any HH:MM:SS in the filenames were not ingested.
Similarly, the shapefile created only lists the YYYY-MM-DD in the "ingestion" field and not any HH:MM:SS
I need to create a regex expression for filenames that have YYYYMMDDTHHZ, but I am not finding any resources on how one would create these regex expressions.


texts:
In QGIS, I added a WFS-Layer from my geoserver, containing few point data with a numeric attribute (Z-values).

This works fine, and all the points display on the map.

Now I just assumed, I would be able to do a quick IDW Interpolation from the Raster -> Interpolation menu. But in the interpolation ui, I can't choose an "Interpolation attribute" from my wfs layer, because the dropdown is empty.

If I do the same from an csv file containing the same data, it obviously works fine.

Has anyone an idea, why this shouldn't work? Do I miss the point here?
Is it because of this operation is not supported for wfs layers or is it about my adjustments?


texts:
I am trying to sum the number of fires in each state of the US. I have a shapefile for the states and I have a shapefile for each fire. Vector->Analysis Tools-> Count points in a polygon isn't working because I keep getting this error but I don't understand it. Does anyone know how to fix this?

line 211, in spatialindex
            idx = QgsSpatialIndex(layer.getFeatures(request))
        Exception: unknown


texts: Adding second ArcGIS Online Basemap to ArcGIS Pro Map?

texts:
In ArcGIS Pro 2.3.0 I would like to be able to switch ArcGIS Online Basemaps at a particular scale.  This should be easy to do by using complementary Visibility Ranges on two basemap layers in the Contents pane.

However, if I use the Basemap button to access the Basemap gallery, and choose another basemap, the one I choose replaces the original.

How can I add a second ArcGIS Online Basemap to an ArcGIS Pro Map?


texts: How to determine the most accurate elevation algorithm

texts:
I use an Android GPS cycling app to track my bike rides and I would like to get a reasonably accurate estimate of the total ascent I did on a ride.  I understand that GPS is very inaccurate at determining elevation, so I have been importing the resulting GPX files into several sites that use mapping data to determine elevation, and I've found they're ridiculously inconsistent.  For example, from a 30-mile ride I get the following climb data from four sites using the same GPX file:

CardioTrainer   980
Strava          1085
MapMyRide       554
RideWithGPS     2131

The results are consistent over multiple rides, and I can add more sites with an equally mixed set of results.

I think I can easily toss RideWithGPS as a badly flawed outlier based on subjective experience alone, but MapMyRide is consistently about half the ascent of the others, so it makes me wonder who is off by 2-to-1.

My question is, how can I determine whose algorithm is the most accurate without spending hundreds on a device with a barometric altimeter?


texts: GDAL reports no bands for a 12-band raster and no projection either

texts:
When opened in ArcGIS, the raster correctly unfolds 12 bands, but when I try to open it with GDAL through Python, it appears as a raster that has zero bands.

Here is the code:

import gdal
from gdalconst import *
filename="Images/Image.hdf"
dataset = gdal.Open( filename, GA_ReadOnly )
print ("Format: ",dataset.GetDriver().ShortName)
print ("Projection: ",dataset.GetProjection())
print ("Bands:", dataset.RasterCount)

And here is the output:

('Format: ', 'HDF4')
('Projection: ', '')
('Bands:', 0)

As you can see, the projection is reported as empty, and the number of bands as zero.
The raster format is HD4 or HDF. You can easily download a similar raster from this FTP site:
ftp://ladsweb.nascom.nasa.gov/allData/5/MOD11A1/2012/193/

Do you have any idea about this issue?


texts: QGIS 2.2: GDAL Grid (Interpolation) Error

texts:
As much as i understand Raster -> Analysis -> Raster (Interpolation) can provide me with a terrain file out of a shape file with a Z coordinate.
I select the input and output file and a Z field.
It gives me this error:
ERROR 6: GDALDriver::Create() ... no create method implemented for this format.

Unable to create target dataset "D:/shp/terrain.dem".
GDALDriver::Create() ... no create method implemented for this format.

Tried some other formats (I'm aiming for either .dem, .hgt or .xyz) but got the same error everytime. What am I missing here?

The Vector layer contains X coord, Y coord, Z coord, number and name fields.


texts: Problem with connecting remote WMS connection to GeoServer

texts:
In my geoserver admin panel I try to connect new stores - remote WMS connection:

http://mapy.geoportal.gov.pl/wss/service/pub/guest/kompozycjaG2_TBD_WMS/MapServer/WMSServer?service=wms&version=1.3.0&request=GetCapabilities

and when I try to save this store, displays errors:

Connection test failed: Content type is required for org.geotools.data.wms.response.WMSGetCapabilitiesResponse

Does anybody knows what could be reason?


texts:
I can't seem to port this line to QGIS 3.0 API

  QObject.connect(self.dlg.ui.chkActivate,SIGNAL("stateChanged(int)"),self.toggle_activate)

Any ideas? I believe QObject is deprecated.


texts: Extract time series values from a 3D (lon,lat,time) NetCDF file using Python

texts:
I want to extract time series from a variable in a 3D (lon,lat,time) netcdf file at specific lon/lat points. The raster is a nc file where:

cell: 1.5°×1.5°
time step: 1 month
lat extension: (0.75,89.25,1.5)
lon extension: (0.75,359.25,1.5)
time: (201001,205012)
step scale: cell center

I have a point shapefile (lat/lon) where a point is (108.3°,58.6°). I just want to extract the cell values of the whole time steps in where points located. Are there directly method to carry out it？

I think a way that round the points number according to lat/lon step: changing the P lat/long to cell lat/lon where point located. For instance, P(108.3°, 58.6°) is in cell (108°,60°) because scale step is 1.5°. Now, the question is how to change P(108.3°,58.6°) to cell(108°,60°) used to slice .nc file.


texts: Changing Mapbox popup based on null values

texts:
I currently have a popup in Mapbox that displays some properties and a photo. Not all of the points have photos though, so the ones that have a null value in the image link field just show a blank box where the photo would be. I want to set it so that any features that have a null value in this field only show the text in the popup and not the image link part.

Here's the popup code I have so far:
  var popup = new mapboxgl.Popup({ offset: [0, 0] })
    .setLngLat(e.lngLat)
    .setHTML('
' + feature.properties.Name + '
' + feature.properties.PopupInfo + '
 ' + '
')
    .addTo(map);

I found this question that seems to address the issue, but I'm fairly new at Javascript, and I'm having trouble understanding how exactly I would implement this in my scenario.
https://stackoverflow.com/questions/55106852/show-not-null-attributes-in-mapbox-gl-js-popups
Here's the code given in that example:
const feature = {
  properties: {
    attr1: 'text',
    attr2: null,
    attr3: 42,
    attr4: null,
    ...
  }
};

const popupContent = Object
  .keys(feature.properties)
  .reduce((acc, property) => {
    const value = feature.properties[property];

    // your condition here
    if (value) {
      acc.push(`${property}: ${value}`);
    }

    return acc;
  }, [])
  .join(', ');

console.log(popupContent) // attr1: text, attr3: 42

I think I'd have to replace the .join part with .concat and add in the HTML text I want and then pass that back into .setHTML in the popup var, but I can't figure out how to pass my properties into this and then how to get it back into the mapbox.Popup part.


texts: Reading big raster getting warnings using GDAL Python

texts:
I am trying to use GDAL Python API to deal with some big raster, and when I read these rasters into Python, I get warnings like that:

Warning 1: TIFFFetchNormalTag:ASCII value for tag "GeoASCIIParams" contains null byte in value; value incorrectly truncated during reading due to implementation limitations
Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.
Warning 1: TIFFFetchNormalTag:ASCII value for tag "GeoASCIIParams" contains null byte in value; value incorrectly truncated during reading due to implementation limitations
Warning 1: TIFFReadDirectory:Sum of Photometric type-related color channels and ExtraSamples doesn't match SamplesPerPixel. Defining non-color channels as ExtraSamples.

Do anybody encounter problems like this? Why did this thing happen?

The raster I use is about 2.2G for one raster, and the raster' format is TIFF.the warnings don't affect the later use of the data, it is still an annoying issue, I want to know how to fix this.


texts: Can I have Points, Lines and Polygon features in one QGIS layer?

texts:
I have a situation where planning is done on a "base map" consisting of rivers, contours, roads etc. There is a need to create a "overlay(s)" over the "base map" consisting of drawings, that could be point(s), line(s) or polygon(s) on one overlay.

Is it possible to create point, line and polygon features on one QGis layer/overlay
I have looked at the vector overlay with a memory data provider but it looked like you can only create a instance of the memory data provider for one type of feature only?

I would appreciate any help, thank you.


texts: Selecting a specific landcover type, filtering date and clipping for MODIS land cover data

texts:
I am trying to select type 5 land cover from the MODIS 500m dataset for the year 2001, and then extract out the projection for that layer. However, my script doesn't seem to work because I still obtain a 'collection' when I am expecting an image (since there is 1 dataset per year and I am selecting a specific layer). I know it doesn't work because when i try to extract the projection using projection(), it gives me the error, 'modis.projection is not a function'. I know it works if i put .first for the image collection, but I want to apply this to 2005 as well (and understand where I went wrong).

Another thing I noticed is that the .clip function and .filterBound doesn't seem to work with MODIS (even when I use the .first function. Why is that so?

Thank you in advance.

// Script that doesn't work and gives the error
var modis = ee.ImageCollection('MODIS/006/MCD12Q1')
           .filterDate('2001-01-01','2001-12-31')
           .select('LC_Type5')
           .clip(geometry);

print('MODIS projection:', modis.projection());

// Script that works but doesn't end up clipping the image
var modis = ee.ImageCollection('MODIS/006/MCD12Q1').select('LC_Type5')
            .first()
            .clip(cprBasin);


texts: Dateline disappears after clipping shapefile to hemisphere?

texts:
I am working with this natural earth shapefile. The dateline disappears after using the mentioned clip to hemisphere plugin. Densifying the geometry does not help.

Weird to me seems that clipping to a hemisphere centered at lat -76,98267 / lon -163,36119 works just fine (dateline present after clip) while after clipping to a center at lat 26,50588 / lon -171,82873 the dateline dissaperes although it should be visible.

Solved as noticed in the comments  by @Jake a bug in the clip to hemisphere plugin is the problem. I quickfixed it and its all good now.


texts: How to know stream order in QGIS on output of r.stream.extract

texts:
I have stream order created over DEM using r.stream.extract , which gives 3 output files:

stream_raster - raster map stores extracted streams.
stream_vector - extracted points for each stream segment
flow accumulation

Can I know which segment comes after each segment in QGIS.
Are these stored with any IDs somewhere so that I can know flow direction or next point from stream_vector for next flow segment


texts: GNSS planning online vs trimble planning 2.9 software

texts:
Why using the same parameters GNSS planning on-line and Planning 2.9 Trimble show me different Dilution of precision (DOP on the graph)?! Which program is better mission planning solution for GNSS measurement?

Link to full size image.


texts: ArcGIS JavaScript API Drawing Graphics to Wrong Places

texts:
I am working on a web application using ArcGIS JavaScript API. I used query and got a group of graphics. However, the map.graphics.add draws the graphics to wrong locations. They should be in Arizona but are located to west coast of Africa. I checked this article: https://stackoverflow.com/questions/17058843/arcgis-javascript-api-need-to-change-wkid-of-the-graphics-layer
and used the method mentioned. It seems the method doesn’t apply to my situation. Even when the map, graphics layer and graphics all have the same spatial reference, the graphics are drew to Africa. Do you know why I get this weird issue? Here is the link of the web application:
http://www.public.asu.edu/~jzhan300/HCIC/HCICInteractiveMap.html
When you click the select button and then hit apply in the lower left corner, it should draw graphics.
The link of the js file:
http://www.public.asu.edu/~jzhan300/HCIC/HCICInteractiveMap.js


texts: Use ArcGIS to classify intersection types

texts:
I am using ArcMap 10.4
I have a shapefile of a road network. Each road has a functional class value in the 'FC' field. I'm trying to determine how often roads of each functional class intersects another.

I Exported a shapefile of each functional class (Definition queries interact oddly with intersect).
As some segments consisted of multiple bits, I Dissolved each, and then ran Multipart to Singlepart to get a file where each segment was one feature.
Some of the roads are divided highways, so I've used 'Merge Divided Roads' and 'Collapse Dual Lines to Centerline' at various distances. Neither has been completely satisfactory, but oh well.
Then I've run Intersect between the two shapefiles, setting the output to point.

But I'm uncertain how to get the data on which lines it intersects back to the point, so as to classify the point by the types of functional class of the roads intersection. Past process I've used has been to use Select, and field calculate four columns, but seems like their ought be a better way.


texts: How to move node position keeping right angles?

texts:
I made a polygon with perpendicular segments using CAD Tools.  I saw that, even if the segments of this polygon are perpendicular, the four points has different values in terms of x, y. I think that`s because of the projection.

The thing is that now I want to expand the polygon but keeping the lower-right point and also, right angles. I duplicated it and I tried to keep the right angles but, manually, it`s impossible. How can I do that ?


texts: raster.kendall function in R messes output of very high resolution raster stack

texts:
I use R in order to calculate the Sen Slope trend for a raster stack (netcdf file) for each pixel, using the following code, from the SpatialEco library:

#
raster.kendall <- function(x, intercept = FALSE, p.value = FALSE, confidence = FALSE, tau = FALSE, ...)
  {
  if(!any(class(x) %in% c("RasterBrick","RasterStack"))) stop("x is not a raster stack or brick object")

  if( raster::nlayers(x) < 5) stop("Too few layers (n<5) to calculate a trend")

  trend.slope <- function(y, p.value.pass = p.value, tau.pass = tau, confidence.pass = confidence,
                          intercept.pass = intercept) {
    options(warn=-1)
    fit <- EnvStats::kendallTrendTest(y ~ 1)
    fit.results <- fit$estimate[2]
    if(p.value.pass == TRUE) { fit.results <- c(fit.results, fit$p.value) }
    if(confidence.pass == TRUE) {
      ci <- unlist(fit$interval["limits"])
      if( length(ci) == 2) {
        fit.results <- c(fit.results, ci)
      } else {
        fit.results <- c(fit.results, c(NA,NA))
      }
    }
    if(intercept.pass == TRUE) { fit.results <- c(fit.results, fit$estimate[3]) }
    if(tau.pass == TRUE) { fit.results <- c(fit.results, fit$estimate[1]) }
    return( fit.results )
  }
  k <- raster::overlay(x, fun=trend.slope, ...)
    names(k) <- c("slope", n)
  return( k )
}

I have tested the above code with many datasets of various resolutions and it worked fine. However, I lately used it with a very high resolution satellite precipitation product (0.0375 x 0.0375) and the results did not make any sense at all.
This is how the raw precipitation data look like:

This is how the output (trend slope from raster.kendall) looks like:

I am reading the netcdf files in R as a raster stack and the code worked fine for all datasets with resolutions up to 0.05 x 0.05, however in this one, results are messy and I am wondering if its something that has to do with memory. A sample .nc file that I use can be found here.
library("raster")
library("maps")
library("rgeos")
library("grid")
library("ncdf4")
library("maptools")
library("mapdata")
library("rgdal")
library("lattice")
library("Kendall")
#
files = list.files(path="/my/path/.../", pattern="*.nc", full.names = T, recursive = T)
#
rast = stack()
trend_rast = stack()
trend_rast_mask = raster()
#
for (i in 1:length(files)) {
  print(files[i])
  rast = raster::stack(files[i])
  print(class(rast))
  file_name = print(files[i])
  trend_rast = raster.kendall(rast, intercept = F, p.value = T, confidence = F, tau = T)
  names(trend_rast) = c("slope","p.value","tau")
  trend_rast_sen = trend_rast[[1]]
  writeRaster(trend_rast_sen, paste("/my/output/.../", file_name, "_SEN_SLOPE", ".nc", sep=""), format="CDF", overwrite=TRUE)
}


texts:
I try to use OSM power lines poles vector data to render 3d power line poles into 3D gis program osgearth.

So in OSM data there are power line poles as points and power lines as lines. So I use point coordinates to place poles, but I don't know orientation of pole. Also I have lines, but theese lines are not bounded to poles. Another way is use OSM lines data to populate poles across line, but there is disadvantage that pole can be placed in river or house.

So I think of using for example QGIS to bound line data to point data, or calculate closest point to specific point and write it to point attributes  so if I have two points I can calculate azimuth of pole. How can I use QGIS to achive this?


texts: Invalid character in UTF-8 when geocoding in QGIS

texts: I am using QGIS 2.18. As suggested, I installed the MMQGIS plugin to help me geocode the csv file containing a list of addresses. I converted the file to utf-8 encoding format. But when i was trying to import the file, it said my file contains invalid characters. I checked my file and could't spot any.


texts: Not possible to print GeoServer layers with CQL filters

texts:
I use geoexplorer together with GeoServer and mapfish.
I make a WMS request passing in it two layers with corresponding CQL filters. I make this request through Geoexplorer to GeoServer.
This works normally if the CQL filters are different to each other. But if apply the same CQL filter to both layers then its not possible the layers to get printed (on mapfish).

For example if I have the following CQL filter for both layers:

iso3=NGA

the layers won't even printed on the mapfish output (the pdf).

If my CQL filter for one layer is:

iso3=NGA

and for the other

iso3=COD

then it will work and the layers will be printed correctly in the pdf.

When I check the geoserver log I see the following error:

  org.geoserver.platform.ServiceException: Layers and filters are mismatched, you need to provide one filter for each layer
      at org.geoserver.wms.map.GetMapKvpRequestReader.getFilter(GetMapKvpRequestReader.java:507)
      at org.geoserver.wms.map.GetMapKvpRequestReader.read(GetMapKvpRequestReader.java:362)

For some reason geoserver or geoexplorer or mapfish is merging the filters, which are identical. Or at least I assume that.

Does anyone has any idea what's going on?


texts: Transform RGB-Map to CMYK (qgis, opensource tools)

texts: I want to print a map I made with qgis. For the print I need CMYK-colors but the print-composers pdf-output is RGB. How can I convert the RGB-map/-pdf to CMYK?


texts:
I have created a map from geojson file of administrative boundaries (polygons layer).

Similar to this example. When the mouse is over the polygon feature, the cursor changes from hand to finger. In the finger mode, map does not zoom on double click.

In my case, at higher zoom, whole map is made of adjacent polygon features so the mouse is always in the finger mode. I need to be able to zoom in on double click.


texts: How can I use py2exe with the OSGeo4W Python installation?

texts: I have a default Python installation in c:\Python27. I have also installed Python through the OSGeo4W installer, which has created a "sandboxed" Python installation with all the geospatial bindings I need. Now, I need to package a script using py2exe, but I can't figure out how to install or configure py2exe to work with the OSGeo4W installation. I have configured OSGeo4W to install external packages with pip, but py2exe doesn't seem to install with pip. If I try to use the downloaded installer, it only finds the system install in the registry. Any ideas how to get it working?


texts: Large Label Expression Causes Application Error

texts:
I'm having trouble with a large label expression and I keep crashing with an application error. I'm trying to display sample results as stacked:

Label ID

Sample depth 1

Sample Depth 2

Sample Depth 3

Now, I also want each of those depths to have three criterias. (Regular, Bold, or RED) Depending on what the value exceeds.
So for example:

SB-1

5

10.5

100(red)

So I've joined an excel table with the sample results

As you see I have a depth field followed by the next three cells being the results based on their criteria. Could not figure out a way to color them based on their value. (<9.5 , >9.5, >45)

This is the VBScript I am currently using.

[MAP_ID] &vbnewline & [F5] & "<BOL>" & [F6] & "</BOL>" & "<CLR red = '255'>" &  [F7] & "</CLR>" &vbnewline & [F9] & "<BOL>" & [F10] & "</BOL>" & "<CLR red = '255'>" &  [F11_1] & "</CLR>" &vbnewline & [F13] & "<BOL>" & [F14] & "</BOL>" & "<CLR red = '255'>" &  [F15] & "</CLR>"

I'm sure there is an easier way to do this...

I did not have a problem with the same expression with only one depth as you can see here

I've tried changing blank values from null to " ", thinking this could be the issue, but no different result, still crashed..

Update, I have exported to shapefile so the attributes are no longer pulling from excel. Also, I have found that if i replace the blank spaces with a value (i.e. "-" or "NA") I do not have this same issue. Except the problem with this is that now I get labels like:

SB-1

5 - -

-10.5 -

--100(Red)

It seems that the spaces tend to be the problem. Maybe an If Then statement is necessary?


texts: How to write syntax to sum values for unique IDs in Survey123 Summary Report?

texts:
How do I get a sum of values for a question for each uniquely named pump station in a Survey123 summary report?
I've developed a workflow for my client to do Pump Station inspections using Survey123.  They can do as many inspections per pump station as they want in a given month.  Each month they compile the inspection data in a monthly report (See screenshot below)

 The report compiles the monthly inspection data for each pump station and provides a sum of the values collected for the specified date range (the date range is chosen when running the report in Survey123)
To recreate this report in Survey123 using the Survey123 inspection data, it seems I need to create a summary report. But I'm struggling with the syntax especially when there are multiple inspections for one pump station in a month.  Using online forums etc.  I was able to get the correct syntax to only show the unique pump station names, but I'm unsure what syntax I should be using to say, for example, "sum all the pump 1 run hours for each unique pump station name."
My working syntax is below.  Again, I was able to get the Pump Station name to show only the unique names and to not duplicate the names when there are more than one inspections done on a pump station.  But I can't figure out to sum each question for each pump station. In addition to that, I also only want to show the most recent inspection date for each pump station.
Pump Station Name:
${#PumpStationLogs | stats:"count,FACILITYID":"PSNAME": "" | where:"1=1" | orderByFields:"PSNAME ASC" }${PSNAME} ${/}

Last Inspection Date:
${#PumpStationLogs | stats:"count,FACILITYID":"PSNAME": "" | where:"1=1" | orderByFields:"PSNAME ASC" }${REPDATE | format:"MM/DD/YYYY"} ${/}

Sum of Pump 1 Run Hours:
${#PumpStationLogs} ${P1_RNHrs} ${/}


texts:
I have two lake boundary shapefile from two different years. The size of the lake decreased over the years and I would like to generate the boundaries of the years that are between what I have.

To achieve it I need a technique that looks like the following:

The black and green polygons represents the ones I have. The red is the what I want first, the Letters should indicate, that the red polygon has the same distance from the black and green.
If possible I'd like the solution in ArcMap (I have advance license).


texts: Adding table to existing map layer in ArcGIS Online

texts: I have a map service that is published to a local ArcGIS Server. I've added it as a service to ArcGIS Online and added it to a map and made a lot of changes to the searching, pop-up's, etc. Later, I went back to the MXD and added a new table that is referenced by another layer. When I go back to the map in ArcGIS Online, it does not show the table I added. However, if I start a new map and add the service, the new table is picked up. How can I add that one table to the existing map without having to delete the layer and add it back? I really don't want to have to set up the map over again every time I need to add a layer.


texts: How To GetFeatureInfo From a GeoWebCached Layer in Geoserver

texts: I'm really really struggling with this.  I have created a map using GeoServer and OpenLayers.  Things have been working out here and there, as I am new to all this so a lot has been just trial and error.  I noticed my map started loading slower and slower so I have decided to use GeoWebCache to cache my layers.  The problem is that now I cannot get the feature info to display in a pop-up as it once was when I was getting the layers straight from the server/wms.  The part that even confuses me further is that when I provide the pop-up with the URL 'http://localhost:8080/geoserver/gwc/service/wms' it actually attempts to send a request.  However when I provide the right URL from the server where the data actually is (http://localhost:8080/geoserver/wms), it's like if no request is being made.  (I am seeing this in firebug).  I've also created a geowebcache.xml files and have listed several of my layers as <queryable>true</queryable>, but no differences there either :(  Any suggestions you can offer me will be greatly appreciated.  I am truly very stuck.


texts: Converting Esri ASCII raster to Shapefile for CartoDB?

texts:
I have an ArcGIS  grid ascii file (Ramankutty Cropland at 5 minute resolution) that begins as follows:

ncols         4320
nrows         2160
xllcorner     -180
yllcorner       -90
cellsize      0.08332825
NODATA_value  -9999
-9999 -9999 -9999 -9999 -9999 -9999 -9999
...

Someone gave me this file but did not provide a projection and could not recall what it was. However, I have some reason to believe the projection is WGS1984:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

I would like to be able to convert this to a format that CartoDB could use--an ESRI shapefile. Another person who did manage this, but who also could not recall the procedure, created a polygonal coverage of grid squares (closed oriented loops). The values at the row and column numbers were assigned to the polygonal squares (somehow--I do not know the data structures involved).

Is there a way to create some kind of polygonal grid square coverage without using ArcGIS (I'm guessing this was used, but no details were provided to me)? Open source tools would be nice, but if ArcGIS is necessary, so be it.

As an alternative, I tried
    gdal_translate -of GTIFF -a_srs WGS1984.prj ...
and uploaded the result to CartoDB, but CartoDB complained that the result was not georeferenced (the_geom lat/lon GeoJSON was empty).

(I'm not sure how to apologize for asking such vague and uninformed questions.)

Some more info: here's part of a row from the ESRI shapefile of this data that was imported into CartoDB

**id**  1
**the_geom** {"type":"MultiPolygon","coordinates":[[[[-8.7649,37.34027], [-8.87019,37.36377],[-8.79447,37.42393],[-8.7649,37.34027]]]]}
**shape_area** 63274611.152
**shape_leng** 36270.0549322

Another row:

**id**  2
**the_geom** {"type":"MultiPolygon","coordinates":[[[[-8.59219,37.34027],[-8.67854,37.34027],[-8.67854,37.40889],[-8.59219,37.40889],[-8.59219,37.34027]]]]}
**shape_area** 92417860.2
**shape_leng** 38453.683350300001621


texts: Changing point icon (geoJSON) on click and changing back when clicking another point using Leaflet?

texts:
I have a geoJSON (10 points) that is drawn using an icon. I want to change the icon when I click on one of the points and then I want to change it back when clicking another point. I don't know how to change it back.

Can anyone help me?

I have this piece of code:

var icon1 = L.icon({
    iconUrl: 'img1.png',
    iconSize: [20,20]
});

var icon2 = L.icon({
    iconUrl: 'img2.png',
    iconSize: [20,20]
});

function clickFeature(e) {
   var layer = e.target;
   layer.setIcon(layer.options.icon = icon2);
}

function onEachFeature (feature, layer){
    layer.on({
        click: clickFeature,
    })
};

var pointsJSON = L.geoJson(points,{
        pointToLayer: function (feature, latlng) {
        return L.marker(latlng, {icon: icon1});
    },
    onEachFeature: onEachFeature
    }
).addTo(map);


texts:
Does anybody have any idea about SuperMap GIS?

What are its limitations and strengths?
How does it compare to other GIS (especially ArcGIS)?
Can a local Government Department rely on it?


texts:
I have a problem wit my current project. When I export my QGIS Project via qgis2web there is a strange behaviour of the map pop ups which are not fitting the content.

It sometimes cuts the container size, but I can not find the reason. It happens randomly. Sometimes it is displayed properly sometimes it looks like on the attached picture.


texts: Loss of synchronization between map and markers (Leaflet + Google Maps)

texts:
I have a problem using the markers with Google Maps with Leaflet. In this example you can see the difference between Google Maps or OSM.

Using OSM, the marker moves with the map but when using Google Maps, it seems that the movement is not synchronized.

Another example of the same problem

Is there any known solution to this issue?


texts: ST_Area yields total area encompased by a polygon. How do I get just the area covered by a polygon?

texts:
Here's my dataset: https://aaron-bauman.carto.com/dataset/curbs_1

I want to calculate the area covered by the purple only.
ST_Area seems to be yielding all the area within the purple as well because when I sum it up the value is greater than the entire area of Philadelphia.

How do I get the area of just the purple?

The dataset bunch of polygons outlining street surfaces. The polygons cover huge areas. Here's a small snippet of one polygon for example:

Here's the query i'm running:

select sum(shape__area) from public.curbs_1;


texts: Accessing geometry with SearchCursor in ArcPy?

texts:
I'm trying to write a script that clips all feature classes in fcList by the successive rows of a clip feature. I'm using arcpy with ArcGIS 10.0

The problem is that I need to verify whether the current fc of the fcList overlaps the current row of the clip rows. Here is the code (I commented the lines with my doubts, questions, remarks):

import arcpy, os
arcpy.env.workspace = r'Z:\Documents\SIG\tests\synthese_tocorrect.gdb'
fcList = arcpy.ListFeatureClasses() #the feature class to be clipped by each of the polygon of the clip layer
pochoir = r'Z:\Documents\SIG\tests\decoupe.shp' #the clip layer with multiple polygons
rows = arcpy.SearchCursor(pochoir)
count = 0
outFolder = r'Z:\Documents\SIG\tests\sorties'
for fc in fcList:
    for row in rows:
        p = row.Shape #I'm trying to get only the geometry of the current row but I thing it's not working
        arcpy.MakeFeatureLayer_management(p,"ptmp") #...the output is empty while it should be the first polygon corresponding to the first row of the clip layer.
        arcpy.SelectLayerByLocation_management("ptmp", "INTERSECT", fc)
        arcpy.CopyFeatures_management("ptmp","ptmpselect")
        nbrow = arcpy.GetCount_management("ptmpselect")

        if nbrow > 0:
            secteur = row.secteur #The integer ID of each polygon of the clip layer
            out_poly = (outFolder, fc + "_" + str(secteur) + "_" + string(count))
            arcpy.ClipAnalysis(fc, "ptmpselect", out_poly)
            arcpy.DeleteManagement("ptmp")
            arcpy.DeleteManagement("ptmpselect")
            count = count + 1
        else:
            secteur = row.secteur # I know I repeat myself. Any idea to avoid it ?
            out_poly = (outFolder, fc + "_" + str(secteur) + "_" + string(count))

del row, rows, pochoir, count, outFolder, fc, fcList


texts:
Why are my labels placing outside my polygons?

My other labels are placed inside their polygons.

I am using ArcGIS 10.4 for Desktop.


texts: GeoPandas Overlay Intersect not intersecting some areas

texts:
I have a GeoDataFrame where I want to get all intersection between the polygons comprised in it. The cope i use for this is this:
intersection_list = []
for index in input_gdf.index.to_list():
    # first we create the intersection between the current row and a slice of the dataframe ommiting the current row
    input_without_current_gdf = input_gdf[input_gdf.index!=index]
    current_row_gdf = input_gdf[input_gdf.index == index].drop(['area','bb_id','id'],axis=1)

    intersection_gdf = gpd.overlay(input_without_current_gdf,current_row_gdf,how='intersection')

    # then we see if there are any intersection generated
    if intersection_gdf.shape[0] > 0:
            intersection_gdf['area'] = intersection_gdf['geometry'].area
            # and check if those intersection are larger than 0.0001
            found_intersections = intersection_gdf.loc[intersection_gdf['area']>=0.0001]['geometry'].to_list()
            intersection_list += found_intersections

This mostly works but in some cases no intersections are created even though there are clearly some for example here:

I visualized the overlapping area between area 417 and 301 by making the the layer semi transperent the darker gray area should be generated as an intersect but it isnt. Instead 613 is create as quasi empty area. Now if i do the following:
snapped_gdf.loc[417].geometry.intersection(snapped_gdf.loc[301].geometry)

the correct intersection is created
is it that area 277 is interfering and if more than two areas are overlapping an incorrect intersection is created?


texts:
I can't get my exported map to show my raster layers as they appear in the layout.  Specifically, my layout includes four data frames, each with its own raster showing areas stippled with colored pixels, but the exported map makes it look as though much more of the area is colored.  I've tried this in both ArcMap 10.6 and ArcGIS Pro 2.1.2 and have tried exporting in every format and at nearly ever setting (even ridiculously high dpi) but can't seem to achieve a WYSIWYG.  C and D seem okay, but A and B (much more diffuse) renders very wrong and certainly doesn't accurately portray the data.

These are simply rasters of unsuitable/suitable (0/1).  DRA is not checked.  I've tried every color combination (light/dark, dark/light), tried adjusting every setting under Image Analysis (contrast, brightness, etc.) and all the resampling methods (the default being nearest neighbor).  Of note, setting resample = majority makes for much fewer colored areas (=1) appearing on the layout, but this doesn't seem to cause any change in the exported map.  Deleting the ovr files (as suggested by @johns) made no discernible difference. It is still a poor representation of the data and what I see in the layout.  Looking at the exported B above, most think that 75-85% of the area =1, when in fact it is only 38%.


texts:
I am using the popular QSpatiaLite plugin for QGIS. Since I am creating and modifying SpatiaLite tables directly with Python (using pyspatialite), I wonder whether I can automate tasks such as converting from a Shapefile to SpatiaLite directly from the Python code (similar to the plugin's "Import OGR files").

Thus, does the plugin provide some sort of an API to do that kind of stuff? Or can I do this from within pyspatialite without using any QGIS-related API?

I know there is a QGIS Python API but I have never dealt with it, that's why I am kind of confused here.


texts:
I have connected to the correct folder on ArcCatalog.

Only the .tif files show on my ArcCatalog tree, but the .tlw files are in that folder too. How do I create a raster layer with this data?

I have tried to drag these over to ArcMap but it says the data is not georeferenced.

There are hundreds of tiles, is having them both in the same folder enough? How else do I get the program to recognise them as the georeferences for my raster .tif tiles?


texts:
How might I go about finding the following points of a MULTIPOLYGON in PostGIS?

Northeasternmost
Northwesternmost
Southeasternmost
Southwesternmost

Can't seem to find a ST_ function to accomplish this.


texts: How to calculate point from lat/long, angle, and distance?

texts:
I am trying to calculate the point from a given lat/long, a bearing, and a distance.

This is for plotting weather given a aeronautical way point.

For example:

Given a waypoint with a known lat/long (LAX airport control tower), with a bearing of North by Northwest and a distance of 45 nautical miles, where is the point?

I could treat the world as flat, but I think the resulting calculations would be inaccurate the further north the point is since longitudes get more narrow near the poles.

I am mainly interested in plotting these points in the US.  Is there a projection that would help me do this?

This is for a Python application I'm writing, and the libraries I have access to are GDAL and Proj4 for projections. So no fancy tools like ArcGIS are available to me.

Ultimately I want to combine 4 or these calculated points to form a polygon that should roughly be about 100-200 miles in width and height.

Thanks!


texts: I'm trying to convert some vector layers(format shp) to raster, to do so I used the following method raster-->conversion-->(vector to raster) and the processing works without any problem. However the result of the conversion is not satisfying at all I got just a black image. I don't know from where the problem is.


texts: Interpolating LiDAR: same steps but different results?

texts:
I'm new to eCognition, and today I tried to use LiDAR data in the program.

First, I converted LiDAR to new image layers: 1) Elevation (first return) and 2) Number of Returns (all returns). I got some pixels with 0 value so I tried to interpolated these two new layers.

First I did 'Interpolating LiDAR Number of Returns', everything worked fine, I got a new image layer called 'Number of Returns Interpolated' (where formally 0 values are updated and contain now the mean of the mean
values of its neighbors).

Then I started 'Interpolating LiDAR Elevation First' with exact the same steps:

First I carried out 'chessboard segmentation' and created a new level called '_temp Level' then I carried out 'assign class' so the Mean Elevation First > 0 will be assigned to 'Objects with value'.  However, a warning showed up and said 'Domain is empty'.  How could this happen? I mean I indeed followed the same step!

The rule set was made by my institution, if I just execute the whole rule set then there are no problem, but if I re-write rule set and execute each individual child process, the problem I described above happens. Any tips?